This is an HTML version of an attachment to the Freedom of Information request 'Statistics: Past Exam Papers'.

Tuesday, 1st May 2018
2.00 pm - 3.30 pm
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics – STATS 4024
Stochastic Processes
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
This paper consists of 4 pages and 3 questions, with the following marks:
Question
Marks
1
25
2
25
3
10
Total
60
Candidates should attempt ALL questions. No statistical tables or probability
formula sheets are permitted.
1. Consider a homogeneous Markov chain in discrete time defined over a finite state
space. The chain is defined by the stochastic matrix A with elements Aik ∈ R, which
define the probability of a transition from state i to state k in one time step. Assume
that at time 0, the system is in state i.
(a) Which conditions do the matrix elements Aik have to satisfy for A to be a valid
stochastic matrix?
[1 MARK]
CONTINUED OVERLEAF/
1

(b) Let g(m) denote the probability that the Markov chain returns to state i in m
i
steps. Express g(2) in terms of A
i
ik .
[1 MARK]
(c) Let f (m) denote the probability of a first return to state i in m steps. Use the law
i
of total probability to derive an expression for f (2) in terms of g(2), f (1) and A
i
i
i
ik .
[2 MARKS]
(d) Let fi denote the probability of a return to state i. Explain why fi = P∞
g(m)
m=1
i
is wrong and why you have to use fi = P∞
f (m) instead.
[2 MARKS]
m=1
i
(e) Consider a homogeneous discrete-time Markov chain defined by the following
stochastic matrix:

0
1
1 
2
2
A =
1
1
0

2
2

0
0
1
i. Draw the state transition diagram.
[1 MARK]
ii. Use this state transition diagram and the result from part (d) of the question
to decide which states are transient and which states are persistent. Justify
your answer!
[6 MARKS]
iii. Find the stationary distribution of this Markov chain. Justify your answer!
[2 MARKS]
iv. Does this Markov chain have a unique limiting distribution? Find the limiting
distribution if it exists. Justify your answer!
[4 MARKS]
(f) Let xt denote the state of the Markov chain at time t. Show that if (x1, x2, . . . , xm)
is a homogeneous Markov chain that has converged to its limiting distribution,
then the reverse process running backwards in time, (xm, xm−1, . . . , x1), is also a
homogeneous Markov chain.
[4 MARKS]
(g) Assume you want to sample from a discrete probability distribution p(i) = φ(i) ,
Z
where φ(i) is easy to compute, but the normalisation constant Z is intractable.
Use the result from the previous part of the question to discuss how that can be
done.
[2 MARKS]
2. (a) Let N (t) ∈ Z be a discrete time-varying random variable representing the number
of events that have occurred in interval [0, t], where t denotes time and N (0) = 0.
Define pn(t) = P [N (t) = n]. Recall that the probability generating function for a
CONTINUED OVERLEAF/
2

discrete stochastic process is defined as:

X
G(s, t) =
pn(t)sn
n=0
Assume that, for λ > 0, the differential equation for the probability generating
function is given by:
∂G(s, t)
= λ(s − 1)G(s, t)
∂t
i. Solve the differential equation under the boundary condition N (0) = 0.
[3 MARKS]
ii. Use this result to prove that the stochastic process defined by the differential
equation above is the Poisson process with rate λ.
[3 MARKS]
iii. Show that the inter-arrival times, i.e. the time intervals between successive
events, have an exponential distribution.
[3 MARKS]
iv. Derive an expression for the mean or expected inter-arrival time.
[3 MARKS]
v. Patients suffering from flu-like symptoms are admitted to the A&E depart-
ment of a hospital at a rate of 10 per day. Assume that admissions form a
Poisson process, and that doctors at the department work in 12-hour shifts.
From the beginning of a shift, what is the expected time a doctor has to wait
until the first call, and what is the probability that exactly 6 patients are
admitted during a shift?
[3 MARKS]
(b) Let T be a nonnegative random variable which is the time to failure for a given
device. The distribution function of T is F (t) = P (T ≤ t). The probability
density function is denoted by f (t). The reliability function R(t) = P (T ≥ t) is
the probability that the device is still operational after time T . The failure rate
function or hazard function is defined as
P (t < T ≤ t + δt|T > t)
r(t) =
lim
δt→0
δt
i. Express dR in terms of f (t).
[1 MARK]
dt
ii. Give an interpretation of the hazard function r(t). What does it measure?
[1 MARK]
iii. Show that
f (t)
r(t) = R(t)
[3 MARKS]
CONTINUED OVERLEAF/
3

iv. Show that

Z
t

R(t) = exp −
r(u)du
0
[3 MARKS]
v. Whilst in use, a printer is observed to have a hazard function of r(t) = 2λt
per hour where λ = 0.0001 hours−2. What is the probability that the printer
is still operational after 100 hours?
[2 MARKS]
3. In the lectures we discussed a queue with a single server. Let the random variable N (t)
denote the number of individuals in the queue at time t (including the person being
served). New customers arrive independently, and these arrivals form a Poisson process
with intensity λ > 0. Service times have an exponential distribution with parameter
µ > 0, that is the distribution with probability density function f (t) = µ exp(−µt), t ≥
0. Now consider a baulked queue where no more than m people (including the person
being served) are allowed to form a queue. If there are m people in the queue, then
any further arrivals are turned away. It can be shown that the probability distribution
for this queue is defined by the following differential equations:
dpm(t) = λpm−1(t) − µpm(t)
(1)
dt
dpn(t) = λpn−1(t) + µpn+1(t) − [λ + µ]pn(t); 1 ≤ n ≤ (m − 1)
(2)
dt
dp0(t) = µp1(t) − λp0(t)
(3)
dt
(a) How do you get the limiting distribution limt→∞ pn(t) = pn, where pn is time
invariant, from these differential equations?
[1 MARK]
(b) Define ρ = λ . Show that p
µ
n = 1 and pn = ρn are solutions to equation (2).
[2 MARKS]
(c) Assuming that ρ 6= 1, find the general solution for pn, 1 ≤ n ≤ m. Hint: Make
use of the mathematical techniques that you have learned in the first lectures on
finite difference equations and treat equations (1) and (3) as boundary conditions.
[6 MARKS].
(d) Assume that the service rate µ is twice as high as the arrival rate λ and that the
maximum permissible length of the queue is m = 10. What is the probability that
n = 5 people are in the queue?
[1 MARK].
END OF QUESTION PAPER.
4

Thursday, 2nd May 2019
09.30 a.m. - 11.00 a.m.
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics – STATS 4024
Stochastic Processes
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
This paper consists of 5 pages and 5 questions, with 12 marks per question.
Candidates should attempt ALL questions.
1. Suppose that deuce (equal score) has been reached in a tennis game between Adrian
and Marian. The player winning the next point has advantage. On the following
point, the player having advantage either wins the game or the game returns to deuce.
Marian has probability p = 1/3 of winning the next point and Adrian has probability
q = 2/3. Assume that the game can be modelled with a homogeneous Markov chain,
with the following five states: k = 0: match point Adrian, k = 1: advantage Adrian,
k = 2: deuce, k = 3 advantage Marian, k = 4: match point Marian. What is the
probability that Adrian will win the game (or equivalently, that Marian will lose the
game)? Break the problem up into the following subproblems:
(a) Let Ck denote the event that Adrian will win (or, equivalently, that Marian will
lose) the game when the game is in state k. Apply the law of total probability to
prove the following homogeneous difference equation for uk = P (Ck):
uk = uk+1p + uk−1q
[2 MARKS]
CONTINUED OVERLEAF/
1

(b) What are the boundary conditions for the difference equation in Part (a)?
[2 MARKS]
(c) Find the general solution of the homogeneous difference equation in Part (a).
[4 MARKS]
(d) Find the specific solution of the homogeneous difference equation subject to the
boundary conditions from Part (b).
[3 MARKS]
(e) What is the probability that Adrian will win the game?
[1 MARK]
2. A simple model for molecular diffusion is an unrestricted symmetric random walk on
a straight line starting at the origin: X0 = 0, where Xn denotes the position at step
n. In each step, n, the molecule either moves one unit to the right, with probability
P (Xn+1 = m + 1|Xn = m) = p = 0.5, or one unit to the left, with probablity
P (Xn+1 = m − 1|Xn = m) = q = 1 − p = 0.5. Let An denote the event that the
molecule is at the origin at step n: Xn = 0; let Bn denote the event that the first
return to the origin occurs in the n-th step. Use the following abbreviations for the
probabilities of An and Bn: vn = P (An) and fn = P (Bn).
(a) Define V (s), the probability generating function for vn, and F (s), the probability
generating function for fn. (Note: Do not copy the specific expressions given in
the next part of the question. What is needed here is a general definition.)
[2 MARKS]

(b) It can be shown that V (s) =
1

, and F (s) = 1 −
1 − s2. Find lim V (s) and
1−s2
s→1
lim F (s).
[1 MARK]
s→1
(c) What does this imply for vn and fn? Are vn and fn proper probability distribu-
tions? Choose the correct option from the following table and write it down in
your script book (only 1 option is correct):
Option
1) Both fn and vn are proper probability distributions.
2) Only fn is a proper probability distribution, but vn is not.
3) Only vn is a proper probability distribution, but fn is not.
4) Neither vn nor fn are proper probability distributions.
5) To become proper probability distributions, vn and fn have to be
normalised by dividing them by lim V (s) and lim F (s), respectively.
s→1
s→1
6) vnV (s) and fnF (s) are proper probability distributions.
[2 MARKS]
CONTINUED OVERLEAF/
2

(d) Select the option that best explains the reason for the difference between
lim V (s) and lim F (s).
s→1
s→1
There is no first return to the origin at step n = 0. For that
1)
reason the events {Bn} are not exhaustive and
fn is not a proper probability distribution.
The events {An} are exhaustive and mutually exclusive
2)
and therefore define a proper partition of the event space,
but the events {Bn} do not.
The events {Bn} are exhaustive and mutually exclusive
3)
and therefore define a proper partition of the event space,
but the events {An} do not.
The events Bn are constrained to include a 1st return to the origin
4)
at step n. For that reason, they do not cover the whole event space, and
fn is not a proper probability distribution. To get a proper probability
distribution we need to multiply fn by lim F (s).
s→∞
The events Bn are constrained to include a 1st return to the origin
5)
at step n, whereas there is an infinite number of realisations for the
events An. For that reason, lim F (s) is well-defined, but
s→1
for V (s), the limit has to be taken at infinity: lim V (s).
s→∞
[3 MARKS]
(e) Use the expression of F (s) to obtain the probability of a first return to the origin
after 2 steps and after 3 steps. Justify your answer!
[4 MARKS]
3. In a certain town, there are four entertainment venues. Both Marian and Adrian visit
one of these venues every weekend, independently of each other. Each of them visits
the venue of the week before with probability 0.4 and chooses otherwise at random
one of the other three venues. Treat this problem as a two-state Markov chain, where
state 0 means that Marian and Adrian are at different venues and state 1 means that
they are at the same venue.
(a) Apply the law of total probability to obtain the one-step transition probability P01,
by deriving the probabilities of three mutually exclusive and jointly exhaustive
events: that Adrian does not change venue and Marian goes to Adrian’s venue,
that Marian does not change venue and Adrian goes to Marian’s venue, and that
both Marian and Adrian change venue and go to the same venue. Apply a similar
method to obtain the probability P11, by deriving the probabilities of two mutually
exclusive and jointly exhaustive events: that neither Marian nor Adrian change
venue, and that both Marian and Adrian change venue and go to the same venue.
[5 MARKS]
CONTINUED OVERLEAF/
3

