通过迭代重加权最小二乘回归得到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() 将模型拟合到此数据集时，我使用具有对数 link 的泊松分布：

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

我想为迭代重新加权最小二乘回归编写一个函数，以获得相同的估计值。到目前为止，我已经能够使用身份 link 来执行此操作，但不能像我在 glm.

中那样使用日志 link

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' 需要在您的循环中，因为您的 'betas' 每次迭代都会更新，因此您的 'z' 也应该如此，因为它基于这些值。

除此之外，实施看起来是正确的。

通过迭代重加权最小二乘回归得到beta的MLE

Obtain the MLE of betas through iterative re-weighted least squares regression

r

regression

glm

log-likelihood