MATH38172 Generalised Linear Models
Computer Lab 2
Inference for logistic and Poisson regression
Wald inference
For logistic and Poisson regression, it is easy to obtain Wald tests and intervals from the output of the R
command summary(), passing the fitted model object as an argument. As an example consider the Challenger
data.
Recall that the fitted model is
6Yi ～ Binomial(6, μi)
log μi1 μi = β0 + β1xi
where Yi is the proportion of O-rings damaged on the ith launch, μi is the probability of an individual O-ring
being damaged, and β0, β1 are unknown parameters.
Now consider the following R code (assuming you have already loaded the orings data):
chal1 <- glm( damage/6 ~ temp, weights=rep(6,nrow(orings)), family=binomial, data=orings)
summary(chal1)
The following table is given as part of the output:
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 11.66299 3.29626 3.538 0.000403 ***
temp -0.21623 0.05318 -4.066 4.78e-05 ***
—
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The output shows for each parameter
the maximum likelihood estimate β j
the standard error, s. e.(β j) = [I(β ) 1]jj
the test statistic (z-value) and p-value for a Wald test of the hypothesis
H0 : βj = 0 vs H1 : βj 6= 0 .
The MLEs an
d standard errors can easily be used to conduct other Wald inferences about the parameter βj .
For example, the end-points of a 95% Wald confidence interval for βTemp can be found using the formula
β j ± z0.025 s. e.(β j) as follows
-0.21623 + c(-1,1)*qnorm(0.975)*0.05318
## [1] -0.3204609 -0.1119991
Recall that the qnorm function gives quantiles of a normal distribution. For example:
1
(i) to compute the 0.95 quantile of the standard normal (i.e. N(0, 1)) distribution use qnorm(0.95).
(ii) to compute the 0.95 quantile of a normal distribution with mean 3 and variance 4 (i.e. N(3, 22)), use
qnorm(0.95,mean=3,sd=2).
Profile likelihood intervals
Profile likelihood intervals for each parameter βj can be easily obtained by passing the fitted model object as
an argument to the confint() command, e.g for the Challenger model:
confint(chal1, level=0.9)
## Waiting for profiling to be done…
## 5 % 95 %
## (Intercept) 6.5165681 17.5027878
## temp -0.3120442 -0.1347199
Note that we can sometimes find profile likelihood intervals for functions of the parameter by fitting a GLM
with different explanatory variables. For example, suppose we want to find a 95% profile likelihood confidence
interval for the probability of O-ring damage at 30 F , which we denote μ . Note that
log μ