(b) Find the transition matrix of the Markov chain.
[1 MARK]
(c) Is the Markov chain ergodic? What does that imply for the limiting distribution?
[2 MARKS]
(d) What is the long-run (i.e. limiting) fraction of weekends that Marian and Adrian
visit the same venue?
[2 MARKS]
(e) What is the limiting probability that Marian and Adrian visit the same venue two
weekends in a row?
[2 MARKS]
4. (a) Consider a time-homogeneous reversible Markov chain defined over discrete states
i, and let T denote its stochastic transition matrix. Show that the stationary
distribution P ∗(i) satisfies the following equation:
P ∗(i)Tik = P ∗(k)Tki
[3 MARKS]
(b) Assume you want to sample from a discrete probability distribution P (i) = φ(i)
Z
defined over the discrete state space, where φ(i) is easy to compute, but the
normalisation constant Z = P φ(k) is intractable. Assume that the stochastic
k
matrix T is ergodic. Use the result from the previous part of the question to
describe a method for approximately sampling from P (i) = φ(i) .
[3 MARKS]
Z
(c) Consider a mouse that is trapped in a closed maze. The maze has 15 rooms. The
rooms are positioned according to the following three-by-five array:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Doors connect any two neighbouring rooms. For example, room 1 is connected
with rooms 2 and 6, room 2 with rooms 1, 3 and 7, room 7 with rooms 2, 6, 8
and 12, and so on. The mouse moves randomly from room to room, where at each
time one of the available doors is chosen with equal probability.
i. It can be shown that there exists a positive integer n > 0 such that Tn is
a strictly positive matrix, meaning that all its elements are strictly positive.
What does that imply for the limiting distribution?
[1 MARK]
CONTINUED OVERLEAF/
4

ii. Use the equation from Part (a) to find the stationary (= limiting) distribution.
Hint: Make use of the symmetry of the maze. There are four corner rooms
(1, 5, 11, and 15), each of which has two neighbouring rooms; eight outer
rooms (2, 3, 4, 6, 10, 12, 13, 14), each of which has three neighbouring rooms;
and three inner rooms (7, 8, 9), each of which has four neighbouring rooms.
[5 MARKS]
5. Consider a queue with a single server. Let the random variable N (t) denote the
number of individuals in the queue at time t (including the person being served).
New customers arrive independently, and these arrivals form a Poisson process with
intensity λ > 0. Service times have an exponential distribution with parameter µ > 0,
that is the distribution with probability density function f (t) = µ exp(−µt), t ≥ 0. It
can be shown that the evolution of the probability pn(t) = P [N (t) = n] is described
by the following differential equations:
dpn(t) = λpn−1(t) + µpn+1(t) − [λ + µ]pn(t); n ≥ 1
dt
dp0(t) = µp1(t) − λp0(t)
dt
(a) How do you get the limiting distribution limt→∞ pn(t) = pn, where pn is time
invariant, from the differential equations?
[2 MARKS]
(b) In the lectures we showed that the limiting distribution is given by
 λ n
pn =
p0
µ
Find the value of p0 for the cases (i) λ < µ and (ii) λ ≥ µ. Interpret the result.
[5 MARKS]
(c) In the lectures we derived the following expression for the expected length of the
queue including the person being served:

X
λ
E(N ) =
npn = µ − λ
n=1
Use this expression to find the expected length of the queue excluding the person
being served:

X
E(M ) =
(n − 1)pn = ?
n=2
[5 MARKS]
END OF QUESTION PAPER.
5

Wedneday, 25th April 2018
09.30 am - 11.00 am
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics - Advanced Bayesian Methods
This paper consists of 5 pages and contains 3 questions.
Candidates should attempt all questions.
Question 1
20 marks
Question 2
20 marks
Question 3
20 marks
Total
60 marks
The following material is made available to you:
Statistical tables∗
Probability formula sheet
Discrete univariate distributions
Distribution
Probability mass function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.m.f.) p(x)
function M(t)
 
Binomial1,2
n
n ∈ N
θx(1

nθ(1
Bi
− θ)n−x
x ∈ {0, 1, . . . , n}
− θ)
(1 − θ + θ exp(t))n
No. of successes in n trials
(n, θ)
x
0 < θ < 1
θ – probability of success
Geometric
1
θ
(1 − θ) exp(t)
No. of trials until (and including)
first failure
Geo
θx−1(1 − θ)
x ∈ N
0 < θ < 1
(θ)
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success


N
No. of type I objects in a sample of
Hypergeometric
M
N −M
− n
x
n−x
N, n ∈ N
nθ(1 − θ) ·
—3
size n, drawn without replacement
HyGe

x ∈ {max{0, n−(N −M)},

N
(n, N, M )
N
− 1
. . . , min{n, M}}
M ∈ {0, . . . , N}
from a population of size N, con-
n
(with θ = M )
N
taining M type I objects.




No. of trials until (and including)
Negative Binomial
x − 1
k
k ∈ N
k

(1 − θ) exp(t)
kth failure
NeBi
θx−k(1 − θ)k
x ∈ {k, k + 1, . . .}
(k, θ)
k − 1
0 < θ < 1
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success
NeBi(1, θ) ≡ Geo(θ)
Poisson
λx
Poi
exp(−λ)
x ∈ N
(λ)
x!
0
λ > 0
λ
λ
exp(λ(exp(t) − 1))
Statistical Tables
1 Bi(n, θ) can be approximated by Poi(nθ), if n large, θ small and nθ moderate.
2 Bi(n, θ) can be approximated by N(nθ, nθ(1 − θ)), if n large and θ not too close to 0 or 1.
3 No simple closed form expression exists.
Continuous univariate distributions
Distribution
Probability density function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.d.f.) f(x)
function M(t)
Beta
xα1−1(1 − x)α2−1
α
X1 ∼ Ga(α1, θ)
1 > 0
α1
α1α2
—3
X
Be
0

≤ x ≤ 1
2 ∼ Ga(α2, θ) independent
1, α2)
B(α1, α2)
α2 > 0
α1 + α2
(α1 + α2)2(α1 + α2 + 1)

X1
∼ Be(α
X
1, α2)
1+X2
Cauchy
1
—4
—4
—4
Ca
Ca


x ∈ R
η ∈ R
(0, 1) ≡ t(1)
(η, γ)
πγ 1 + (x−η)2
γ > 0
γ2

Chi-Squared
ν
x 2 −1 exp − x
1
2

x > 0
ν ∈ N
ν

Xi ∼ N(0, 1) independent
χ2(ν)
2 ν2 Γ ν
(1 − 2t) ν2
⇒ Pν
X2
2
i=1
i ∼ χ2(ν)
Exponential
1
1
1
Expo
θ exp(−θx)
x > 0
θ > 0
(θ)
θ
θ2
1 − tθ
ν1
ν2
ν
F
ν
2
2 ν2
ν 2
2
1
2 (ν1 + ν2 − 2)
X1 ∼ χ2(ν1)
1
ν2
x 2 −1
—4
X
F

x > 0
ν
ν
2 ∼ χ2(ν2) independent

ν
1, ν2 ∈ N
2 − 2
ν1(ν2 − 2)2(ν2 − 4)
1, ν2)
B ν1 , ν2
1+ν2
2
2
(ν1x + ν2) 2
(for ν
⇒ X1/ν1 ∼ F(ν
2 > 2)
(for ν
1, ν2)
2 > 4)
X2/ν2
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
∗Note for invigilators: will be delivered to the exam venue by the School
CONTINUED OVERLEAF/
1

1. A computer administrator has collected data on web server response times over a
period of 90 days from September to December, and your job is to determine if there
are specific periods of high server activity over a semester, which could lead to higher
response times during those periods. The data set is given by y = (y1, . . . , y90), where
you can assume that each response time yi, (i = 1, . . . , 90) comes from a mixture of
two Gamma distributions, with the pdf
f (yi|α1, θ1, α0, θ0, π) = πf1(yi|α1, θ1) + (1 − π)f0(yi|α0, θ0),
where f1 and f0 denote the pdfs of Gamma distributions with parameters (α1, θ1) and
(α0, θ0), respectively, corresponding to periods of high and normal server activity.
(a) Write down a mathematical expression for the likelihood p(y1, . . . , y90|α1, θ1, α0, θ0, π),
assuming that all the measured times are independent.
[3 MARKS]
(b) For i = 1, . . . , 90, assume that Zi is a latent (unknown) variable which equals 1 if
the i-th measured time corresponds to a period of high server activity, and 0 if not
(that is, Zi ∼ Bernoulli(π)). Assuming independent Gamma(β1, β2) priors for both
θ1 and θ0, independent Uniform(10, 100) priors for α1 and α0, and a Beta (a, b)
prior for π, write down an expression proportional to the joint posterior distribu-
tion of all parameters (α1, θ1, α0, θ0, π) and the latent variables, z = (z1, . . . , z90).
[6 MARKS]
(c) Next, you plan to construct a Gibbs sampler to sample parameters from their
joint posterior distribution. At step t of the sampler, suppose you have sampled
values of α(t), α(t), θ(t), π(t), and z(t) = (z(t), . . . , z(t)). Write down the expression
0
1
0
1
90
for the full conditional distribution you would need to sample the t-th value of θ1,
conditional on all the other parameters and latent variables, and identifying it by
name.
[6 MARKS]
(d) Now, using a Gamma(c, d) proposal density for α0, write down a Metropolis-
Hastings sampling step (with an expression for the acceptance ratio) to sample
the t + 1-th iterate of α0, conditional on the values of all other parameters and
latent variables at that stage.
[5 MARKS]
2. (a) The editor of a college newspaper wants to improve its quality of printing and
enlists your help in investigating the average rate λ of misprints per sheet. A
trainee journalist has randomly selected 20 newsprint sheets from past newspapers
and counted the total number of errors Y on all of them. You may assume that
Y ∼ Poisson(20λ), and the observed number of misprints is y = 44. You are asked
to test the hypothesis:
H0 : λ ≤ 2
against H1 : λ > 2.
CONTINUED OVERLEAF/
2

It is believed that on average one misprint occurs on every 5 sheets, and you
thus choose a Gamma(0.10, 0.50) prior on the Poisson rate parameter λ. For this
problem, you may also make use of the following output from R:
> pgamma(2,.1,.5)
[1] 0.9758727
> pgamma(2,44.1,20.5)
[1] 0.3344346
i. Write down the mathematical expressions, and the numeric values, of the
prior and posterior odds in favour of H1 versus H0. Then, calculate the Bayes
Factor and state your conclusions.
[8 MARKS]
ii. Suppose you had chosen your prior as λ ∼ Gamma(0, 0) instead, because
you had doubts about the hyperparameter values stated above. Would it be
possible to compute the Bayes factor in that case? Explain why or why not.
[2 MARKS]
iii. The editor wants to review the print quality in a month’s time, at the end of
which the trainee will be asked again to take another random sample of 20
newsprint sheets and count the total number of errors, Z, in all of them. Write
down a mathematical expression for the density of the posterior predictive
distribution of Z|Y .
[4 MARKS]
(b) Suppose X is a random variable with a pdf given by
f (x|β) = 2β−2xe−(x/β)2,
x > 0.
where β > 0. Write down the steps of a Monte Carlo method to generate 100
samples from f , when the only density available to you to sample from is a U(0, 1)
(i.e. Uniform) distribution.
[6 MARKS]
3. The High School and Beyond survey conducted by the US Department of Educa-
tion collected data on students from approximately 1,100 high schools in the United
States, recording student scores on a certain achievement test along with a set of var-
ious ethnic, demographic and social characteristics of the students, their families, and
the regions in which schools were located. The variable score was the score on the
achievement tests given to high school seniors in the sample. The rest of the variables
are: unemp: county unemployment rate; wage: state hourly wage in manufacturing;
distance: distance from 4-year college (in 10 miles); tuition: average state 4-year
college tuition (in 1000 USD); education: number of years of education.
Overleaf, you are given sections of some R code and output from a Bayesian regression
analysis of a subset of the data, followed by some questions.
CONTINUED OVERLEAF/
3

