CODEWITHSUNDEEP
r
John Schroeder
John Schroeder

Reputation: 85

Standard error in glm output

I am using r glm to model Poisson data binned by year. So I have x[i] counts with T[i] exposure in each year, i. The r glm with poisson family log link output produces model coefficients a, b for y = a + bx.

What I need is the standard error of (a + bx) not the standard error of a or the standard error of b. The documentation describing a solution I am trying to implement says this should be calculated by the software because it is not straightforward to calculate from the parameters for a and b. Perhaps SAS does the calc, but I am not recognizing it in R.

I am working working through section 7.2.4.5 of the Handbook of Parameter Estimation (NUREG/CR-6823, a public document) and looking at eq 7.2. Also not a statistician so I am finding this is very hard to follow.

The game here is to find the 90 percent simultaneous confidence interval on the model output, not the confidence interval at each year, i.

Adding this here so I can show some code. The first answer below appears to get me pretty close. A statistician here put together the following function to construct the confidence bounds. This appears to work.

# trend line simultaneous confidence intervals
# according to HOPE 7.2.4.5
HOPE = function(x, ...){
t = data$T
mle<-predict(model, newdata=data.frame(x=data$x), type="response") 
se = as.data.frame(predict(model, newdata=data.frame(x=data$x), type="link", se.fit=TRUE))[,2]
chi = qchisq(.90, df=n-1)
upper = (mle + (chi * se))/t
lower = (mle - (chi * se))/t
return(as.data.frame(cbind(mle, t, upper, lower)))

}

Upvotes: 1

Views: 2399

Answers (1)

Julian Wittische
Julian Wittische

Reputation: 1237

I think you need to provide the argument se.fit=TRUE when you create the prediction from the model:

hotmod<-glm(...)
predz<-predict(hotmod, ..., se.fit=TRUE)

Then you should be able to find the estimated standard errors using:

predz$se.fit

Now if you want to do it by hand on this software, it should not be as hard as you suggest:

covmat<-vcov(hotmod)
coeffs<-coef(hotmod)

Then I think the standard error should simply be:

sqrt(t(coeffs) %*% covmat %*% coeffs)

The operator %*% can be used for matrix multiplication in this software.

Upvotes: 1

Related Questions