STAT 720 STATISTICAL MODELLING

Assignment 1 Solutions

R code for this solution in a plain text file is also available separately so students can
run the R code for themselves.

Q. 1 a) First we note that the total sample size n = n1 +n2 is fixed so we are only talking about
how to allocate the subjects to the two treatments. For convenience let us suppose that
n1 = cn and so n2 = n − n1 = (1 − c)n. The allocation is then completely determined
by the value of c.
First assume that σ12 = σ22 and we want to choose c which minimizes the variance of
Y 1 − Y 2.
That variance is
σ12 σ22 σ12 1
 
1
Var(Y 1 − Y 2 ) = + = +
n1 n2 n c 1−c
Since c is the only variable under our control in this situation we need only minimize
the function
1 1
g(c) = +
c 1−c
for c ∈ (0, 1). The first and second derivatives of g are
1 1 2 2
g ′ (c) = − 2
+ g ′′ (c) = +
c (1 − c)2 c 3 (1 − c)3
Now we note that g ′′ (c) > 0 for any c ∈ (0, 1) and so any point in this interval satisfying
g ′ (c) = 0 must minimize g and so minimize the variance.

g ′ (c) = 0 ⇒ c2 = (1 − c)2 ⇒ 1 − 2c = 0 ⇒ c = 0.5

Hence we have that when σ12 = σ22 the allocation which minimizes Var(Y 1 − Y 2 ) is
n1 = n2 = 0.5n [5 marks]

Now suppose that Var(Y 1 − Y 2 ) is minimized at n1 = n2 or, equivalently, c = 0.5. As

before, we have
σ2 σ22
Var(Y 1 − Y 2 ) = 1 +
cn (1 − c)n
but we do not know the relation between σ12 and σ22 . The derivatives of the variance
with respect to c are

d Var(Y 1 − Y 2 ) σ2 σ22 d2 Var(Y 1 − Y 2 ) 2σ12 2σ22

= − 21 + = +
dc c n (1 − c)2 n dc2 c3 n (1 − c)3 n

1
 The second derivative is again positive for all c ∈ (0, 1) and all positive σ12 and σ22 . Thus
the variance is minimized at c = 0.5 if
d Var(Y 1 − Y 2 ) 4σ12 4σ22
= − + = 0
dc c=0 n n

If this derivative equals 0 then we must have that σ12 − σ22 = 0 and so, if c = 0.5
minimizes the variance of the difference in means then it must be true that σ12 = σ22 .
[5 marks]

b) Once again let us set n1 = cn so that

σ2
 
1 k
Var(Y 1 − Y 2 ) = 1 +
n c 1−c
where k is a fixed constant not under our control. As before we need only minimize the
function
1 k
g1 (c) = +
c 1−c
This has derivatives
1 k 2 2k
g1′ (c) = − + g1′′ (c) = +
c2 (1 − c)2 c3 (1 − c)3
Once again the second derivative is always positive since k > 0 and c ∈ (0, 1) so the
variance is minimized when

g1′ (c) = 0 ⇒ kc2 = (1 − c)2 ⇒ (k − 1)c2 + 2c − 1 = 0

This is a quadratic equation in c and two solutions given by

p √
−2 ± 4 + 4(k − 1) 1± k
c = =
2(k − 1) 1−k
We now need to find which of these is the correct solution.
√
If k < 1 then 1 + √k > 1 > 1 − k > 0 which leads to an invalid solution of c > 1. Also,
if k > 1 then 1 +√ k > 0 > 1 − k which results in a negative solution for c. Therefore
taking c = (1 + k)/(1√ − k) results in invalid solutions for any k ̸= 1. Conversely, if
k < 1 then 0 < 1 − k <√1 − k and so this gives us a solution for c ∈ (0, 1). Similarly
if k > 1 then 1 − k < 1 − k < 0 which again results in a solution for c in (0, 1). Hence
the correct value of c is
√
1− k 1 + σσ12 σ1 (σ1 − σ2 ) σ1
c = = 2 = 2 2
=
1−k σ
1 − 22 σ1 − σ2 σ1 + σ2
σ1

which gives us the required value of n1 to be

σ1
n1 = n
σ1 + σ2
[6 marks]

2
 If σ12 > σ22 then we have that the variance is minimized when
σ1 n σ1 n n
n1 = > =
σ1 + σ2 σ1 + σ1 2

Hence σ12 > σ22 ⇒ n1 > 0.5n ⇒ n1 > n2 .