> V <- 1/1000
> X <- as.matrix(collegedata[,c("unemp","wage","distance","tuition","education")])
> reg1 <- MCMCregress(score ~ X, data = collegedata, B0 = V, marginal.likelihood="Laplace" )
> effectiveSize(reg1)
(Intercept)
Xunemp
Xwage
Xdistance
Xtuition
Xeducation
sigma2
10000.00
10000.00
10000.00
10000.00
10000.00
10790.25
10000.00
>
> X <- as.matrix(collegedata[,c("unemp","wage","tuition","education")])
> reg2 <- MCMCregress(score ~ X, data = collegedata, B0 = V, marginal.likelihood="Laplace" )
> effectiveSize(reg2)
(Intercept)
Xunemp
Xwage
Xtuition
Xeducation
sigma2
10000.000
10000.000
9698.952
9708.361
10000.000
9534.389
>
> X <- as.matrix(collegedata[,c("wage","tuition","education")])
> reg3<- MCMCregress(score ~ X, data = collegedata, B0 = V, marginal.likelihood="Laplace" )
> effectiveSize(reg3)
(Intercept)
Xwage
Xtuition
Xeducation
sigma2
10000.000
10000.000
10578.708
10000.000
9620.113
> summary(reg3)
...
Sample size per chain = 10000
...
> BF <- BayesFactor(reg1, reg2, reg3)
...
reg1 :
log marginal likelihood =
-16379.98
reg2 :
log marginal likelihood =
-16373.57
reg3 :
log marginal likelihood =
-16377.28
> bvs1 <- Bvs(formula="score ~.", data=collegedata, prior.betas = "gZellner", n.keep=5)
> bvs1$modelsprob
unemp wage distance tuition education
prob
1
*
*
*
* 9.304640e-01
2
*
*
*
*
* 6.844099e-02
3
*
*
* 1.025943e-03
4
*
*
*
* 6.906445e-05
5
*
*
* 6.170775e-09
> summary(bvs1)
Call:
Bvs(formula = "score ~.", data = collegedata, prior.betas = "gZellner", n.keep = 5)
Inclusion Probabilities:
Incl.prob. HPM MPM
unemp
0.9989
*
*
wage
1
*
*
distance
0.0685
tuition
1
*
*
education
1
*
*
CONTINUED OVERLEAF/
4

---
Code: HPM stands for Highest posterior Probability Model and
MPM for Median Probability Model.
(a) Based on the R output above, answer the following questions.
i. Write down the full mathematical form of the Bayesian model being fitted, for
the model “reg3”. Also, write down the form of the prior on the regression
coefficients and the values set for the prior hyperparameters.
[2 MARKS]
ii. Comment on the MCMC convergence (as far as you can tell from the output)
for all three models and state whether there could be any cause for concern
regarding any of the parameters.
[1 MARKS]
iii. Using any information you need from the output, calculate the log-Bayes
Factors for each pair of models. Which model, out of “reg1”, “reg2” or
“reg3”, would you choose for this data set, and why?
[4 MARKS]
iv. Explain briefly what the command “Bvs” does. Why might this be a more
convenient way to choose the best model, for this example?
[2 MARKS]
v. Explain clearly what the term “Inclusion Probabilities” refers to, and what
the number “0.0685” (for distance) stands for in this context. [2 MARKS]
vi. Giving numerical evidence, state which model you would choose for this data
set using “Bvs”. Does this model match the one you chose in part (iii)?
[2 MARKS]
(b) Suppose you want to fit a Bayesian multiple regression model with data (Xi, Yi),
(i = 1, . . . , n), where each Xi is a vector of p measurements. You assume that
Y = Xβ + ,
 ∼ N (0, σ2I) (where Y = (Y1, . . . , Yn)0 and X is a matrix with
rows corresponding to X1, . . . , Xn). Next, you assume, a priori, β ∼ N (0, σ2V 0)
(V 0 being a matrix of order p × p), while σ2 ∼ Inv-Gamma(c, d), where if Z ∼
Inv-Gamma(a, b) it means that p(z) ∝ z−a−1e−b/z.
i. Write down an expression proportional to the joint posterior distribution of
the parameters given the data, p(β, σ2|X, Y ).
[3 MARKS]
ii. From the joint posterior distribution in (i), deduce the form of the conditional
posterior distribution of σ2|X, Y , β, identifying it by name, and stating its
parameters.
[4 MARKS]
END OF QUESTION PAPER.
5

Tuesday, 7th May 2019
2.00 p.m. - 3.30 p.m.
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Advanced Bayesian Methods
This paper consists of 6 pages and contains 3 questions.
Candidates should attempt all questions.
Question 1
22 marks
Question 2
18 marks
Question 3
20 marks
Total
60 marks
The following material is made available to you:
Probability formula sheet
Discrete univariate distributions
Distribution
Probability mass function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.m.f.) p(x)
function M(t)
 
Binomial1,2
n
n ∈ N
θx(1

nθ(1
Bi
− θ)n−x
x ∈ {0, 1, . . . , n}
− θ)
(1 − θ + θ exp(t))n
No. of successes in n trials
(n, θ)
x
0 < θ < 1
θ – probability of success
Geometric
1
θ
(1 − θ) exp(t)
No. of trials until (and including)
first failure
Geo
θx−1(1 − θ)
x ∈ N
0 < θ < 1
(θ)
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success


N
No. of type I objects in a sample of
Hypergeometric
M
N −M
− n
x
n−x
N, n ∈ N
nθ(1 − θ) ·
—3
size n, drawn without replacement
HyGe

x ∈ {max{0, n−(N −M)},

N
(n, N, M )
N
− 1
. . . , min{n, M}}
M ∈ {0, . . . , N}
from a population of size N, con-
n
(with θ = M )
N
taining M type I objects.




No. of trials until (and including)
Negative Binomial
x − 1
k
k ∈ N
k

(1 − θ) exp(t)
kth failure
NeBi
θx−k(1 − θ)k
x ∈ {k, k + 1, . . .}
(k, θ)
k − 1
0 < θ < 1
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success
NeBi(1, θ) ≡ Geo(θ)
Poisson
λx
Poi
exp(−λ)
x ∈ N
(λ)
x!
0
λ > 0
λ
λ
exp(λ(exp(t) − 1))
1 Bi(n, θ) can be approximated by Poi(nθ), if n large, θ small and nθ moderate.
2 Bi(n, θ) can be approximated by N(nθ, nθ(1 − θ)), if n large and θ not too close to 0 or 1.
3 No simple closed form expression exists.
Continuous univariate distributions
Distribution
Probability density function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.d.f.) f(x)
function M(t)
Beta
xα1−1(1 − x)α2−1
α
X1 ∼ Ga(α1, θ)
1 > 0
α1
α1α2
—3
X
Be
0

≤ x ≤ 1
2 ∼ Ga(α2, θ) independent
1, α2)
B(α1, α2)
α2 > 0
α1 + α2
(α1 + α2)2(α1 + α2 + 1)

X1
∼ Be(α
X
1, α2)
1+X2
Cauchy
1
—4
—4
—4
Ca
Ca


x ∈ R
η ∈ R
(0, 1) ≡ t(1)
(η, γ)
πγ 1 + (x−η)2
γ > 0
γ2

Chi-Squared
ν
x 2 −1 exp − x
1
2

x > 0
ν ∈ N
ν

Xi ∼ N(0, 1) independent
χ2(ν)
2 ν2 Γ ν
(1 − 2t) ν2
⇒ Pν
X2
2
i=1
i ∼ χ2(ν)
Exponential
1
1
1
Expo
θ exp(−θx)
x > 0
θ > 0
(θ)
θ
θ2
1 − tθ
ν1
ν2
ν
F
ν
2
2 ν2
ν 2
2
1
2 (ν1 + ν2 − 2)
X1 ∼ χ2(ν1)
1
ν2
x 2 −1
—4
X
F

x > 0
ν
ν
2 ∼ χ2(ν2) independent

ν
1, ν2 ∈ N
2 − 2
ν1(ν2 − 2)2(ν2 − 4)
1, ν2)
B ν1 , ν2
1+ν2
2
2
(ν1x + ν2) 2
(for ν
⇒ X1/ν1 ∼ F(ν
2 > 2)
(for ν
1, ν2)
2 > 4)
X2/ν2
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
CONTINUED OVERLEAF/
1

1. (a) Tickets to a famous rock band’s concert in Glasgow next year have been put into
an unspecified number of boxes of a particular cereal. Three friends, A, B and C,
decide to buy boxes of cereal every week until they get enough tickets to attend
the concert. The number of boxes, X, they need to buy to get 3 tickets, is given
by a Negative Binomial distribution NeBi(3, 1 − π), where π, the probability of
getting a concert ticket in a box, is unknown. They manage to get 3 tickets by the
time they buy 20 boxes. Friend A believes that probably less than 1 out of a 1000
boxes contain tickets, or, π < 0.001. However, friend B thinks it likely that many
more tickets for the concert have been put into boxes, and instead says π > 0.02.
Friend C, being a Statistics student, decides to carry out a Bayesian hypothesis
test to find who of A and B is more plausible, a posteriori.
i. Write down the two hypotheses, H0 and H1, for the testing problem.
[2 MARKS]
ii. Assuming a Be(2, 200) prior distribution for π, find the posterior distribution
for π, given X, identifying it by name.
[4 MARKS]
iii. Calculate the prior and posterior odds for H1 vs. H0, assuming that, a priori,
each of the hypotheses are equally likely to be true. Then, write down the
expression for the Bayes factor for testing H1 vs. H0 and explain what decision
you would make regarding H1 and H0. For the purposes of this problem, you
may make use of the following output from R:
> pbeta(0.02, 5, 217)
[1] 0.4534573
> pbeta(0.001, 5, 217)
[1] 3.507448e-06
> pbeta(0.02,2,200)
[1] 0.9120603
> pbeta(0.001,2,200)
[1] 0.0176214
[6 MARKS]
(b) The owner of a computer shop has just received word that an unknown number
of laptops out of 30 delivered from a certain company are defective, leading the
battery to run out after a very short period of time. The owner then decides to
collect data, yi, on how long it takes for laptop i’s battery to run out (i = 1, . . . , 30),
and use a mixture model to determine which laptops were indeed faulty. yi is
assumed to be a realization from a mixture of two exponential distributions with
parameters α and β, with
f (yi|α, β, π) = παe−αyi + (1 − π)βe−βyi,
i = 1, . . . , 30,
and π, (0 < π < 1) is the proportion of defective laptops made by the company.
CONTINUED OVERLEAF/
2

