R：对数似然估计中生成的二元运算符的非数字参数

问题描述

我想尝试通过在 R 中使用 optim() 函数来找到最大似然估计量。我首先使用 GLM 模型对数据进行建模，以将估计值与 optim() 进行比较。

这是我的代码

 d <- read.delim("http://dnett.github.io/S510/disease.txt")

d$disease=factor(d$disease)
d$ses=factor(d$ses)
d$sector=factor(d$sector)
#str(d)

oreduced <- glm(disease~age+sector,family=binomial(link=logit),data=d)
summary(oreduced)

y<-as.matrix(d$disease)
x1<-as.matrix(d$age)
x2<-as.matrix(d$sector)

nlldbin=function(param){
  eta<-param[1]+param[2]*x1+param[3]*x2
  p<-1/(1+exp(-eta))
  -sum(y*log(p)+(1-y)*log(1-p),na.rm=TRUE)
}
MLE_estimates<-optim(c(0.1,0.1,0.1),nlldbin,hessian=TRUE)
MLE_estimates

结果是

Error in param[3] * x2 : non-numeric argument to binary operator

有人可以解决这个问题吗？谢谢。

解决方法

在使用因子变量时需要非常小心。有时，它们像字符串一样，有时像数字一样。对 as.matrix(d$disease) 的调用将 d$disease 视为一个字符串，为您提供以下无意义的输出：

> head(y)
     [,1]
[1,] "0" 
[2,] "0" 
[3,] "0" 
[4,] "0" 
[5,] "1" 
[6,] "0"

你想要的可能是这样的：

y<-as.numeric(as.character(d$disease))
x1<-as.numeric(as.character(d$age))
x2<-as.numeric(as.character(d$sector))

> head(y)
[1] 0 0 0 0 1 0

有关此主题的更多信息，我建议参阅 R Inferno 的第 8.2 节。它有点过时了，但它会有所帮助。

顺便说一句，如果我没有弄错你的意图，你可能想用这个替换倒数第二行：

MLE_estimates<-optim(c(Intercept=0.1,age=0.1,sector2=0.1),nlldbin,hessian=TRUE)

这将使输出更易于阅读。

glm mle r r