Conversely if we the optimal allocation has n1 > n2 then it must be true that σ1 /(σ1 +
σ2 ) > 0.5 which implies that σ1 > σ2 and so σ12 > σ22 . [4 marks]

c) To conduct the pilot study I would take a sample of np = 100 individuals (this is
arbitrary but should be large enough that n0 /2 will give a good estimate of a variance)
and randomly allocate n1p = n2p = 50 to each of the two treatments.
I would then measure the response variable for each of these individuals and calculate
the sample standard deviations s1 and s2 .
A reasonable value of n1 would then be
1000s1
n1 =
s1 + s2
.
In the main study, I would get another n − np = 900 observations. Provided that
50 < n1 < 950, I would round n1 to the nearest integer and randomize n1 − 50 to
treatment A and 950 − n1 to treatment B. If the calculated n1 < 50 then I would assign
all of the new observations to treatment B and, similarly, if the calculated n1 > 950
then I would assign all new observations to treatment A. I would then combine the
new observations with the n0 from the pilot study and analyze these two groups for
differences in means. [5 marks]

3
Q. 2 a) To answer this question we need to examine any interaction between sex and genotype
in their effect on systolic blood pressure. The emmeans package (Lenth, 2025) has
a function [Link] although it is not difficult to use base R functions to
construct such a plot as I do below. Either solution is perfectly acceptable.

BP = [Link]("Data/[Link]", stringsAsFactors=TRUE)
[Link] = with(BP,tapply(SBP, list(Genotype, Sex), mean))
plot([Link]~col([Link]), xlab="Sex", pch=16,
ylab="Mean SBP", xaxt="n", xlim=c(0.8,2),
col=c("Red","Green", "Blue"))
axis(1,c(1,2),levels(BP$Sex))
segments(1,[Link][,1], 2, [Link][,2],
col=c("Red","Green", "Blue"),lty=1, pch=16)
legend("topleft", legend=levels(BP$Genotype),
col=c("Red","Green", "Blue"),
lty=1, pch=16)
122

aa
Aa
AA
Mean SBP

118
114

F M

Sex

The three lines in this plot are very close to perfectly parallel so it looks like the ef-
fect of genotype on systolic blood pressure is the same for males and females. [4 marks]

A formal test of that can be done by examining the interaction term in a full model for
this data.

[Link] = lm(SBP~Sex*Genotype, data=BP)

drop1([Link], test="F")

## Single term deletions

##

4
 ## Model:
## SBP ~ Sex * Genotype
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 620.85 102.897
## Sex:Genotype 2 10.743 631.59 99.411 0.2076 0.8139

From this we see that the test statistic for the interaction term is F = 0.2076 which cor-
responds to a p-value of p = 0.8139. Thus we conclude that the data gives no evidence
at all of a different effect of genotype on systolic blood pressure for different biological
sexes confirming what we say in the plot above. [3 marks]

b) Since we have already seen that the main effect of these two factors is the same regardless
of the level of the other a main effects model is appropriate to examine the effects of
the factors on systolic blood pressure. Since the data is unbalanced (this is unavoidable
here due to natural genotype distributions) we should not use the anova function which
calculates sequential (Type I) sums of squares. In an unbalanced design, it is best to
use the drop1 function to conduct the tests since that function uses hierarchical (Type
II) sums of squares which are the correct ones for the tests of main effects conditional
on the other effects.

[Link] = lm(SBP~Sex+Genotype, data=BP)

drop1([Link], test="F")

## Single term deletions

##
## Model:
## SBP ~ Sex + Genotype
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 631.59 99.411
## Sex 1 229.69 861.29 106.717 9.4555 0.004902 **
## Genotype 2 117.06 748.65 100.512 2.4094 0.109663
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From this we see that the test statistic for the effect of sex on systolic blood pressure is
F = 9.4555 which corresponds to a p-value of p = 0.0049. Thus we conclude that there
is fairly strong evidence of different mean systolic blood pressure in the two biological
sexes.
For genotype we get a test statistic of F = 2.4094 which corresponds to a p-value of
p = 0.1097. Thus the data does not give any compelling evidence that the mean sys-
tolic blood pressure is different for different genotypes at this genetic marker. [6 marks]

c) Age is a continuous covariate so the correct model is an Analysis of Covariance model.

That, however, just means that we include the covariate Age in our linear model.

[Link] = lm(SBP~Age+Sex+Genotype, data=BP)

drop1([Link], test="F")

5
 ## Single term deletions
