r语言 - 通过迭代加权最小二乘回归获得 beta 的 MLE

Question

我有以下数据集：

y <- c(5,8,6,2,3,1,2,4,5)
x <- c(-1,-1,-1,0,0,0,1,1,1)
d1 <- as.data.frame(cbind(y=y,x=x))

当我使用glm()将模型拟合到该数据集时，使用带有日志链接的泊松分布：

model <- glm(y~x, data=d1, family = poisson(link="log"))
summary(model)

我得到以下输出：

Coefficients:
Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.3948     0.1671   8.345   <2e-16 ***
x            -0.3038     0.2250  -1.350    0.177   

我想为迭代重新加权最小二乘回归编写一个函数，以获得相同的估计值。到目前为止，我已经能够使用身份链接来做到这一点，但不能像在 glm 中那样使用日志链接来做到这一点。

X <- cbind(1,x)
#write an interatively reweighted least squares function with log link
glmfunc.log <- function(d,betas,iterations=1)
{
X <- cbind(1,d[,"x"])
z <- as.matrix(betas[1]+betas[2]*d[,"x"]+((d[,"y"]-exp(betas[1]+betas[2]*d[,"x"]))/exp(betas[1]+betas[2]*d[,"x"])))

for(i in 1:iterations) {
W <- diag(exp(betas[1]+betas[2]*d[,"x"]))
betas <- solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z
}
return(list(betas=betas,Information=t(X)%*%W%*%X))
}
#run the function
model <- glmfunc.log(d=d1,betas=c(1,0),iterations=1000)

这给出了输出：

#print betas
model$betas
[,1]
[1,]  1.5042000
[2,] -0.6851218

有谁知道我在编写自定义函数时哪里出错了，以及如何更正它以复制glm()函数的输出

Answer 1

看起来你的"z"需要在你的循环内，因为你的"beta"每次迭代都会更新，因此你的"z"也应该基于这些值。

否则，实现看起来是正确的。