は、私は、データセットを持って反復再重み付き最小二乗回帰

Question

を通じてベータのMLEを取得します。

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'はあなたのループの中にある必要があります。あなたの 'ベータ'が各繰り返しを更新するので、それらの値に基づいてあなたの 'z'を更新する必要があります。

そうでなければ、実装は正しく見えます。