##
## Model:
## SBP ~ Age + Sex + Genotype
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 461.79 92.017
## Age 1 169.81 631.59 99.411 9.1929 0.0055908 **
## Sex 1 269.25 731.03 103.798 14.5763 0.0007897 ***
## Genotype 2 167.39 629.18 97.296 4.5311 0.0209340 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In this model we see that, as expected, age is related to systolic blood pressure but
adjusting for age has also had an impact on the effect of the other two factors. Sex is
still highly significant but the p-value is now much smaller at p = 7.9 × 10−4 . When
we adjust for Age, however, we see that the genotype term is now also significant with
a p-value of p = 0.020934. Adjusting for the known effect of age on systolic blood
pressure has allowed us to find evidence for this effect which was not apparent when
we did not adjust for age. [6 marks]

d) There a number of assumptions we need to test for the model conducted in part (c).
First we shall examine the relationship between age and systolic blood pressure which
should be linear with the same slope at all levels of the other two factors. To examine
this we look at a plot of age against blood pressure differentiating for the different
treatment levels as well as a formal test for parallel lines. Since we are not including
an interaction between Sex and Genotype, we should also not include the three-way
interaction but only interactions between Age and each of the categorical covariates.

[Link] = lm(SBP~Age*(Sex+Genotype), data=BP)

anova([Link], [Link])

## Analysis of Variance Table

##
## Model 1: SBP ~ Age + Sex + Genotype
## Model 2: SBP ~ Age * (Sex + Genotype)
## [Link] RSS Df Sum of Sq F Pr(>F)
## 1 25 461.79
## 2 22 458.66 3 3.1227 0.0499 0.9849

plot(SBP~Age, data=BP, pch=c(15,17)[Sex],

col=c("red","green","blue")[Genotype],
xlim=range(Age)-c(5,0))
legend("topleft",
c("F aa", "F Aa", "F AA", "M aa", "M Aa", "M AA"),
cex=0.75, bty="n", pch=rep(c(15,17),each=3),
col=rep(c("red","green","blue"),2))

6
 F aa
F Aa

120
F AA
SBP

M aa
M Aa
M AA
105

15 20 25 30 35 40 45 50

Age

From the plot there is no major evidence of non-linearity in the relationship between
age and systolic blood pressure at the different levels of gender (plotting character) or
genotype (colour) although it must be said that low numbers of the aa genotype as well
as females with the Aa genotype mean that we cannot really tell for sure. The test
for unequal slopes is also rejected as the residual sum of squares almost stays the same
even when we allow for unequal slopes. [2 marks]

The other assumptions we make are to do with the error terms in the model. To examine
these we create the standardized residuals and plot them against the fitted values, the
treatment factors and also a normal quantile-quantile plot.

BP$fitted = fitted([Link])
BP$Hat = [Link]([Link])$hat
[Link] = summary([Link])$sigma
BP$sresid = residuals([Link])/([Link]*sqrt(1-BP$Hat))
par(mfrow=c(2,2))
BP$sresid = residuals([Link])/summary([Link])$sigma
plot(sresid~[Link](Genotype), data=BP, pch=16,
ylab="Std. Residuals", xlab="Genotype", xaxt="n")
axis(1,1:3,levels(BP$Genotype))
abline(h=0,lty=2,col="red")
plot(sresid~[Link](Sex), data=BP, pch=16,
ylab="Std. Residuals", xlab="Sex", xaxt="n")
axis(1,1:2,levels(BP$Sex))
abline(h=0,lty=2,col="red")
plot(sresid~fitted, data=BP, pch=16,
ylab="Std. Residuals", xlab="Fitted Values")
abline(h=0,lty=2,col="red")
qqnorm(BP$sresid, pch=16)
abline(0,1, col="red")

7
 2

2
Std. Residuals

Std. Residuals
1

1
0

0
−1

−1
aa Aa AA F M

Genotype Sex

Normal Q−Q Plot

2

2
Sample Quantiles
Std. Residuals

1
0

0
−1

−1

110 115 120 125 −2 −1 0 1 2

Fitted Values Theoretical Quantiles

The first two plots show no evidence of violations of the assumption of equal variances
since the vertical spreads are roughly constant in each plot. There is no evidence of a
pattern in the plot of residuals against fitted values showing the model appears to be
appropriate. Finally the normal quantile plot is remarkably linear so there is also no
indication of a violation of the normality assumption. Overall all of the assumptions
on the ANCOVA model seem to be plausible based on these diagnostics. [4 marks]