i. For i = 1, . . . , 30, assume that zi is a latent variable which equals 1 if laptop
i is defective, and 0 if not. Assuming Gamma(a, b) priors for both α and β,
and a Beta(c, d) prior for π (a, b, c, d are known numbers), all independent,
write down an expression proportional to the joint posterior distribution of
all parameters (α, β, π) and latent variables z = (z1, . . . , z30). [4 MARKS]
ii. Derive the form of the posterior full conditional distribution of zi|π, α, β, yi
for i = 1, . . . , 30, and describe how it can be used to sample a set of values
of (z(t), . . . , z(t)) in the t-th iteration of a Gibbs sampler, given the values of
1
30
the other parameters from the (t − 1)-th iteration, i.e., (α(t−1), β(t−1), π(t−1)).
[6 MARKS]
2. You suspect that your younger brother is spending too much time on his phone nowa-
days. However, he thinks that he spends a similar amount of time on it as most of
his friends. You therefore decide to carry out a Bayesian model assessment to check
whether his assertion is true. You ask his 10 friends to each record yi, the total time
(in hours) last week they spent on their phones. To this data, you add y11, the hours
spent by your brother on his phone, in the same week (which is the largest number in
the set). Next, you assume that
yi|µ, σ2 ∼ N (µ, σ2),
i = 1, . . . , 11,
and also that p(µ) ∝ 1 and σ2 ∼ Gamma(α, β).
(a) Show that the full conditional posterior distribution p(µ|σ2, y1, . . . , y11) is a normal
distribution, and derive its parameters.
[4 MARKS]
(b) Find an expression proportional to the full conditional posterior distribution of
σ2|µ, y1, . . . , y11 and write down a single step of sampling from it using a Metropolis-
Hastings algorithm, including the expression for the acceptance ratio.
[7 MARKS]
(c) Defining T (y) to be max{yi}, describe the steps you would take to carry out a
posterior predictive check of the validity of the assumption that the data y1, . . . , y11
come from a normal distribution, mentioning in detail any distributions you would
use, and any datasets you simulate in the process.
[7 MARKS]
CONTINUED OVERLEAF/
3

3. (a) Suppose you have data (Xi, Yi), (i = 1, . . . , n), where each Xi is a vector of length
p, to which you want to fit a Bayesian linear regression model,
Y = Xβ + ,
 ∼ N (0, σ2I),
where I denotes the identity matrix, Y = (Y1, . . . Yn)0 and X is a matrix composed
of rows given by X1, . . . , Xn. You choose to use a prior
p(β, σ2) = p(β|σ2)p(σ2),
where β|σ2 ∼ N (0, σ2B0), B0 being a known non-singular matrix of dimension p.
You also choose a prior for σ2, as p(σ2) ∝ (σ2)−a−1e−b/σ2.
i. Write down an expression proportional to the joint posterior distribution of
the parameters given the data, p(β, σ2|X, Y ).
[4 MARKS]
ii. From the joint posterior distribution in (i), you need to find the form of
the conditional posterior distribution p(β|σ2, X, Y ). During your derivation,
suppose you arrive at a point where

1 

p(β|σ2, Y , X) ∝ exp

β0B−1β + Y 0Y − 2Y 0Xβ + β0X0Xβ
.
2σ2
0
By completing the next few steps, clearly showing your reasoning, identify
by name the conditional posterior distribution of β|σ2, Y , X and state its
parameters.
[5 MARKS]
(b) Below, you are provided with a section of the R code and output from an analysis
comparing several Bayesian regression models, fitted with an MCMC sampler
(10,000 iterations in each case) to a dataset relating the fuel consumption of cars
to a number of aspects of their design and performance. The predictors are cyl
(Number of cylinders), hp (horsepower), drat (Rear axle ratio), wt (Weight in 1000
lbs), qsec (1/4 mile time), and gear (Number of forward gears). The response
variable is mpg (Miles per gallon). 9 quantities in the output below are missing,
and are denoted by “(x?)” (“x” being a number between 1 to 9).
> BayesFactor(m1,m2,m3,m4)
...
The matrix of the natural log Bayes Factors is:
m1
m2
m3
m4
m1
0.00
(1?)
(2?) -4.5875
m2
(3?)
0.0000
-6.66
(4?)
m3
(5?)
(6?)
0.00
6.7119
m4
4.59
(7?)
(8?)
0.0000
CONTINUED OVERLEAF/
4

m1 :
call =
MCMCregress(formula = mpg ~ cyl + hp + drat + wt + qsec + gear,
data = cardata, B0 = 0.01, marginal.likelihood = "Laplace")
log marginal likelihood =
-105.3858
m2 :
call =
MCMCregress(formula = mpg ~ wt + hp + qsec, data = cardata, B0 = 0.01,
marginal.likelihood = "Laplace")
log marginal likelihood =
-100.7424
m3 :
call =
MCMCregress(formula = mpg ~ wt + qsec, data = cardata, B0 = 0.01,
marginal.likelihood = "Laplace")
log marginal likelihood =
-94.08635
m4 :
call =
MCMCregress(formula = mpg ~ wt, data = cardata, B0 = 0.01,
marginal.likelihood = "Laplace")
log marginal likelihood =
(9?)
Based on the R output above, answer the following questions.
i. Write down the full mathematical form of the Bayesian model being fitted,
for the model “m3”. Also, write down the form of the prior on the regression
coefficients and the values set for the hyperparameters of this prior.
[2 MARKS]
ii. Using your knowledge of Bayes Factors, and by making any necessary calcu-
lations, write down the values of all the missing quantities denoted as “(x?)”
in the output above.
[5 MARKS]
iii. Which model, out of the four, would you choose for the data, and why?
[2 MARKS]
CONTINUED OVERLEAF/
5

iv. Suppose you additionally get to see the following output, which you had missed
earlier. Would this change your conclusions in the previous part? Give reasons
for your answer.
[2 MARKS]
> effectiveSize(m1)
(Intercept)
cyl
hp
drat
wt
qsec
9623.421
10000.000
10556.160
9555.263
10000.000
10000.000
gear
sigma2
10000.000
5796.755
> effectiveSize(m2)
(Intercept)
wt
hp
qsec
sigma2
8781.439
9922.985
9608.375
8442.052
7221.439
> effectiveSize(m3)
(Intercept)
wt
qsec
sigma2
8484.696
10000.000
9011.905
7653.456
> effectiveSize(m4)
(Intercept)
wt
sigma2
9636.247
9714.818
8017.129
END OF QUESTION PAPER.
6





Friday, 18th May 2018
09.30 am - 11.00 am
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics – Flexible Regression
This paper consists of 6 pages and contains 4 questions.
Candidates should attempt only 3 out of the 4 questions. If all 4 questions
are attempted, candidates should clearly indicate which questions they wish
to be marked. Otherwise, only the first 3 questions in the script book will be
marked.
Every question is worth 20 marks.
The following material is made available to you:
Statistical tables∗
Probability formula sheet
Discrete univariate distributions
Distribution
Probability mass function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.m.f.) p(x)
function M(t)
 
Binomial1,2
n
n ∈ N
θx(1

nθ(1
Bi
− θ)n−x
x ∈ {0, 1, . . . , n}
− θ)
(1 − θ + θ exp(t))n
No. of successes in n trials
(n, θ)
x
0 < θ < 1
θ – probability of success
Geometric
1
θ
(1 − θ) exp(t)
No. of trials until (and including)
first failure
Geo
θx−1(1 − θ)
x ∈ N
0 < θ < 1
(θ)
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success


N
No. of type I objects in a sample of
Hypergeometric
M
N −M
− n
x
n−x
N, n ∈ N
nθ(1 − θ) ·
—3
size n, drawn without replacement
HyGe

x ∈ {max{0, n−(N −M)},

N
(n, N, M )
N
− 1
. . . , min{n, M}}
M ∈ {0, . . . , N}
from a population of size N, con-
n
(with θ = M )
N
taining M type I objects.




No. of trials until (and including)
Negative Binomial
x − 1
k
k ∈ N
k

(1 − θ) exp(t)
kth failure
NeBi
θx−k(1 − θ)k
x ∈ {k, k + 1, . . .}
(k, θ)
k − 1
0 < θ < 1
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success
NeBi(1, θ) ≡ Geo(θ)
Poisson
λx
Poi
exp(−λ)
x ∈ N
(λ)
x!
0
λ > 0
λ
λ
exp(λ(exp(t) − 1))
Statistical Tables
1 Bi(n, θ) can be approximated by Poi(nθ), if n large, θ small and nθ moderate.
2 Bi(n, θ) can be approximated by N(nθ, nθ(1 − θ)), if n large and θ not too close to 0 or 1.
3 No simple closed form expression exists.
Continuous univariate distributions
Distribution
Probability density function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.d.f.) f(x)
function M(t)
Beta
xα1−1(1 − x)α2−1
α
X1 ∼ Ga(α1, θ)
1 > 0
α1
α1α2
—3
X
Be
0

≤ x ≤ 1
2 ∼ Ga(α2, θ) independent
1, α2)
B(α1, α2)
α2 > 0
α1 + α2
(α1 + α2)2(α1 + α2 + 1)

X1
∼ Be(α
X
1, α2)
1+X2
Cauchy
1
—4
—4
—4
Ca
Ca


x ∈ R
η ∈ R
(0, 1) ≡ t(1)
(η, γ)
πγ 1 + (x−η)2
γ > 0
γ2

Chi-Squared
ν
x 2 −1 exp − x
1
2

x > 0
ν ∈ N
ν

Xi ∼ N(0, 1) independent
χ2(ν)
2 ν2 Γ ν
(1 − 2t) ν2
⇒ Pν
X2
2
i=1
i ∼ χ2(ν)
Exponential
1
1
1
Expo
θ exp(−θx)
x > 0
θ > 0
(θ)
θ
θ2
1 − tθ
ν1
ν2
ν
F
ν
2
2 ν2
ν 2
2
1
2 (ν1 + ν2 − 2)
X1 ∼ χ2(ν1)
1
ν2
x 2 −1
—4
X
F

x > 0
ν
ν
2 ∼ χ2(ν2) independent

ν
1, ν2 ∈ N
2 − 2
ν1(ν2 − 2)2(ν2 − 4)
1, ν2)
B ν1 , ν2
1+ν2
2
2
(ν1x + ν2) 2
(for ν
⇒ X1/ν1 ∼ F(ν
2 > 2)
(for ν
1, ν2)
2 > 4)
X2/ν2
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
∗Note for invigilators: will be delivered to the exam venue by the School
CONTINUED OVERLEAF/
1


