0
$\begingroup$

The difference of Standard Error between glm(y~x, family=poisson(link=identity)) and optim() in R

Following the above question...

And actually, I have an another question now. I could cross validation for glm(poisson, log) to understand Jacobian matrix from log-likelihood by the following program . But I confuse... I can't do with “link=identity”. I wonder how should I change the program for glm(poisson, identity)... Someone, please give me some advice.

x<-c(1,2,3,4)
y<-c(2,3,5,4)
f1<-expression(sum(y-exp(a+b*x)))#simultaneous equations1
f2<-expression(sum(x*(y-exp(a+b*x)))) #simultaneous equations2
g11<-expression(-sum(exp(a+b*x)))#first‐order differentiation for f1 by a
g12<-expression(-sum(x*exp(a+b*x)))#first‐order differentiation for f1 by b
g21<-expression(-sum(x*exp(a+b*x)))#first‐order differentiation for f2 by a
g22<-expression(-sum(x**2*exp(a+b*x)))#first‐order differentiation for f2 by b
a<-1#initial value
b<-0.5#initial value
for(i in 1:10){
  m<-matrix(c(a,b),2,1)#a,b to matrix
  f<-matrix(c(eval(f1),eval(f2)),2,1)#f1,f2 to matrix
  j<-matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)#Jacobian matrix
  m<-m-solve(j)%*%f
  a<-m[1]
  b<-m[2]
  print(sprintf("[%d] (a,b)=(%f,%f)",i,a,b))
}
#cut
#[1] "[7] (a,b)=(0.639265,0.232040)"
#[1] "[8] (a,b)=(0.639265,0.232040)"
#[1] "[9] (a,b)=(0.639265,0.232040)"
#[1] "[10] (a,b)=(0.639265,0.232040)"
(SEs<-sqrt(abs(diag(solve(j)))))
#[1] 0.7317073 0.2445158

summary(glm(y~x,family=poisson(link=log)))
#            Estimate Std. Error z value Pr(>|z|)
#(Intercept)   0.6393     0.7317   0.874    0.382
#x             0.2320     0.2445   0.949    0.343
$\endgroup$
0
$\begingroup$

I got the value. But it still didn’t have validation of SE. I don’t know why…

x<-c(1,2,3,4)
y<-c(2,3,5,4)
f<-expression((exp(-a-b)*(a+b)^2/2)*(exp(-a-2*b)*(a+2*b)^3/6)*(exp(-a-3*b)*(a+3*b)^5/120)*(exp(-a-4*b)*(a+4*b)^4/24))
f1<-D(f,"a")
f2<-D(f,"b")
g11<-D(f1,"a")
g12<-D(f1,"b")
g21<-D(f2,"a")
g22<-D(f2,"b")
a<-1.2
b<-0.7
for(i in 1:10){
  m<-matrix(c(a,b),2,1)
  f<-matrix(c(eval(f1),eval(f2)),2,1)
  j<-matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)
  m<-m-solve(j)%*%f 
  a<-m[1]
  b<-m[2] 
  print(sprintf("[%d] (a,b)=(%f,%f)", i, a, b))
}
#cut
#[1] "[7] (a,b)=(1.278347,0.888661)"
#[1] "[8] (a,b)=(1.278347,0.888661)"
#[1] "[9] (a,b)=(1.278347,0.888661)"
#[1] "[10] (a,b)=(1.278347,0.888661)"
(SEs<-sqrt(abs(diag(solve(j)))))
#[1] 50.43801 20.72464

model<-glm(y~x, family=poisson(link="identity"))
#            Estimate Std. Error z value Pr(>|z|)
#(Intercept)   1.2784     1.9766   0.647    0.518
#x             0.8887     0.8141   1.092    0.275
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.