0

我有以下数据集：

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()函数的输出

4

1 回答 1

1

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

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

于 2017-09-20T02:10:16.020 回答