1. In order to investigate changes in the air quality, the mean sulphur dioxide level for
each month was recorded over several years at a monitoring station in Europe. It
was of interest to investigate the relationship between sulphur dioxide (SO2), the air
temperature, oC, (Temp) and Humidity, %. The following model was fitted
log(SO2)i = β0 + f1(Temp ) + f
) + 
i
2(Humidityi
i
for i = 1, . . . , 471, smooth functions fj(), j = 1, 2, and i ∼ N (0, σ2).
(a) Describe how the backfitting algorithm can be used to fit the model above. Include
in your answer the constraint that is required for the model terms to be identifiable.
[8 MARKS]
(b) After fitting the model within R, using the gam function in the mgcv package, the
following summary output and plots were produced:
Formula:
log(SO2) ~ s(Temp) + s(Humidity)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.11606
0.03889
-2.984
0.003
Approximate significance of smooth terms:
edf Ref.df
F p-value
s(Temp)
3.339
4.214 24.810
<2e-16
s(Humidity) 3.427
4.312
2.663
0.0259
R-sq.(adj) =
0.218
Deviance explained = 22.9%
GCV = 0.72447
Scale est. = 0.71252
n = 471
CONTINUED OVERLEAF/
2

i. Use the summary output and associated plots to comment on the fitted rela-
tionships between log SO2, temperature and humidity.
[6 MARKS]
ii. In the model output above the smoothing parameters have been selected au-
tomatically. One approach for this is by using cross-validation. State the
cross-validation criterion and use the fact that:
y
y
i − ˆ
yi
i − ˆ
y−i ≈ 1 − tr(P/n)
to derive the GCV criterion for an additive model, where ˆ
y−i is the estimate of
yi from fitting the additive model without data point i, and P is a projection
matrix such that ˆ
y = Py.
[4 MARKS]
iii. Explain the bias-variance tradeoff when smoothing is used.
[2 MARKS]
2. The following model for a response y and explanatory variable x uses a truncated
power basis for a polynomial spline of degree p:
K
X
Yi = β0 + β1xi + · · · + βpxp +
β
+ 
i
pk (xi − κk )p
+
i,
i = 1, . . . , n
(1)
k=1
where n is the number of data points, p ≥ 1 is the degree of the polynomial, κ1, . . . , κK
denote the positions of the knots, and (.)p = max {(.)p, 0}.
+
(a) Use equation (1) to write down a model with a cubic spline basis that has 7 degrees
of freedom for the model overall.
[3 MARKS]
(b) Use vector-matrix notation E(y) = Xβ, where X is a matrix of basis functions,
β is a vector of basis coefficients and y is a response vector, to write down the
model fitting criterion used to estimate E(y) using a regression spline, and hence
derive that ˆ
β = (XT X)−1XT y.
[5 MARKS]
(c) Derive an expression for tr(S) where tr is trace and S is a smoothing matrix such
that ˆ
y = Sy, and hence show that tr(S) = dfmod, where dfmod is the degrees of
freedom for the model.
[5 MARKS]
(d) A penalised regresssion spline would enable the level of smoothing to be altered
by adding a penalty term with smoothing parameter λ to the regression spline.
Therefore, the fitting criterion is to minimise:
||y − Xβ||2 + λβT Dβ
CONTINUED OVERLEAF/
3

with respect to β, where D is the penalty matrix, such that
ˆ
β = (XT X + λD)−1XT y.
i. Derive the condition required on the smoothing parameter λ for ˆ
β to be an
unbiased estimator of β, the vector of basis coefficients.
[3 MARKS]
ii. In the method of p-splines the penalty matrix is based on differences of the
coefficients. State the form of the penalty used for first order differences and
derive the penalty matrix, D, for the case of 3 basis coefficients, βj, j = 0, 1, 2.
[4 MARKS]
3. Consider the following nonparametric regression model:
M 1 : Yi = f (xi) + i,
i = 1, . . . , n, i ∼ N (0, σ2), and f () is a smooth function.
(a) Describe the model fitting criterion for estimating the function f (x) in M 1 at a
target point x using local linear regression.
[4 MARKS]
(b) Assuming ˆ
y = Sy, where S is a smoothing matrix, explain how a variability band
is produced for ˆ
f (x) at a target point x.
[4 MARKS]
(c) An alternative way to fit the smooth function could be to use a smoothing splines
approach. Write down the fitting criterion for smoothing splines and hence explain
the difference between how the level of smoothing is determined in local linear
regression and for smoothing splines.
[4 MARKS]
(d) State two other parametric models that it could be of interest to compare M 1
with to test whether or not the smooth function provides useful information (i.e.
is there an ‘effect’) and, whether or not the smooth function should be a linear
relationship, respectively.
[2 MARKS]
(e) Explain how you would compare M 1 to each of the two models that you have
stated in part (d) using an approximate F-test. In particular:
i. state the hypotheses of interest in terms of the models being compared;
ii. state the F-statistic;
iii. explain how the F-statistic is compared to an appropriate F-distribution and
how the results are interpreted.
[6 MARKS]
CONTINUED OVERLEAF/
4

4. (a) Environmental regulators were interested in assessing if the water quality in river
water was above or below the environmental quality standard. In order to inves-
tigate this the water quality was recorded at n different sites and each site was
classified as ‘pass’ or ‘fail’. It was of interest to model the response variable of
pass/fail as binomial and to explore relationships with a smooth function of total
monthly rainfall and a smooth function of distance to the nearest sewage works.
State the equation for a generalised additive model for this context, with a response
of ‘pass’ or ‘fail’ (0, 1), a smooth covariate for rainfall and a smooth covariate for
distance.
[2 MARKS]
(b) The full model described in part (a) was fitted in R. The results from fitting the
full model in part (a), Model A, and two reduced models are contained in the
table below. Specifically, the table displays the AICs, degrees of freedom for error
(dferr), deviance for each model, and a p-value for one model comparison.
Model (covariates)
AIC
dferr
Deviance
p-value
A (s(rain), s(distance)):
1186.7
462
330.07
B (s(rain)):
1277.97
466
407.87
0.011 (A vs B)
C (s(distance)):
1193.9
466
340.91
? (A vs C)
i. The model comparison for Models A and C is missing in the table above.
Use a chi-squared test to compare the deviances and hence formally compare
Models A and C. (Note: you do not need to compute the p-value).
[4 MARKS]
ii. Use your answer to part (i) along with the rest of the table above to comment
on the most appropriate model. Comment on both the model comparisons
and the AIC values.
[4 MARKS]
(c) The investigation team then decided to fit an alternative model with a response
of the phosphorus concentration, the same two smooth covariates and assuming
gaussian errors. State an equation for this new model.
[2 MARKS]
(d) State the fitting criterion for the model in part (c) in vector matrix form, E(y) =
Xβ, using penalised regression splines and hence derive that the estimates of the
basis coefficients can be estimated as:
ˆ
β = (XT X + R)−1XT y
where X is a matrix of basis functions and β is a vector of basis coefficients, for
a response y.
Additionally, state the form of R in the above expression.
[5 MARKS]
CONTINUED OVERLEAF/
5

(e) Define the partial residuals for covariate 1, rainfall, in the context of an additive
model, and explain how they are used and why they are useful in this context.
[3 MARKS]
END OF QUESTION PAPER.
6

Tuesday, 23 April 2019
2.00p.m. - 3.30 p.m.
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics – Flexible Regression
This paper consists of 6 pages and contains 4 questions.
Candidates should attempt only 3 out of the 4 questions. If all 4 questions
are attempted, candidates should clearly indicate which questions they wish
to be marked. Otherwise, only the first 3 questions in the script book will be
marked.
Question 1
20 marks
Question 2
20 marks
Question 3
20 marks
Question 4
20 marks
Total
80 marks
The following material is made available to you:
Statistical tables∗
Probability formula sheet
Discrete univariate distributions
Distribution
Probability mass function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.m.f.) p(x)
function M(t)
 
Binomial1,2
n
n ∈ N
θx(1

nθ(1
Bi
− θ)n−x
x ∈ {0, 1, . . . , n}
− θ)
(1 − θ + θ exp(t))n
No. of successes in n trials
(n, θ)
x
0 < θ < 1
θ – probability of success
Geometric
1
θ
(1 − θ) exp(t)
No. of trials until (and including)
first failure
Geo
θx−1(1 − θ)
x ∈ N
0 < θ < 1
(θ)
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success


N
No. of type I objects in a sample of
Hypergeometric
M
N −M
− n
x
n−x
N, n ∈ N
nθ(1 − θ) ·
—3
size n, drawn without replacement
HyGe

x ∈ {max{0, n−(N −M)},

N
(n, N, M )
N
− 1
. . . , min{n, M}}
M ∈ {0, . . . , N}
from a population of size N, con-
n
(with θ = M )
N
taining M type I objects.




No. of trials until (and including)
Negative Binomial
x − 1
k
k ∈ N
k

(1 − θ) exp(t)
kth failure
NeBi
θx−k(1 − θ)k
x ∈ {k, k + 1, . . .}
(k, θ)
k − 1
0 < θ < 1
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success
NeBi(1, θ) ≡ Geo(θ)
Poisson
λx
Poi
exp(−λ)
x ∈ N
(λ)
x!
0
λ > 0
λ
λ
exp(λ(exp(t) − 1))
Statistical Tables
1 Bi(n, θ) can be approximated by Poi(nθ), if n large, θ small and nθ moderate.
2 Bi(n, θ) can be approximated by N(nθ, nθ(1 − θ)), if n large and θ not too close to 0 or 1.
3 No simple closed form expression exists.
Continuous univariate distributions
Distribution
Probability density function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.d.f.) f(x)
function M(t)
Beta
xα1−1(1 − x)α2−1
α
X1 ∼ Ga(α1, θ)
1 > 0
α1
α1α2
—3
X
Be
0

≤ x ≤ 1
2 ∼ Ga(α2, θ) independent
1, α2)
B(α1, α2)
α2 > 0
α1 + α2
(α1 + α2)2(α1 + α2 + 1)

X1
∼ Be(α
X
1, α2)
1+X2
Cauchy
1
—4
—4
—4
Ca
Ca


x ∈ R
η ∈ R
(0, 1) ≡ t(1)
(η, γ)
πγ 1 + (x−η)2
γ > 0
γ2

Chi-Squared
ν
x 2 −1 exp − x
1
2

x > 0
ν ∈ N
ν

Xi ∼ N(0, 1) independent
χ2(ν)
2 ν2 Γ ν
(1 − 2t) ν2
⇒ Pν
X2
2
i=1
i ∼ χ2(ν)
Exponential
1
1
1
Expo
θ exp(−θx)
x > 0
θ > 0
(θ)
θ
θ2
1 − tθ
ν1
ν2
ν
F
ν
2
2 ν2
ν 2
2
1
2 (ν1 + ν2 − 2)
X1 ∼ χ2(ν1)
1
ν2
x 2 −1
—4
X
F

x > 0
ν
ν
2 ∼ χ2(ν2) independent

ν
1, ν2 ∈ N
2 − 2
ν1(ν2 − 2)2(ν2 − 4)
1, ν2)
B ν1 , ν2
1+ν2
2
2
(ν1x + ν2) 2
(for ν
⇒ X1/ν1 ∼ F(ν
2 > 2)
(for ν
1, ν2)
2 > 4)
X2/ν2
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
∗Note for invigilators: will be delivered to the exam venue by the School
CONTINUED OVERLEAF/
1

1. Consider the model, for data (xi, yi), i = 1, . . . , n:
Yi = f (xi, β) + i,
i ∼ N(0, σ2)
where f is a non-linear vector valued function of the parameters β.
(a) Derive the iterative steps that would be used to fit this non-linear model using
non-linear least squares.
[7 MARKS]
(b) Alternatively a nonparametric regression model could be fitted to the data (xi, yi),
i = 1, . . . , n i.e.
Yi = f (xi) + i,
i ∼ N(0, σ2).
Assuming local quadratic regression is the smoother being used,
i. State the model fitting criterion that is used to estimate a function f (x) at
the target point x, and explain how an estimate of the smooth function is
obtained at x.
[3 MARKS]
ii. Write the local quadratic fitting criterion, from part (i), in vector-matrix no-
tation for target point x, where y is the response data, X is the design matrix,
β is the vector of parameters and W is a weight matrix, and hence state an
expression for the estimator of f (x) at the point x in vector-matrix notation.
[4 MARKS]
iii. For a covariate vector, x, give the formulation of the smoothing matrix, S,
that results in ˆ
y = f (x) = Sy.
[2 MARKS]
iv. Assuming ˆ
y = Sy, where S is the smoothing matrix, derive variability bands
for f (x) at the point x.
[4 MARKS]
2. Consider the following nonparametric regression model:
Yi = f (xi) + i,
i ∼ N(0, σ2)
for data (xi, yi), i = 1, . . . , n where f is a smooth function of the covariate x.
(a) One way of estimating f () would be through some spline-based function. Two
commonly used methods are smoothing splines and penalised regression splines,
which both have a penalty term, often expressed as:
Z
x(n)
λ
[f 00(x)]2 dx.
x(1)
CONTINUED OVERLEAF/
2

