Practical Exercise: The Wool exercise - Solutions

This exercise is useful for the following courses:

Aims of this practical exercise

In this exercise you will:

  • rework the example given in ELMER
  • get some practice using R markdown
  • fit a variety of models to the Wool data involving a variety of transformations.

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

Get the data

data(Wool, package="ELMER")
str(Wool)
'data.frame':   27 obs. of  5 variables:
 $ X        : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Length   : int  -1 -1 -1 -1 -1 -1 -1 -1 -1 0 ...
 $ Amplitude: int  -1 -1 -1 0 0 0 1 1 1 -1 ...
 $ Load     : int  -1 0 1 -1 0 1 -1 0 1 -1 ...
 $ Cycles   : int  674 370 292 338 266 210 170 118 90 1414 ...

The Exercise

Re-create the graph for the residual sum of squares resulting from fitting first and second order models based on differing selections of \(\lambda\). N.B. To make a comparison using RSS, you will need to use the scale invariant form of the transformation.

Refine the search space for \(\lambda\), and evaluate the benefit of choosing a convenient value for \(\lambda\) as against the optimal value.

    Lambda= (-200:200)/100
    DotY = exp(mean(log(Wool$Cycles)))
    RSS =Lambda
    RSS2=Lambda
    for(i in 1:length(Lambda)){
    if(Lambda[i]==0){NewResp= DotY*log(Wool$Cycles)}
    else{NewResp=(Wool$Cycles^Lambda[i] -1)/(Lambda[i]*DotY^(Lambda[i]-1))}
    RSS[i] <- anova(aov(NewResp~Length+Amplitude+Load, data=Wool))["Residuals", "Sum Sq"]
    RSS2[i] <- anova(aov(NewResp~Length*Amplitude*Load-Length:Amplitude:Load, data=Wool))["Residuals", "Sum Sq"]
    }
    RSS[201]
[1] 251900.4
    RSS2[201]
[1] 226814.5
    min(RSS)
[1] 243414.4
    which.min(RSS)
[1] 195
    min(RSS2)
[1] 226742.9
    which.min(RSS2)
[1] 200
    plot(RSS~Lambda, xlab=expression(lambda), ylim=range(c(RSS,RSS2)), ylab="Residual SS", type="l", lty=2)
    lines(Lambda, RSS2, col=2)
    legend("topleft",  c("First order model","Second order model"), lty=c(2,1), col=c(1,2))

to fix

    Lambda= (-150:100)/1000
    RSS =Lambda
    RSS2=Lambda
    for(i in 1:length(Lambda)){
    if(Lambda[i]==0){NewResp= DotY*log(Wool$Cycles)}
    else{NewResp=(Wool$Cycles^Lambda[i] -1)/(Lambda[i]*DotY^(Lambda[i]-1))}
    RSS[i] <- anova(aov(NewResp~Length+Amplitude+Load, data=Wool))["Residuals", "Sum Sq"]
    RSS2[i] <- anova(aov(NewResp~Length*Amplitude*Load-Length:Amplitude:Load, data=Wool))["Residuals", "Sum Sq"]
    }
    plot(RSS~Lambda, xlab=expression(lambda), ylim=range(c(RSS,RSS2)), ylab="Residual SS", type="l", lty=2)
    lines(Lambda, RSS2, col=2)
    legend("topleft",  c("First order model","Second order model"), lty=c(2,1), col=c(1,2))

to fix

    min(RSS)
[1] 243413.3
    which.min(RSS)
[1] 92
    min(RSS2)
[1] 226742.9
    which.min(RSS2)
[1] 141