8
Q. 3 a) We know from Page 2-10 of my notes that the least squares estimator of β satisfies

−1 2 
t
β̂ ∼ Normalp β, X X σ

Furthermore we can use standard results about the multivariate normal that were given
in class to see that

−1 t 2 
Rβ̂ ∼ Normalq Rβ, R X t X Rσ

since R is a matrix of constants.

If H0 is true then Rβ = 0 and so, when H0 is true we have

−1 t 2 
t
Rβ̂ ∼ Normalq 0, R X X Rσ

[5 marks]

b) Using the distribution of Rβ̂ when H0 is true we have the following quadratic from
1  t 
t

−1 t −1
Q = 2
R β̂ R X X R Rβ̂ ∼ χ2q
σ
[3 marks]

c) From Slide 2-10 of my class notes we know that the least squares estimators are given
by
−1 t
β̂ = X t X XY
and so we can rewrite Q as
1  t

−1 t t  t

−1 t −1 t

−1 t
Q = R X X X Y R X X R R X X XY
σ2
 
t 1 t

−1 t  t

−1 t −1 t

−1 t
= Y X XX R R XX R R XX X Y
σ2

= Y t AY

Where we have
1 t

−1 t  t

−1 t −1
−1 t
A = 2X X X R R XX R R X tX X
σ
[5 marks]

d) We saw in class that the sum of squared residuals can be written as the quadratic form

SSE = Y t (Ip − H) Y

where we define
−1

H = X X tX Xt

9
 Furthermore we know that Var(Y ) = σ 2 In and hence from the statement of the
question SSE is independent of the quadratic form Q if, and only if, A (Ip − H) = 0.
Next we note that
 
1 t

−1 t  t

−1 t −1 t

−1 t h t

−1 t i
AH = X XX R R XX R R XX X X XX X
σ2
1 t

−1 t  t

−1 t −1 t

−1 t t

−1 t
= X X X R R X X R R X X X X X X X
σ2
1 t

−1 t  t

−1 t −1 t

−1 t
= X X X R R X X R R X X X
σ2
= A

and so we have that A (Ip − H) = A − AH = A − A = 0 and hence Q and SSE

are independent quadratic forms. [6 marks]

e) From (b) above we have

1  t 
t

−1 t −1
Q = 2
R β̂ R X X R Rβ̂ ∼ χ2q
σ
We also saw in class that
SSE
∼ χ2n−p
σ2
and from (d) above we have that these two random variables are independent.
Hence from the definition of the F distribution we have
Q/q
F =
SSE / ((n − p)σ 2 )
 t  −1
t −1 t
Rβ̂ R (X X) R Rβ̂/ (qσ 2 )
=
MSE /σ 2
 t  −1
−1
Rβ̂ R (X t X) Rt Rβ̂
= ∼ Fq,n−p
q MSE

To test the null hypothesis of interest (for a fixed known matrix R) we calculate the
observed value (Fobs ) of this quantity and hence we can get p-values as P(Fq,n−p > Fobs ).
[6 marks]
NB: You were not asked to do so but it can be shown that the Wald Test is exactly
equivalent to the Likelihood Ratio test in the normal linear model although they are
different in other situations such as generalized linear models

10
Q. 4 Note: Even though the dataset does include student ages, we do not include it in the model.
The main reason for this is that there is very little variability in the ages since these were
all second-year undergraduate students (mostly 19 or 20) and so it is not a useful predictor.
Although a formal test does agree with this, we are basing the decision to exclude Age based
on prior knowledge of the population and so this does not cause any variable selection issues.

a) Here is the code to fit the model and get a useful summary of that fitted model

Students = [Link]("Data/[Link]",
stringsAsFactors = TRUE)
Stud.model1 = lm(Weight~Height, data=Students)
Stud.summ1 = summary(Stud.model1)
Stud.summ1

##
## Call:
## lm(formula = Weight ~ Height, data = Students)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.209 -15.436 -4.599 9.672 61.271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -85.4990 52.5563 -1.627 0.109
## Height 3.3480 0.7878 4.250 7.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.89 on 58 degrees of freedom
## Multiple R-squared: 0.2375, Adjusted R-squared: 0.2243
## F-statistic: 18.06 on 1 and 58 DF, p-value: 7.861e-05

From the summary output, we see that there is strong evidence that weight depends on
height (p ≈ 7.86 × 10−5 ). We estimate that each 1 inch increase in height is associated
with a 3.35 pound increase in weight on average. From the value of R2 we further see
that taking weight alone into account explains approximately 23.75% of the variability
in height among these students. [2 marks]
For the diagnostic plots, I will first find the standardized residuals and use these in the
plots.

