In this exercise you will:
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 ...
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.
<- lm(JawLength~log(Midage), data=Jaw)
Jaw.lm1U <- lm(JawLength~log(Midage), weight=Number/SD^2,
Jaw.lm1W 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))
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.
|> ggplot(aes(y=JawLength, x=Midage)) + ylab("Jaw Length") + xlab("Age") + geom_point() +
Jaw 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.
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.