Explain the role of the penalty term detailed above and give one reason why
penalised regression splines could be preferred to smoothing splines.
[3 MARKS]
(b) For penalised regression splines, f (x) is expressed in terms of basis functions and
basis coefficients, as:
X
f (x) =
βjbj(x)
j
where βj is the basis function corresponding to the jth basis function, bj().
Use this to show how the penalty term, given in part (a), can be expressed as
λβ>Dβ where β is a vector of basis coefficients, λ is the smoothing parameter
and D a penalty matrix.
[4 MARKS]
(c) Using part (b), write down, in vector-matrix form, the model fitting criterion for
a penalised regression spline, and hence minimise this criterion with respect to the
basis coefficients, β, to obtain an expression for the estimated basis coefficients ˆ
β.
[6 MARKS]
(d) A common choice for the basis functions are B-splines due to their computational
efficiency, explain the property of B-splines which makes them computationally
efficient.
[2 MARKS]
(e) The following three models were fitted to a dataset containing 100 observations:
M1 :
Yi = β0 + i
M2 :
Yi = β0 + β1xi + i
M3 :
Yi = f (xi) + i
where i = 1, . . . , n and i ∼ N (0, σ2). Model 3 was fitted using cubic B-splines
and 7 degrees of freedom.
i. Using these models, state the set of hypotheses that could be compared to
determine whether a nonparametric relationship provided more information
than a linear relationship.
[2 MARKS]
CONTINUED OVERLEAF/
3

ii. Running the anova function in R to compare a linear relationship with a
nonparametric relationship gives the following incomplete R output:
Analysis of Variance Table
Model A: y ~ x
Model B: y ~ bs(x, df = 6)
Res.Df
RSS
A
C
74.498
B
93 70.903
State the value of C and hence use an approximate F-test, the incomplete R
output and your hypotheses from part (i) to investigate if a nonparametric
relationship provides more information than a linear relationship.
[3 MARKS]
3. Estate agents in Chicago were interested in assessing whether socio-economic charac-
teristics were related to median house prices across the city. For each of Chicago’s
76 communities the following variables were recorded: Median house price (MEDHP),
Number of homicides per 1000 residents (HOMICIDE), Percentage of families below
the poverty line (POVERTY) and Percentage of houses whose highest earner was in a
managerial role (MANAGERIAL).
The following additive model was fitted:
log (
e MEDHPi) = β0 + f1(HOMICIDEi) + f2(POVERTYi) + f3(MANAGERIALi) + i
where i = 1, . . . , 76, i ∼ N (0, σ2) are independent and fj() is a smooth function
j = 1, . . . , 3.
(a) Assuming the model is fitted using the penalised regression spline approach, state
the fitting criterion in this context.
[2 MARKS]
(b) Using your answer to part (a), derive the projection matrix (P) for the penalised
regression spline approach in this context.
[5 MARKS]
(c) State a general expression for the degrees of freedom of each of the covariates and
hence state an expression for the overall degrees of freedom of the model.
[2 MARKS]
(d) After fitting the model in R using the gam function in the mgcv library, the output
overleaf was obtained:
CONTINUED OVERLEAF/
4


Family: gaussian
Link function: identity
Formula:
log(MEDHP) ~ s(MANAGERIAL) + s(POVERTY) + s(HOMICIDE)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.19067
0.02443
499
<2e-16
---
Approximate significance of smooth terms:
edf Ref.df
F
p-value
s(MANAGERIAL) 1.752
2.177 45.195 1.18e-15
s(POVERTY)
1.000
1.000
0.143
0.706
s(HOMICIDE)
2.367
2.885
9.568 2.96e-05
---
R-sq.(adj) =
0.771
Deviance explained = 78.7%
GCV = 0.049328
Scale est. = 0.045356
n = 76
i. Use the summary output and plots above to comment on the fit of the model
and the relationships between each of the covariates and the response.
[6 MARKS]
CONTINUED OVERLEAF/
5

ii. The smoothing parameters were selected here via generalised cross validation
(GCV), alternatively cross validation (CV) could have been used. Explain an
advantage of using GCV instead of CV for smoothing parameter selection and
give an alternative method for selecting the smoothing parameters.
[3 MARKS]
iii. Here penalised regression splines have been used and the level of smoothness
has been controlled by the smoothing parameters, give two further quantities
that can be altered to control the smoothness of the function when penalised
regression splines are being used.
[2 MARKS]
4. Colic is a gastrointestinal condition in horses which is generally caused by an ob-
struction of the intestine and is potentially fatal. Colic can sometimes be fixed by
surgical intervention, however veterinarians believe there are other characteristics of
the horses health which are related to whether or not the horse eventually dies due to
the condition.
To assess their beliefs a generalised additive model was fitted with a binomial response
(Outcome) of ‘died’ (including euthanasia) (1) or ‘lived’ (0), smooth functions for
temperature in degrees Celsius (Temperature) and pulse in beats per minute (Pulse)
and a factor effect for if a surgical intervention was performed (Surgery) with levels
‘yes’ (0) and ‘no’(1).
(a) State the equation for a generalised additive model in this context. [3 MARKS]
(b) A reduced model could also be fitted which only contains smooth terms for Pulse
and Temperature, assume this model is called ‘Model 2’. Assuming the model
stated in the problem description is referred to as ‘Model 1’, use the table below
of deviances and effective degrees of freedom for error to compare Models 1 and 2
via a Chi-squared test.
[4 MARKS]
Model
dferr
Deviance
Model 1
220
222.3031
Model 2
221
222.6807
(c) Why is it appropriate to compare the change in deviances to a χ2 distribution
here?
[2 MARKS]
(d) One way of fitting these models would be by the local scoring procedure. Describe
the steps involved in using the local scoring procedure for Model 2 i.e. a generalised
additive model with 2 smooth covariates.
[8 MARKS]
(e) Explain what is meant by the term concurvity and describe an issue that it may
pose when fitting a generalised additive model.
[3 MARKS]
END OF QUESTION PAPER.
6





May 2018
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics – Linear Mixed Models
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
Candidates should attempt any three questions.
NOTE: If all four questions are attempted, candidates should clearly indi-
cate which questions they wish to be marked. Otherwise, only the first three
questions in the script book will be marked
1. A commercial baker wants to develop the perfect lemon sponge cake. The perfect cake
should be fluffy and moist, and should produce a nice clean slice when you cut into it.
A cake which is too dry will crumble when cut, while a cake which is too moist will
stick to the knife, both of which result in the loss of cake. The quality of the cake can
therefore be measured by the mass of cake loss (in grams) when cut, with a loss of 0
grams representing a perfect cake. She wishes to choose between two recipes (A and
B) to identify which produces the lowest amount of cake loss. She randomly selects
4 of her chefs to take part in the experiment, and each chef is asked to cook 5 cakes
using each recipe. Each of the 40 cakes is then cut and the cake loss is measured.
The dataset contains the variables loss (g of cake lost), recipe (Recipe A or B) and
chef (chef who cooked the cake, taking values 1, 2, 3, 4).
CONTINUED OVERLEAF/
1

(a) Write down a suitable model for these data, clearly stating all assumptions.
[6 MARKS]
(b) This model was fitted, and we obtained the table below.
Source
DF
SS
recipe
I-1
5831
chef
J-1
434.58
recipe*chef
(I-1)(J-1)
401.29
error
IJ(K-1)
1032
Test the hypothesis that there is no difference between her two recipes in terms
of cake loss. State your conclusion in terms of the model parameters, and then
relate this back to the practical example.
One or more of the following distributional results may be useful:
F (3, 16; 0.95) = 3.24
F (3, 1; 0.95) = 215.71
F (2, 4; 0.95) = 6.94
F (1, 3; 0.95) = 10.13
F (24, 1; 0.95) = 249.05
F (4, 2; 0.95) = 19.25
[5 MARKS]
(c) Provide method of moments point estimates for each of the three variance com-
ponents in the model. The following equations may be useful:
E(M SB) = Var(error) + 5 Var(recipe*chef) + 10 Var(chef)
E(M SAB) = Var(error) + 5 Var(recipe*chef)
[4 MARKS]
(d) Let δ represent the mean difference in cake loss between Recipe A and Recipe B.
Derive an expression for the variance of δ.
[3 MARKS]
(e) With reference to this example, explain the difference between crossed and nested
factors.
[2 MARKS]
2. Scientists carried out an experiment to test whether a new protein supplement could
improve athletic performance. A group of 50 amateur runners were recruited; 25 were
given a regular four week course of the new supplement, while the remaining athletes
were given a four week course of a protein-free supplement (placebo). Prior to starting
the experiment, each athlete was asked to complete a 5km run, and their baseline time
was recorded. Each week during the experiment, they did another timed 5km run to
monitor their progress.
CONTINUED OVERLEAF/
2

(a) A model was fitted to this data using the following SAS code.
proc mixed;
class id week supp;
model runtime = basetime week supp week*supp;
repeated week / type=un subject=id;
run;
Write down the mean model corresponding to the above code.
[5 MARKS]
(b) The errors in this model are assumed to follow a multivariate normal distribution
with mean 0 and variance-covariance matrix R. Give the form that R takes when
each of the following covariance structures are used.
i. unstructured (type=un)
ii. AR(1) (type=ar)
iii. compound symmetry (type=cs)
[8 MARKS]
(c) The mean model in part (a) was fitted with the three different covariance structures
mentioned. Based on the fit statistics below, which covariance structure(s) are the
best candidates? If you have more than one candidate structure, carry out a formal
test to compare them.
un
ar(1)
cs
-2 Res Log Likelihood
762.1
784.2
791.7
AIC (Smaller is Better)
782.1
786.2
793.7
AICC (Smaller is Better)
783.4
786.3.
793.8
BIC (Smaller is Better)
793.1
790.2.
791.8
One or more of the following distributional results may be useful.
χ2(8; 0.95) = 15.5
χ2(9; 0.95) = 16.9
χ2(10; 0.95) = 18.3.
[4 MARKS]
(d) The general linear mixed model takes the form
y = Xβ + Zu + e
where X and Z are given matrices and
 u 

θ 
 u 
 G 0 
E
=
, Var
=
.
e
λ
e
0
R
CONTINUED OVERLEAF/
3