Hii = [Link](Stud.model1)$hat
[Link] = Stud.summ1$sigma
Stud.e = residuals(Stud.model1)/([Link]*sqrt(1-Hii))
par(mfrow=c(3,1))
plot(Weight~Height, data=Students, pch=16, col="blue")
abline(Stud.model1, col="red")
plot(fitted(Stud.model1), Stud.e, xlab="Fitted Values",

11
 ylab="Standarized Residuals", pch=16, col="blue")
abline(h=0, lty=2, col="red")
qqnorm(Stud.e, ylab="Standardized Residuals", main="",
pch=16, col="blue")
abline(0,1, lty=2, col="red")
180
Weight

140
100

55 60 65 70

Height
3
Standarized Residuals

2
1
0
−1

100 110 120 130 140 150 160

Fitted Values
Standardized Residuals

3
2
1
0
−1

−2 −1 0 1 2

Theoretical Quantiles

The first plot shows that there is one point that is quite separate from the rest of
the data in terms of Height. This point also has a very large residual as the person
in question weighed quite a bit more than would be predicted by the model despite
only being 55 inches in height. It is possible that this individual comes from a distinct
population or there was an error in the measurements. In any case we should be very
wary of points like this and look into them further. Another potential issues is that there

12
 appears to be a tendency for the line to overestimate Weight for shorter individuals and
underestimate it for taller individuals suggesting that there may be an issue with the
fitted line. This could be a result of the point identified above or something else (such
as excluding gender).
From the second plot we also see that there may be some issues with heteroscedasticity
since (ignoring the two outliers with large residual) the vertical spread seems to increase
with the fitted values.
The normal quantile plot shows some evidence of skewness with a shorter lower tail and
longer upper tail than would be expected but it does not seem to be a very big issue
here. [6 marks]
b) For this plot we are simply repeating the second plot above but now seeing if gender
may be a cause of the heteroscedasticity and/or the issues with the model.
Note: When doing this part of the question, I realized that I made an error when
switching the numeric coding the original dataset to ”F” and ”M” so these were reversed!
I apologise for that. No marks were deducted from anyone but I gave bonus marks if
you noticed and commented on this issue. I have fixed the dataset on Avenue and used
the corrected dataset here. [2 BONUS marks]

plot(fitted(Stud.model1), Stud.e, xlab="Fitted Values",

ylab="Standarized Residuals", pch=c(15,17)[Students$Sex],
col=c("green","blue")[Students$Sex])
legend("bottomleft", c("Female","Male"), pch=c(15,17),
col=c("green", "blue"), bty="n", cex=0.8)
abline(h=0, lty=2, col="red")
3
Standarized Residuals

2
1
0
−1

Female
Male

100 110 120 130 140 150 160

Fitted Values

We see that there is a definite difference in the residuals by sex. Most of the residuals
for females are negative and most for males are positive. This means that the line

13
 is overestimating weights for females and underestimating weights for males which is
likely a consequence of ignoring sex in the model. We also note that there seems to be
more variability in the residuals for males. [3 marks]
c) In this part of the question we want to have separate intercepts for males and females
which we can achieve by including a dummy variable for Sex in the model. We also
want separate slopes for the two biological sexes which we get by also including an
interaction between Sex and Height.

Stud.model2 = lm(Weight~Height*Sex, data=Students)

beta2 = coef(Stud.model2)
Stud.summ2 = summary(Stud.model2)
Stud.summ2

##
## Call:
## lm(formula = Weight ~ Height * Sex, data = Students)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.827 -12.679 -1.912 10.133 57.896
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -32.650 84.416 -0.387 0.7004
## Height 2.384 1.298 1.837 0.0716 .
## SexM 103.591 106.649 0.971 0.3356
## Height:SexM -1.172 1.610 -0.728 0.4697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.77 on 56 degrees of freedom
## Multiple R-squared: 0.4588, Adjusted R-squared: 0.4299
## F-statistic: 15.83 on 3 and 56 DF, p-value: 1.426e-07

