Generalised Linear Models Computer Lab 2

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 µi
1 − µ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 and 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, 2
2
)), 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 µ
∗
1 − µ∗
= β0 + β1.30 = λ
The logistic model can be rewritten equivalently (i.e. reparameterised) as
6Yi ∼ Binomial(6, µi)
log µi
1 − µi
= β0 + β1xi = β0 + β130 + β1(xi − 30)
= λ + β1(xi − 30)
This is now a Binomial GLM with explanatory variables 1 and (x − 30), and parameters λ and β1. This
reparameterised model can be fitted as follows:
chalR <- glm(damage/6~I(temp-30), family=binomial, weights=rep(6,23),data=orings)
coef(chalR)
## (Intercept) I(temp - 30)
## 5.1759798 -0.2162337
This gives λˆ = 5.176. We can obtain a profile likelihood CI for λ via:
confint(chalR)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 1.938259 8.799512
## I(temp - 30) -0.332657 -0.120179
Using the fact that µ
∗ = e
λ/(1 + e
λ
), we can transform the interval to get a 95% profile likelihood confidence
interval for µ
∗
:
cil <- c( 1.938259, 8.799512 )
exp(cil)/(1+exp(cil))
## [1] 0.8741608 0.9998492
2
Estimated variance-covariance matrix
Recall that for a model with vector parameter β = (β1, . . . , βp)
T
the approximate distribution of the maximum
likelihood estimator is β ∼ N(β, I(β)
−1
), with asymptotic variance matrix I(β)
−1 depending on the true
values of the parameters.
The estimated asymptotic variance matrix of βˆ is obtained by plugging in the MLE of β:
Var( d βˆ) = I(βˆ)
−1
.
For logistic and Poisson regression, the estimated asymptotic variance matrix can easily be obtained in R by
passing the fitted model object to the command vcov(). For example, for the Challenger model:
vcov(chal1)
## (Intercept) temp
## (Intercept) 10.865351 -0.174240974
## temp -0.174241 0.002827797
Delta method
Recall that the estimated asymptotic variance matrix can be used in the delta method. In the lecture notes
we showed that h(β) can be estimated by h(βˆ), and a 100κ% confidence interval for h(β) is given by the end
points
h(βˆ) ± z(1−κ)/2 s. e.(βˆ) = h(βˆ) ± z(1−κ)/2
q
∇h(βˆ)
T I(βˆ)−1∇h(βˆ)
For example suppose that x
∗
is the temperature at which there is a 50% probability of O-ring damage. Note
that x
∗ = h(β) = −β0/β1, and this can be estimated by h(βˆ) = −βˆ
0/βˆ
1:
betah <- coef(chal1)
betah
## (Intercept) temp
## 11.6629897 -0.2162337
xh <- -betah[1]/betah[2]
xh
## (Intercept)
## 53.93697
Note that βˆ
0 is betah[1] and βˆ
1 is betah[2], due to the fact that R indexes the elements of a vector
beginning with a 1.
To compute the 95% CI for x
∗ note that ∇h = (−1/β1, β0/β2
1
)
T and use the following
gradh <- c( -1/betah[2], betah[1]/(betah[2]ˆ2) )
gradh
## temp (Intercept)
## 4.624627 249.438380
se <- sqrt( t(gradh)%*% vcov(chal1)%*%gradh)
se
## [,1]
## [1,] 2.515673
xh + c(-1,1)*qnorm(0.975)*c(se)
## [1] 49.00635 58.86760
3
This agrees with the results in the lecture – however this time we didn’t have to do the calculations by hand.
Note that in the computation of the standard error we used the following commands:
• t gives the transpose of a vector/matrix
• %*% corresponds to matrix multiplication.
Exercises 1
1 For the Challenger example:
a) Use the output of vcov to compute the standard error of βˆ
0 and βˆ
1, and compare the answer to the
output of summary.
b) Calculate 95% Wald and profile likelihood intervals for the intercept parameter.
c) Compare the width and symmetry of the Wald and profile likelihood intervals for both parameters.
What do you notice? Which do you expect to be more accurate?
d) Conduct Wald and profile likelihood tests of the hypothesis H0 : βT emp = −0.33 vs H1 : βT emp =6 −0.33
at a 5% significance level.
2 Data have been collected in each of n = 31 areas. In the ith area, we have:
• xi
: the number of Eucalypt trees (named Eucs in R)
• yi
: the number of noisy miner birds (also called the abundance, named Minerab in R).
The data are given in the file nminer2.csv. Using R:
a) Fit a Poisson regression model with abundance as the response variable and the number of Eucalypt
trees as the explanatory variable, and write down the fitted model.
b) Using three different methods, find a 95% confidence interval for φ, the expected abundance in an area
with 6 Eucalypt trees. Comment on which method you think is likely to be best.
c) Using three different methods, find a 95% confidence interval for λ, the percentage increase in expected
abundance associated with one additional Eucalypt tree. Comment on which method is likely to be
best.
Fisher scoring
Recall from Lecture 6 that for a statistical model with parameter vector θ the maximum likelihood estimates
θˆ can be approximated numerically via an iterative procedure known as Fisher scoring, which proceeds as
follows:
1. Set initial approximate values for the estimates, θˆ
0
2. For k = 0, 1, 2, . . . until convergence set
θˆ
k+1
= θˆ
k
+ I(θˆ
k
)
−1u(θˆ
k
)
Recall further from Lectures 5-6 that the score vector and Fisher information matrix for simple logistic
regression are
u(β) =  P
P
i mi(yi − µi)
i
ximi(yi − µi)

, I(β) =  P
i Wi
P
P
i Wixi
i Wixi
P
i Wix
2
i

,
where Wi = miµi(1 − µi), and that these can be calculated using the R functions u and FIM given in the file
Lab2-functions.R.
4
Exercises 2
3 Using these functions with the Challenger data:
a) calculate the estimated asymptotic variance matrix of βˆ and compare your result to the output of vcov.
b) without using the glm function, carry out 10 iterations of Fisher scoring starting with initial vector
βˆ
0
= (0, 0)
4 Recall from Example Sheet 3 that for simple Poisson regression the score vector and information matrix
are
u(β) = " P
i
(yi − µi)
P
i
xi(yi − µi)
#
, I(β) =  P
i µi
P
P
i µixi
i µixi
P
i µix
2
i

where µi = e
β0+β1xi
.
a) by adapting the given functions u and FIM, write functions to calculate the score vector and Fisher
information matrix for simple Poisson regression
For the noisy miner model in Question 2:
b) use your function to calculate the estimated asymptotic variance matrix of βˆ, and compare the result
to vcov
c) use your functions to carry out 15 iterations of Fisher scoring, and verify that the results coincide with
those from glm.
5