###############################################

## Joshua M. Tebbs

## Date: 22 May 2010

## Update: 9 August 2011

## STAT 513 course notes: R Code

## Chapter 11

###############################################

 

# Example 11.1.

# Botany data

# Page 69-72

# Enter the data

time=c(16.4, 18.2, 21.6, 22.3, 24.1, 29.7, 34.6, 35.2, 56.5, 58.7, 

    65.5, 68.6, 75.4, 76.3, 88.0, 92.0, 96.6, 98.1, 103.9, 115.9,

    121.6, 121.8, 122.4, 124.4, 128.0, 128.0, 131.2, 140.7, 145.8, 149.5)

absorp=c(5.2, 1.0, 4.8, 2.7, 1.1, 3.5, 8.7, 10.1, 11.4, 10.8,

    15.3, 11.2, 16.9, 12.3, 15.3, 19.9, 21.1, 19.5, 20.7, 22.4,

    23.0, 22.3, 24.6, 22.4, 28.1, 20.5, 26.5, 31.3, 29.1, 32.6)

## Fit the model

fit <- lm(absorp~time)

summary(fit)

plot(time,absorp,xlab="Time (x)",ylab="Absorption rate (y)",pch=16)

abline(fit)

## Confidence/Prediction intervals

predict(fit,data.frame(time=80),level=0.95,interval="confidence")

predict(fit,data.frame(time=80),level=0.95,interval="prediction")

 

# Example 11.2.

# Cheese data

# Page 93-96

# Import the data from a directory that you know

cheese<-read.table("C:\\texfiles\\Classes\\USC\\stat513\\f10\\R\\cheese.txt",header=TRUE)

## Define variables

taste<-cheese[,1]

acetic<-cheese[,2]

h2s<-cheese[,3]

lactic<-cheese[,4]

## Fit the model

fit<-lm(taste~acetic+h2s+lactic)

summary(fit)

## Confidence/Prediction intervals

predict(fit,data.frame(acetic=5.5,h2s=6.0,lactic=1.4),level=0.95,interval="confidence")

predict(fit,data.frame(acetic=5.5,h2s=6.0,lactic=1.4),level=0.95,interval="prediction")

## ANOVA table

# Page 101-102

anova.fit<-anova(lm(taste~acetic+h2s+lactic))

anova.fit

## Reduced versus full model test

# Page 106-107

fit.full<-lm(taste~acetic+h2s+lactic)

fit.reduced<-lm(taste~acetic)

anova(fit.reduced,fit.full,test="F")

 

################################################################################

# Example 11.2.

# Cheese data

##############

## SAS code ##

##############

# This is how you could code this example in SAS

# I use SAS to get the "conventional ANOVA table"

options ps=60 ls=80 pageno=1 formdlim='_' nodate;

data cheese;

input taste acetic h2s lactic;

cards;

12.3 4.543 3.135 0.86

20.9 5.159 5.043 1.53

39.0 5.366 5.438 1.57

47.9 5.759 7.496 1.81

5.6  4.663 3.807 0.99

25.9 5.697 7.601 1.09

37.3 5.892 8.726 1.29

21.9 6.078 7.966 1.78

18.1 4.898 3.850 1.29

21.0 5.242 4.174 1.58

34.9 5.740 6.142 1.68

57.2 6.446 7.908 1.90

0.7  4.477 2.996 1.06

25.9 5.236 4.942 1.30

54.9 6.151 6.752 1.52

40.9 6.365 9.588 1.74

15.9 4.787 3.912 1.16

6.4  5.412 4.700 1.49

18.0 5.247 6.174 1.63

38.9 5.438 9.064 1.99

14.0 4.564 4.949 1.15

15.2 5.298 5.220 1.33

32.0 5.455 9.242 1.44

56.7 5.855 10.20 2.01

16.8 5.366 3.664 1.31

11.6 6.043 3.219 1.46

26.5 6.458 6.962 1.72

0.7  5.328 3.912 1.25

13.4 5.802 6.685 1.08

5.5  6.176 4.787 1.25

;

/* xpx = X'X matrix

i = inverse of X'X

covb = MSE * inverse of X'X (i.e., the estimated var-cov. matrix of beta-hat)

clb = gives confidence intervals for each regression parameter */

 

/* FULL model */

proc reg;

model taste = acetic h2s lactic/xpx i covb clb;

run;

 

/* REDUCED model */

proc reg;

model taste = acetic;

run;

#################################################################################