When looking at this output we see that none of the three non-intercept parameters
seem to be significant even though height was highly significant before. This can be
explained by noting that Height and Sex are themselves highly related and recalling
that the parameters are for the added effect of each covariate given the others in the
model. Also, since there is an interaction in the model, we should not really interpret
the main effects. We do, however, not that the value of R2 has almost doubled to
45.88% which suggests that adding Sex to Height does increase the predictive accuracy
of the model for the given dataset. We will formally test this in the next part of the
question. [2 marks]
The two lines fitted by this model are
Height = -32.65 + 2.384 Height for females
Height = 70.94 + 1.213 Height for males
[2 marks]

14
We will now repeat the three plots from part (a) but this time I will use different
plotting symbols for males and females and plot the two different lines.

Hii2 = [Link](Stud.model2)$hat
sigma.hat2 = Stud.summ2$sigma
Stud.e2 = residuals(Stud.model2)/(sigma.hat2*sqrt(1-Hii2))
par(mfrow=c(3,1)) # put 3 plots in a single column
plot(Weight~Height, data=Students,
pch=c(15,17)[Sex], col=c("green","blue")[Sex])
abline(beta2[1:2], col="green")
abline(beta2[1:2]+beta2[3:4], col="blue")
legend("topleft", c("Female","Male"), pch=c(15,17),
col=c("green", "blue"), bty="n", cex=0.8)
plot(fitted(Stud.model2), Stud.e2, xlab="Fitted Values",
ylab="Standarized Residuals",pch=c(15,17)[Students$Sex],
col=c("green","blue")[Students$Sex])
legend("bottomleft", c("Female","Male"), pch=c(15,17),
col=c("green", "blue"), bty="n", cex=0.8)
abline(h=0, lty=2, col="red")
qqnorm(Stud.e2, ylab="Standardized Residuals", main="",
pch=c(15,17)[Students$Sex],
col=c("green","blue")[Students$Sex])
abline(0,1, lty=2, col="red")
legend("topleft", c("Female","Male"), pch=c(15,17),
col=c("green", "blue"), bty="n", cex=0.8)

15
 Female
Male

180
Weight

140
100

55 60 65 70

Height
Standarized Residuals

3
2
1
0

Female
Male
−2

110 120 130 140 150 160

Fitted Values
Standardized Residuals

Female
Male
2
1
0
−2

−2 −1 0 1 2

Theoretical Quantiles

We see from these plots that the issue with the model is now fixed since residuals for
both sexes are evenly distributed. The major outlier from before is now much closer
to the line although may still be of concern due to the unusual value of Height. There
is one new outlier now which is a female whose weight is badly underestimated by the
line for females and seems closer to the males. In practice we would want to examine
this individual record to ensure that the biological sex was properly recorded or if there
are other issues with the datum or the individual in question.
Some issues do still remain with the model. Even though variability seems quite stable
within sex, there does seem to be a difference between the sexes. Also the skewness of
the residuals, although attenuated, is still present. [4 marks]

16
d) The null and alternative hypotheses that we wish to test here are

H0 : β2 = β3 = 0 V H1 : β2 ̸= 0 and/or β3 ̸= 0

where we number the parameters starting at β0 for the intercept.

The easiest way to test this in R is to use the anova function to conduct the likelihood
ratio F-test.

anova(Stud.model1, Stud.model2)

## Analysis of Variance Table

##
## Model 1: Weight ~ Height
## Model 2: Weight ~ Height * Sex
## [Link] RSS Df Sum of Sq F Pr(>F)
## 1 58 27796
## 2 56 19726 2 8069.9 11.455 6.754e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From this we see that the F statistic to compare these two models is quite large and the
p-value is very small so there is strong evidence that including Sex and its interaction
with Height in the model does indeed improve the model. [3 marks]
Note that we could have also used the Wald test from Question 3 which gives the exact
same statistic value and p-value since this is normal linear model. I include R code in
the R file to do that if you wish to check
e) The two lines are parallel if the interaction between Sex and Height (β3 ) equals 0.
Therefore we wish to test

H0 : β3 = 0 V H1 : β3 ̸= 0

We could use the same technique as in part (d) to do this if we wish but it is easiest
to note that, since Sex only has 2 levels, this is exactly the same as the t test already
carried out as part of the summary of the model in part (c). The test statistic is
t = −0.7279 which corresponds to a two-sided p-value of p = 0.4697.
Therefore we have almost a 50% chance of getting an interaction estimate exceeding
1.1717 in absolute value even if the lines are parallel. Based on this we can conclude
that there is little or no evidence against the hypothesis of parallel lines for the two
sexes. We do note that, although these results may not be entirely accurate due to
the issues seen in the plots, they are unlikely to be very far off and so the qualitative
conclusion would remain the same. [3 marks]

17