### Additional R code for Chapter 8:


### Special Models:

##### Polynomial regression on the Rabbit Jawbone Data (Table 8.8)


##########
##
# Reading the data into R:

my.datafile <- tempfile()
cat(file=my.datafile, "
0.01 15.5
0.2 26.1
0.2 26.3
0.21 26.7
0.23 27.5
0.24 27
0.24 27
0.25 26
0.26 28.6
0.34 29.8
0.41 29.7
0.83 37.7
1.09 41.5
1.17 41.9
1.39 48.9
1.53 45.4
1.74 48.3
2.01 50.7
2.12 50.6
2.29 49.2
2.52 49
2.61 45.9
2.64 49.8
2.87 49.4
3.39 51.4
3.41 49.7
3.52 49.8
3.65 49.9
", sep=" ")

options(scipen=999) # suppressing scientific notation

rabbit <- read.table(my.datafile, header=FALSE, col.names=c("age","length")) 
  
attach(rabbit)

#########

## Initial scatterplot of data:

plot(age, length)

# Not a linear trend?

# Defining transformed variables for the polynomial terms:

age2 <- age^2
age3 <- age^3
age4 <- age^4

quartic.reg <- lm(length ~ age + age2 + age3 + age4)
summary(quartic.reg)

## Looks like the fourth-degree term is not needed.

cubic.reg <- lm(length ~ age + age2 + age3)
summary(cubic.reg)

# Plotting the fitted curve (this plotting approach only works well with a moderate to large number of observations)

plot(age, length)
lines(age, fitted(cubic.reg))  


##### Multiplicative model with the squid data (Table 8.10):

my.datafile <- tempfile()
cat(file=my.datafile, "
1 1.31 1.07 0.44 0.75 0.35 1.95
2 1.55 1.49 0.53 0.9 0.47 2.9
3 0.99 0.84 0.34 0.57 0.32 0.72
4 0.99 0.83 0.34 0.54 0.27 0.81
5 1.05 0.9 0.36 0.64 0.3 1.09
6 1.09 0.93 0.42 0.61 0.31 1.22
7 1.08 0.9 0.4 0.51 0.31 1.02
8 1.27 1.08 0.44 0.77 0.34 1.93
9 0.99 0.85 0.36 0.56 0.29 0.64
10 1.34 1.13 0.45 0.77 0.37 2.08
11 1.3 1.1 0.45 0.76 0.38 1.98
12 1.33 1.1 0.48 0.77 0.38 1.9
13 1.86 1.47 0.6 1.01 0.65 8.56
14 1.58 1.34 0.52 0.95 0.5 4.49
15 1.97 1.59 0.67 1.2 0.59 8.49
16 1.8 1.56 0.66 1.02 0.59 6.17
17 1.75 1.58 0.63 1.09 0.59 7.54
18 1.72 1.43 0.64 1.02 0.63 6.36
19 1.68 1.57 0.72 0.96 0.68 7.63
20 1.75 1.59 0.68 1.08 0.62 7.78
21 2.19 1.86 0.75 1.24 0.72 10.15
22 1.73 1.67 0.64 1.14 0.55 6.88
", sep=" ")

options(scipen=999) # suppressing scientific notation

squid <- read.table(my.datafile, header=FALSE, col.names=c("obs", "rl", "wl", "rnl", "nwl", "w", "wt")) 
  
attach(squid)

#########

# We will only use weight, rostal length, and width in our model:

# Untransformed regression model:

untran.reg <- lm(wt ~ rl + w)
summary(untran.reg)

plot(fitted(untran.reg), resid(untran.reg), xlab="Fitted Values", ylab="Residuals"); abline(h=0)

# Residual plot shows some violations of model assumptions.

####

# Defining transformed variables for the logarithmic terms:

log.wt <- log(wt)  
log.rl <- log(rl)  
log.w <- log(w)

# transformed regression model:

mult.reg <- lm(log.wt ~ log.rl + log.w)
summary(mult.reg)

plot(fitted(mult.reg), resid(mult.reg), xlab="Fitted Values", ylab="Residuals"); abline(h=0)

# Much better!