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(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 ...
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.
= (-200:200)/100
Lambda= exp(mean(log(Wool$Cycles)))
DotY =Lambda
RSS =Lambda
RSS2for(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))}
<- anova(aov(NewResp~Length+Amplitude+Load, data=Wool))["Residuals", "Sum Sq"]
RSS[i] <- anova(aov(NewResp~Length*Amplitude*Load-Length:Amplitude:Load, data=Wool))["Residuals", "Sum Sq"]
RSS2[i]
}201] RSS[
[1] 251900.4
201] RSS2[
[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))
= (-150:100)/1000
Lambda=Lambda
RSS =Lambda
RSS2for(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))}
<- anova(aov(NewResp~Length+Amplitude+Load, data=Wool))["Residuals", "Sum Sq"]
RSS[i] <- anova(aov(NewResp~Length*Amplitude*Load-Length:Amplitude:Load, data=Wool))["Residuals", "Sum Sq"]
RSS2[i] }
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))
min(RSS)
[1] 243413.3
which.min(RSS)
[1] 92
min(RSS2)
[1] 226742.9
which.min(RSS2)
[1] 141