Obtenir le MLE des bêtas par des carrés repondérées moins itérative régression

Question

je jeu de données:

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)) 

I lorsque je forme un modèle à cet ensemble de données avec glm(), en utilisant une distribution de poisson avec un lien de journal:

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

je reçois la sortie suivante:

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 

Je veux écrire une fonction pour l'il régression rééquilibrée des moindres carrés qui obtiendront les mêmes estimations. Jusqu'à présent, j'ai été capable de le faire en utilisant un lien d'identité, mais pas un lien de journal, comme je le fais dans le 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) 

Ce qui donne une sortie:

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

Est-ce que quelqu'un sait où je vais mal par écrit la fonction personnalisée et comment je corriger cela pour répliquer la sortie de la fonction glm()

Answer 1

Il Apparemment, votre 'z' doit être à l'intérieur de votre boucle, car vos 'betas' sont mis à jour à chaque itération, ainsi votre 'z' devrait être basé sur ces valeurs.

L'implémentation semble correcte.