This can be rewritten as a multivariate normal distribution of the form
y ∼ N (Xβ, V).
i. Complete the above expression by providing the correct values for θ and λ.
ii. Write an expression for the overall variance V in terms of the between-subject
variance G and the error variance R.
[3 MARKS]
3. Suppose we have data (xij, Yij) (where i labels subjects) for which we consider the
model for Yij given xij of the form
M0 : Yij = β0 + β1xij + b0i + b1ixij + eij,
where β0 , β1 are unknown parameters; b0i, b1i are random effects; and eij is the error
term.
(a)
i. Write down the distributional assumptions for eij, b0i and b1i in the context of
a normal linear mixed model with (possibly correlated) random coefficients.
ii. Determine E(Y

ij |xij ) and Var Yij |xij
under these assumptions.
[4 MARKS]
(b) Identify the following three special cases in terms of the covariance parameters in
the model.
i. M1: ‘independent random intercept and slope for each subject’;
ii. M2: ‘random intercept for each subject’; and
iii. M3: ‘no random subject effects’.
[3 MARKS]
(c) Provide a rough sketch of some subject-specific regression lines to illustrate the
differences between models M1, M2 and M3.
[3 MARKS]
(d) The diameters of 12 apples on a particular tree were measured every week over
a six week period as part of an agricultural expirement. Let yij be the diameter
measurement from the ith apple in the jth week. Random coefficient models M0,
M1 and M2 were all fitted in R, and the following model comparisons were carried
out:
CONTINUED OVERLEAF/
4

> anova(m0,m1)
Data: apple
Models:
m1: Diam ~ Time + (1 | AppleID) + (0 + Time | AppleID)
m0: Diam ~ Time + (Time | AppleID)
Df
AIC
BIC
logLik deviance
Chisq Chi Df Pr(>Chisq)
m1
5 121.37 132.75 -55.684
111.37
m0
6 121.78 135.44 -54.887
109.78 1.5934
1
0.2068
> anova(m1,m2)
Data: apple
Models:
m2: Diam ~ Time + (1 | AppleID)
m1: Diam ~ Time + (1 | AppleID) + (0 + Time | AppleID)
Df
AIC
BIC
logLik deviance
Chisq Chi Df Pr(>Chisq)
m2
4 137.91 147.02 -64.957
129.91
m1
5 121.37 132.75 -55.684
111.37 18.546
1
1.659e-05 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ?’ 1
i. Based on these results, which of the three models would you select, and why?
ii. Which assumption is violated when comparing models M1 and M2? Why is
this not a problem in this particular case?
[6 MARKS]
(e) Output from fitting model M1 is shown below.
> summary(m1)
Linear mixed model fit by REML [’lmerMod’]
Formula: Diam ~ Time + (1 | AppleID) + (0 + Time | AppleID)
Data: apple
REML criterion at convergence: 117.8
Scaled residuals:
Min
1Q
Median
3Q
Max
-4.2206 -0.0257
0.0354
0.0907
2.5608
Random effects:
Groups
Name
Variance Std.Dev.
AppleID
(Intercept) 0.03776
0.1943
AppleID.1 Time
0.03028
0.1740
Residual
0.16775
0.4096
Number of obs: 72, groups:
AppleID, 12
CONTINUED OVERLEAF/
5

Fixed effects:
Estimate Std. Error t value
(Intercept)
2.82772
0.12354
22.889
Time
-0.04800
0.05764
-0.833
Correlation of Fixed Effects:
(Intr)
Time -0.393
> ranef(m1)
$AppleID
(Intercept)
Time
1
0.025454395
0.06123113
4
0.018182339
0.07151399
5
-0.007426805
0.04424870
10
0.006071866
0.06182190
11 -0.010692557
0.04920600
13
0.033039132
0.07698327
14 -0.269180731 -0.47939795
15
0.063014383
0.10350599
17
0.149847413 -0.16014983
18 -0.006403647
0.04993039
19
0.003546664
0.06641095
25 -0.005452450
0.05469548
i. Predict the diameter of a new, unobserved apple on this tree in week 4.
ii. Predict the diameter of the apple with AppleID = 10 after 3 weeks.
[4 MARKS]
4. Researchers carried out a study to explore the effects of a new drug on arthritis sever-
ity. A group of 527 patients with arthritis were recruited, and the severity of their joint
pain was observed over a six month period. To simplify the analysis, the severity score
was dichotomised into two categories, “no/minimal joint pain” and “moderate/severe
joint pain”.
This outcome variable, pain, was coded as 0 for “no/minimal pain” and 1 for “mod-
erate/severe pain”. Each patient had their baseline pain level observed and were
randomly assigned either to a placebo or treatment group. The patients were then
CONTINUED OVERLEAF/
6

monitored once a month for 6 months.
The variables id (patient id number), treat (treatment status; 0 for placebo, 1 for
drug) and month (month number; 0 to 6, with 0 indicating baseline values) were
recorded.
The following SAS code was used to fit a GEE model to the data:
proc genmod data=arthritisstudy;
class pain id treat (param=ref ref=’0’);
model pain = treat month treat*month / dist=bin type3;
repeated subject = id/ corrw modelse type= ar(1) printmle;
run;
Selected output from the procedure is shown below:
Analysis Of Initial Parameter Estimates
Parameter
DF
Estimate
Standard
Wald 95%
Wald
Error
Confidence Limits
Chi-Square
Pr> ChiSq
Intercept
1
3.7186
0.4351
2.8658
4.5714
73.04
<.0001
treat 1
1
-0.4091
0.4783
-1.3466
0.5284
0.73
0.3929
month
1
-1.1208
0.2283
-1.5683
-0.6733
24.10
<.0001
month*treat 1
1
-0.4113
0.2581
-0.9172
0.0946
2.54
0.1110
Scale
0
1.0000
0.0000
1.0000
1.0000
Working Correlation Matrix
Col1
Col2
Col3
Col4
Col5
Row1
1.0000
0.3581
0.1282
0.0459
0.0164
Row2
0.3581
1.0000
0.3581
0.1282
0.0459
Row3
0.1282
0.3581
1.0000
0.3581
0.1282
Row4
0.0459
0.1282
0.3581
1.0000
0.3581
Row5
0.0164
0.0459
0.1282
0.3581
1.0000
GEE Fit Criteria
QIC
1372.4128
QICu
1370.3721
CONTINUED OVERLEAF/
7

Analysis Of GEE Parameter Estimates
Empirical Standard Error Estimates
Parameter
Estimate
Standard
95%
Error
Confidence Limits
Z
Pr> |Z|
Intercept
3.7201
0.4805
2.7783
4.6619
7.74
<.0001
treat 1
-0.3912
0.5233
-1.4169
0.6345
-0.75
0.4547
month
-1.1289
0.2497
-1.6183
-0.6395
-4.52
<.0001
month*treat 1
-0.4357
0.2671
-0.9592
0.0878
-1.63
0.1028
Analysis Of GEE Parameter Estimates
Model-Based Standard Error Estimates
Parameter
Estimate
Standard
95%
Error
Confidence Limits
Z
Pr> |Z|
Intercept
3.7201
0.5063
2.7278
4.7124
7.35
<.0001
treat 1
-0.3912
0.5489
-1.4670
0.6846
-0.71
0.4760
month
-1.1289
0.2415
-1.6022
-0.6556
-4.67
<.0001
month*treat 1
-0.4357
0.2673
-0.9596
0.0882
-1.63
0.1031
Scale
1.0000
.
.
.
.
.
Note: The scale parameter was held fixed.
Score Statistics For Joint Tests For GEE
Source
DF
Chi-Square
Pr > ChiSq
treat
1
0.77
0.3795
month
1
19.87
<.0001
month*treat
1
2.17
0.1406
(a) The outcome variable contained 1701 missing values (out of a total of 527×7=3689).
What assumption must we make about this missing data for our GEE analysis to
be valid? Make sure you explain clearly what this means about the relationship
between the missing and observed data.
[3 MARKS]
(b) In the model, what indicates that a person who is designated as having moder-
ate/severe joint pain is likely to continue to remain in that designation?
[2 MARKS]
CONTINUED OVERLEAF/
8

(c) The researchers thought that there might not be a time decaying structure to the
correlation and decided to fit an exchangeable GEE model to the data. The fit
statistics for the exchangeable model are given below. Which model would you
choose and why?
GEE Fit Criteria
QIC
1372.8352
QICu
1370.4420
[2 MARKS]
(d) Does it seem like the probability of having moderate/severe pain over time is
different for different treatment groups?
[2 MARKS]
(e) Partial output from the model without interaction term is given below. Interpret
the month and treat coefficient estimates in terms of odds of having moder-
ate/severe pain.
Analysis Of GEE Parameter Estimates
Empirical Standard Error Estimates
Parameter
Estimate
Standard
95%
Error
Confidence Limits
Z
Pr> |Z|
Intercept
4.4522
0.2797
3.9040
5.0004
15.92
<.0001
treat 1
-1.2481
0.2271
-1.6932
-0.8030
-5.50
<.0001
month
-1.4965
0.0928
-1.6784
-1.3146
-16.13
<.0001
[4 MARKS]
(f) Using the output in part (e), what proportion of patients assigned to the drug are
expected to experience moderate/severe pain in month 3?
[2 MARKS]
(g) Give the name of an alternative model to GEE discussed in this course for analysing
repeated data with non-normal distributions. How does it differ from GEE? In
what situations might we prefer one or the other?
[5 MARKS]
END OF QUESTION PAPER.
9

Tuesday, 7th May 2019
09.30 a.m. - 11.00 a.m.
EXAMINATION FOR THE DEGREES OF M.A., M.SCI. AND B.SC.
(SCIENCE)
Statistics – Linear Mixed Models
This paper consists of 7 pages and contains 4 questions.
Candidates should attempt only 3 out of the 4 questions. If all 4 questions
are attempted, candidates should clearly indicate which questions they wish
to be marked. Otherwise, only the first 3 questions in the script book will be
marked.
Question 1
20 marks
Question 2
20 marks
Question 3
20 marks
Question 4
20 marks
Total
80 marks
The following material is made available to you:
Probability formula sheet
Discrete univariate distributions
Distribution
Probability mass function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.m.f.) p(x)
function M(t)
 
Binomial1,2
n
n ∈ N
θx(1

nθ(1
Bi
− θ)n−x
x ∈ {0, 1, . . . , n}
− θ)
(1 − θ + θ exp(t))n
No. of successes in n trials
(n, θ)
x
0 < θ < 1
θ – probability of success
Geometric
1
θ
(1 − θ) exp(t)
No. of trials until (and including)
first failure
Geo
θx−1(1 − θ)
x ∈ N
0 < θ < 1
(θ)
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success


N
No. of type I objects in a sample of
Hypergeometric
M
N −M
− n
x
n−x
N, n ∈ N
nθ(1 − θ) ·
—3
size n, drawn without replacement
HyGe

x ∈ {max{0, n−(N −M)},

N
(n, N, M )
N
− 1
. . . , min{n, M}}
M ∈ {0, . . . , N}
from a population of size N, con-
n
(with θ = M )
N
taining M type I objects.




No. of trials until (and including)
Negative Binomial
x − 1
k
k ∈ N
k

(1 − θ) exp(t)
kth failure
NeBi
θx−k(1 − θ)k
x ∈ {k, k + 1, . . .}
(k, θ)
k − 1
0 < θ < 1
1 − θ
(1 − θ)2
1 − θ exp(t)
θ – probability of success
NeBi(1, θ) ≡ Geo(θ)
Poisson
λx
Poi
exp(−λ)
x ∈ N
(λ)
x!
0
λ > 0
λ
λ
exp(λ(exp(t) − 1))
1 Bi(n, θ) can be approximated by Poi(nθ), if n large, θ small and nθ moderate.
2 Bi(n, θ) can be approximated by N(nθ, nθ(1 − θ)), if n large and θ not too close to 0 or 1.
3 No simple closed form expression exists.
Continuous univariate distributions
Distribution
Probability density function
Range
Parameters
E(X)
Var(X)
Moment-generating
Comments
(p.d.f.) f(x)
function M(t)
Beta
xα1−1(1 − x)α2−1
α
X1 ∼ Ga(α1, θ)
1 > 0
α1
α1α2
—3
X
Be
0

