Practical Computing Exercise for Week 2 :Reworking of the Deer Jaw example Solutions

Aims of this practical exercise

In this exercise you will:

rework some aspects of the example given in ELMER, and the associated exercise,
get some practice using R markdown

Before you undertake this exercise…

You need to have installed R, RStudio, and the necessary packages for the course, including the ELMER package. See how to get set up for this course

data(Jaw, package="ELMER")
str(Jaw)

'data.frame':   10 obs. of  5 variables:
 $ Age      : Factor w/ 10 levels "<1","1+","2+",..: 1 2 3 4 5 6 7 8 9 10
 $ Midage   : num  0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5
 $ Number   : int  71 250 210 59 44 34 12 9 8 7
 $ JawLength: num  186 232 255 272 272 ...
 $ SD       : num  21.2 16.6 15.3 11.1 12.9 ...

Fit all possible models

Replicate the model fitting done by Brook & Arnold (1985) as shown in ELMER. Make sure you fit jaw length to the log of age, with and without using weights, and provide a graph to show the model.

Hint: Use weight=Number/SD^2 to get the weighted model.

    Jaw.lm1U <- lm(JawLength~log(Midage), data=Jaw)
    Jaw.lm1W <- lm(JawLength~log(Midage), weight=Number/SD^2, 
data=Jaw)
    plot(JawLength~log(Midage), data=Jaw, ylab="Jaw Length", xlab="ln(Age)", pch=19)
    abline(Jaw.lm1U, col=3)
    abline(Jaw.lm1W, lty=3)
    legend("bottomright",c("Unweighted","Weighted"), lty=c(1,3), col=c(3,1))

Model fitted to the deer jaw data by Brook & Arnold (1985)

Exercise

Use the value k=0.67 to fit the weighted and unweighted models taking the form:

\[y = A - B e^{-kt} \]

where t is the age of the deer.

    summary(Jaw.lmU <- lm(JawLength~exp(-0.67*Midage), data=Jaw))


Call:
lm(formula = JawLength ~ exp(-0.67 * Midage), data = Jaw)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.9741 -1.3389  0.2505  3.1821  5.1176 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          281.505      1.665  169.08 1.67e-15 ***
exp(-0.67 * Midage) -134.125      6.323  -21.21 2.57e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.377 on 8 degrees of freedom
Multiple R-squared:  0.9825,    Adjusted R-squared:  0.9803 
F-statistic: 449.9 on 1 and 8 DF,  p-value: 2.565e-08

    summary(Jaw.lmW <- lm(JawLength~exp(-0.67*Midage), weight=Number/SD^2, data=Jaw))


Call:
lm(formula = JawLength ~ exp(-0.67 * Midage), data = Jaw, weights = Number/SD^2)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-2.06816 -1.16613 -0.01042  0.75952  2.35962 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          282.275      1.318  214.21 2.52e-16 ***
exp(-0.67 * Midage) -137.499      4.984  -27.59 3.22e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.59 on 8 degrees of freedom
Multiple R-squared:  0.9896,    Adjusted R-squared:  0.9883 
F-statistic:   761 on 1 and 8 DF,  p-value: 3.216e-09

Use these models to predict jaw length for the range of ages and plot JawLength against Age with the fitted curves added.

Jaw |>    ggplot(aes(y=JawLength, x=Midage))  + ylab("Jaw Length") + xlab("Age") + geom_point() + 
geom_function(fun=function(x) coef(Jaw.lmU)[1] + coef(Jaw.lmU)[2]*exp(-0.67*x), col="red") + 
geom_function(fun=function(x) coef(Jaw.lmW)[1] + coef(Jaw.lmW)[2]*exp(-0.67*x), col="blue")

In spite of the fact that the coefficients for these models look different, the practical outcome is that the lines on the plot are pretty similar; so much so, that we might decide to use the simpler unweighted model in this instance.

Optional extra task

Vary the value of k to see if the model can be made to fit the data better.

    summary(Jaw.lmW2 <- lm(JawLength~exp(-0.6*Midage), weight=Number/SD^2, data=Jaw))


Call:
lm(formula = JawLength ~ exp(-0.6 * Midage), data = Jaw, weights = Number/SD^2)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.2842 -0.8694 -0.5043  0.3828  2.6419 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         284.788      1.354  210.40 2.91e-16 ***
exp(-0.6 * Midage) -131.463      4.647  -28.29 2.63e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.551 on 8 degrees of freedom
Multiple R-squared:  0.9901,    Adjusted R-squared:  0.9889 
F-statistic: 800.3 on 1 and 8 DF,  p-value: 2.634e-09

    summary(Jaw.lmW3 <- lm(JawLength~exp(-0.75*Midage), weight=Number/SD^2, data=Jaw))


Call:
lm(formula = JawLength ~ exp(-0.75 * Midage), data = Jaw, weights = Number/SD^2)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.6743 -1.0749  0.5819  1.9320  2.9640 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          279.810      1.595  175.38 1.25e-15 ***
exp(-0.75 * Midage) -144.672      6.696  -21.61 2.22e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.024 on 8 degrees of freedom
Multiple R-squared:  0.9831,    Adjusted R-squared:  0.981 
F-statistic: 466.8 on 1 and 8 DF,  p-value: 2.22e-08

The selection of k has had little impact really.

Massey University — School of Mathematical and Computational Sciences

161.331 Biostatistics