≤ x ≤ 1
2 ∼ Ga(α2, θ) independent
1, α2)
B(α1, α2)
α2 > 0
α1 + α2
(α1 + α2)2(α1 + α2 + 1)

X1
∼ Be(α
X
1, α2)
1+X2
Cauchy
1
—4
—4
—4
Ca
Ca


x ∈ R
η ∈ R
(0, 1) ≡ t(1)
(η, γ)
πγ 1 + (x−η)2
γ > 0
γ2

Chi-Squared
ν
x 2 −1 exp − x
1
2

x > 0
ν ∈ N
ν

Xi ∼ N(0, 1) independent
χ2(ν)
2 ν2 Γ ν
(1 − 2t) ν2
⇒ Pν
X2
2
i=1
i ∼ χ2(ν)
Exponential
1
1
1
Expo
θ exp(−θx)
x > 0
θ > 0
(θ)
θ
θ2
1 − tθ
ν1
ν2
ν
F
ν
2
2 ν2
ν 2
2
1
2 (ν1 + ν2 − 2)
X1 ∼ χ2(ν1)
1
ν2
x 2 −1
—4
X
F

x > 0
ν
ν
2 ∼ χ2(ν2) independent

ν
1, ν2 ∈ N
2 − 2
ν1(ν2 − 2)2(ν2 − 4)
1, ν2)
B ν1 , ν2
1+ν2
2
2
(ν1x + ν2) 2
(for ν
⇒ X1/ν1 ∼ F(ν
2 > 2)
(for ν
1, ν2)
2 > 4)
X2/ν2
“Hand calculators with simple basic functions (log, exp, square root, etc.) may be used
in examinations. No calculator which can store or display text or graphics may be used,
and any student found using such will be reported to the Clerk of Senate”.
CONTINUED OVERLEAF/
1

1. A statistics lecturer has two pet rabbits who absolutely love eating hay. There are
three main types of hay available for rabbits - alfalfa hay, meadow hay and timothy
hay. He decides to carry out an experiment to investigate which of these types of hay
rabbits prefer. A total of 10 rabbits are recruited to the study, and each rabbit is given
each type of hay for a total of three days. The order in which the rabbits are fed the
types of hay is randomly decided in advance. Each day, we measure the total weight
of hay (in grams) eaten by each rabbit, giving a total of 90 observations (10 rabbits,
9 days).
(a) Write down a suitable model for this experiment, clearly stating all assumptions.
[6 MARKS]
(b) With reference to this example, explain the difference between a nested and a
crossed experiment.
[2 MARKS]
(c) Using the sums of squares below, test the hypothesis that there is no difference
between the types of hay in terms of preference. State your conclusion both in
terms of the model parameters and in relation to the original example.
Source
SS
hay
93192
rabbit
26136
hay:rabbit
67301
error
8641
One or more of the following distributional results may be useful:
F (2, 9; 0.95) = 4.26
F (2, 18; 0.95) = 3.55
F (2, 60; 0.95) = 3.15
F (3, 9; 0.95) = 3.86
F (3, 10; 0.95) = 3.71
F (3, 90; 0.95) = 2.70
[6 MARKS]
(d) Provide unbiased (method of moments) point estimates for each of the variance
components in your model. [Note that it is possible for these estimates to be
negative]
[6 MARKS]
2. Suppose we have data (xij, Yij) where i indicates a grouping/clustering level and j
indicates measurements within groupings. We consider a model of the form:
Yij = β0 + β1xij + b∗ + b∗ x
0i
1i
ij + eij
where β0, β1 are unknown parameters; b∗ , b∗ are random effects; and e
0i
1i
ij is the error
term.
(a) This is known as a random coefficient model. Explain how it differs from an
analysis of covariance (ANCOVA).
[2 MARKS]
CONTINUED OVERLEAF/
2

(b) Write down all of the distributional assumptions for eij, b∗ and b∗ in the context
0i
1i
of this random coefficient model.
[4 MARKS]
(c) Determine E(Y

ij |xij ) and Var Yij |xij .
[2 MARKS]
(d) Creatine is a supplement widely used by bodybuilders and athletes to build muscle
strength. In a study, 42 subjects were were asked to take 7 grams of creatine per
day over a 5-week period. The subjects had their muscle strength measured at the
start of the study and then again every week until the end of the study, giving a
total of 6 measurements per subject.
The following four models were fitted to this study data using the lme4 library in
R.
m1 <- lmer(Muscle ~ Time + (Time|Subject))
m2 <- lmer(Muscle ~ Time + (1|Subject)+(0+Time|Subject))
m3 <- lmer(Muscle ~ Time + (1|Subject))
m4 <- lm(Muscle ~ Time)
Define each of these four models, both in terms of the creatine example and the
normal linear model described in part (a). You may wish to use sketches and/or
make reference to the model parameters.
[6 MARKS]
(e) The following model comparisons were carried out:
anova(m1,m2,m3,m4)
#This compares models sequentially (m4 v m3), (m3 v m2), (m2 v m1)
Models:
m4: Muscle ~ Time
m3: Muscle ~ Time + (1 | Subject)
m2: Muscle ~ Time + (1 | Subject) + (0 + Time | Subject)
m1: Muscle ~ Time + (Time | Subject)
Df
AIC
BIC
logLik deviance
Chisq Chi Df Pr(>Chisq)
m4
3 2041.8 2052.4 -1017.91
2035.8
m3
4 1964.1 1978.2
-978.03
1956.1 79.7474
1
< 2.2e-16 ***
m2
5 1956.1 1973.8
-973.06
1946.1
9.9490
1
0.001609 **
m1
6 1958.0 1979.2
-973.01
1946.0
0.0891
1
0.765376
Based on these results, which model would you select, and why?
[2 MARKS]
CONTINUED OVERLEAF/
3

(f) Selected output from the final correct model is shown below.
Fixed effects:
Estimate Std. Error t value
(Intercept)
49.4249
1.6648
29.688
Time
4.1103
0.4612
8.911
> ranef()
$Subject
(Intercept)
Time
1
1.39676209
1.8205717
2
4.20886279
1.9158831
3
-3.94390531 -2.2711469
4
-6.28764761 -1.1444901
5
12.48438416 -0.2981106
6
7.56926501
1.9019373
7
0.84347632
1.1823382
8
-1.73151497
0.1020777
...
i. Predict the muscle strength of Subject 6 after 4 weeks.
ii. Predict the muscle strength of a new, unobserved subject in week 3.
[4 MARKS]
3. A cancer research charity is interested in investigating some of the factors which de-
termine whether or not cancer treatment will be successful for patients. A researcher
obtained data for 2000 lung cancer patients who had been given treatment in US hos-
pitals - these patients are treated with chemotherapy and have their cancer progress
monitored monthly over a 6-month period. The response variable for the study is
success which is coded as 0 if the patient still has cancer symptoms and 1 if the
patient is in remission (ie has no symptoms of cancer). The dataset also contains a
number of explanatory variables: sex (0=female, 1=male); age - a continuous variable
measuring age in years, tumour - the size of the patient’s cancerous tumour (in cm),
month - the month of the treatment (1-6), id - a unique patient id.
(a) The researcher decides to fit a generalised linear mixed model (GLMM). What
feature of this dataset makes a GLMM a more appropriate choice than a standard
mixed model?
[2 MARKS]
(b) Another alternative model for this type of data would be a generalised estimating
equation (GEE). What are the main differences between a GLMM and a GEE?
[2 MARKS]
CONTINUED OVERLEAF/
4

(c) In a GLMM, the conditional mean relates to the linear predictor through a link
function:
g(E(y|u)) = Xβ + Zu.
What would be an appropriate link function g() for this example? [2 MARKS]
(d) A model was fitted to these data in R. Selected output is shown below.
mod1 <- glmer(success ~ sex + age + tumour + month + ( 1 | id ))
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)
1.52166
1.06991
1.422
0.15496
sexMale
0.66822
0.21086
3.169
0.00153 **
age
0.01449
0.01419
1.021
0.30731
tumour
-1.44060
0.26811
-5.373 7.74e-08 ***
month
0.39948
0.03242
12.322
< 2e-16 ***
Random effects:
Groups Name
Variance Std.Dev.
id
(Intercept) 1.31
1.145
Number of obs: 1046, groups:
id, 356
ranef(mod1)
$id
(Intercept)
1
1.026738316
2
0.045266159
3
0.552361630
4
0.265157840
5
-0.511174643
...
Explain the fixed effects output in terms of each of the covariates in our model.
[6 MARKS]
(e) In a standard mixed model, we would test the significance of the random effect
term by comparing the models with and without the random effect. Why can’t
we do this here?
[2 MARKS]
(f) Subject 1 is a 45 year-old male with a 3cm tumour. What is the probability that
his treatment will have been successful after 5 months?
[6 MARKS]
CONTINUED OVERLEAF/
5

4. A psychologist is interested in how quickly young children can learn to complete cog-
nitive tasks. 50 children were recruited for the study, and she measured how long it
took each child to complete a card matching task. The children were asked to repeat
the task on five consecutive days, to see how much their performance had improved
over time. She was also interested in whether boys and girls performed differently.
The dataset contains the following variables: id - a unique indentifier for each child;
time - the time in seconds for the child to complete the task; day - the day of the
experiment (1-5); sex - the sex of the child (0=female, 1=male).
(a) A model was fitted to these data using the following R code:
mod1 <- gls(time ~ day*sex, data=memory,
correlation = corSymm(form = ~ 1|id),
weights=varIdent(form = ~ 1|day))
Write down the mean model corresponding to this code.
[5 MARKS]
(b) The errors in this model are assumed to follow a multivariate normal distribution
with mean 0 and variance-covariance matrix R. Give the form that R takes when
the following covariance structures are used.
i. unstructured (correlation = corSymm)
ii. AR1 (correlation = corAR1)
[5 MARKS]
(c) The mean model in part (a) was fitted using each of the two covariance structures
above. We would like to carry out a likelihood ratio test to compare these two
structures.
i. Explain why a likelihood ratio test is appropriate for comparing these struc-
tures.
[2 MARKS]
ii. Explain why the asymptotic reference distribution is appropriate for this like-
lihood ratio test.
[2 MARKS]
iii. What is the correct reference distribution for this likelihood ratio test?
[2 MARKS]
CONTINUED OVERLEAF/
6

iv. A likelihood ratio test was carried out, and we obtained the results below.
Note that mod1 is unstructured and mod2 is AR1.
anova(mod1, mod2)
Model df
AIC
BIC
logLik
Test L.Ratio p-value
mod1
1 19 1897.502 1964.104 -929.7512
mod2
2
6 1901.268 1922.300 -944.6341 1 vs 2 29.7659
0.0051
Explain what the result of this likelihood ratio test means in terms of our
covariance structure.
[2 MARKS]
v. The AIC and the BIC appear to disagree with each other in terms of the
preferred model. Explain why this might occur.
[2 MARKS]
END OF QUESTION PAPER.
7





Document Outline