Power Transformation

Toxic Agents Example

> # poisons data set is in package boot

> require(boot) # load package boot

Loading required package: boot
[1] TRUE

> poisons[1:3,] # see first 3 obs.

  time poison treat
1 0.31      1     A
2 0.45      1     A
3 0.46      1     A

> # Both factors are in R factor mode

> is.factor(poisons$poison)

[1] TRUE

> is.factor(poisons$treat)

[1] TRUE

> #

> # Produce full linear model (lm) object and residual

> # plot.

> poisons.lm <− lm(time ~ poison*treat, data=poisons)

> windows(width=4,height=3,pointsize=10)

> par(mar=c(4,4,.5,.5)+.1)

> plot(fitted(poisons.lm), residuals(poisons.lm),

+      pch=21, bg="green3")

> abline(h=0, lty="dotted", col="red")

> dev.off()

null device 
          1 

(poisons01.png here)

> #

> # Now produce a residual plot of cell means modeled

> # against main effects. Need to make a data frame

> # of cell means.

> poisons.mean <− with(poisons,

+     tapply(time,list(poison,treat),mean))

> poisons.mean

       A     B      C      D
1 0.4125 0.880 0.5675 0.6100
2 0.3200 0.815 0.3750 0.6675
3 0.2100 0.335 0.2350 0.3250

> row(poisons.mean)

     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    2    2    2    2
[3,]    3    3    3    3

> factor(row(poisons.mean))

 [1] 1 2 3 1 2 3 1 2 3 1 2 3
Levels: 1 2 3

> col(poisons.mean)

     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4

> factor(LETTERS[col(poisons.mean)])

 [1] A A A B B B C C C D D D
Levels: A B C D

> poisons2 <− data.frame( cells=poisons.mean[1:12],

+   poison=factor(row(poisons.mean)),

+   treat=factor(LETTERS[col(poisons.mean)]) )

> poisons2

    cells poison treat
1  0.4125      1     A
2  0.3200      2     A
3  0.2100      3     A
4  0.8800      1     B
5  0.8150      2     B
6  0.3350      3     B
7  0.5675      1     C
8  0.3750      2     C
9  0.2350      3     C
10 0.6100      1     D
11 0.6675      2     D
12 0.3250      3     D

> poisons2.lm <− lm(cells ~ poison + treat, data=poisons2)

> windows(width=4, height=3, pointsize=10)

> par(mar=c(4,4,.5,.5)+.1)

> plot(fitted(poisons2.lm), residuals(poisons2.lm),

+      pch=21, bg="green3")

> abline(h=0, lty="dotted", col="red")

> dev.off()

null device 
          1 

(poisons02.png here)

> #

> # To decide the power, informal method 1 can be used:

> # Now, produce a plot of ln(σ) vs. ln(μ).

> # A graphical function logslplot can be downloaded from

> #

> # http://www.stat.wmich.edu/wang/R/codes/logslplot.R

> #

> windows(width=4,height=3,pointsize=10)

> logslplot(time ~ poison + treat, data=poisons,

+     pch=21, bg="green3")

   slope 
1.977040

> dev.off()

null device 
          1

(poisons03.png here)

> #

> # Use informal method 2:

> # A graphical function logrfplot can be downloaded from

> #

> # http://www.stat.wmich.edu/wang/R/codes/logrfplot.R

> #

> windows(width=4,height=3,pointsize=10)

> par(mar=c(4,4,4,.5)+.1)

> logrfplot(poisons.lm, pch=21, bg="green3")

   theta 
1.930791 

> dev.off()

null device 
          1

(poisons04.png here)

> #

> # Now use Box-Cox method:

> # A graphical function boxcox can be downloaded from

> #

> # http://www.stat.wmich.edu/wang/R/codes/boxcox.R

> #

> windows(width=4,height=3,pointsize=10)

> par(mar=c(4,4,.5,.5)+.1)

> with(poisons, boxcox(time,poison,treat,

+                      lambda=seq(-1.5,0,.01),identify=T))

Identify confidence limits!

> # Upon submitting above graphical function call with

> # identify. specified as TRUE, you are prompted to

> # identify confidence limit. Click the graphics window

> # to make it current and left click mouse on the two

> # intersecting points of the SSλ curve and the (dotted)

> # horizontal line marking SS.

> dev.off()

null device 
          1

(poisons05.png here)

> # The optimal value for λ is 0.75. As commented in

> # class notes, a practical value of power λ = -1

> # (reciprocal transformation) is considered which is

> # consistent with the results from the above two informal

> # methods.

> # A better graphical function boxcox2 can be downloaded from

> #

> # http://www.stat.wmich.edu/wang/R/codes/boxcox2.R

> #

> # There is no need of interactive identification of the

> # confidence limits. Do this only:

> with(poisons, boxcox2(time,poison,treat,

+                      lambda=seq(-1.5,0,.01),identify=T))

> #

> # With reciprocal transformation, additive model

> # will be appropriate which is evident from inspecting

> # the interaction plots of the transformed data below:

> cols <− paste("dark",c("magenta","green","blue","orange4"),sep="")

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2),mar=c(4,4,.5,3)+.1,oma=c(2,0,3,0),bty="l")

> with(poisons, {

+   interaction.plot(poison,treat,1/time,type="b",

+       col=cols[1:4], pch=1:4)

+   box(which="figure",col="grey80")

+   interaction.plot(treat,poison,1/time,type="b",

+       col=cols[1:3], pch=1:3)

+   box(which="figure",col="grey80")

+ })

> dev.off()

null device 
          1

(poisons06.png here)

> #

> # Now, fit the reciprocal response with additive model:

> poisons3.lm <− lm(1/time ~ poison + treat, data=poisons)

> anova(poisons3.lm)

Analysis of Variance Table

Response: 1/time
          Df Sum Sq Mean Sq F value    Pr(>F)    
poison     2 34.877  17.439  71.708 2.865e-14 ***
treat      3 20.414   6.805  27.982 4.192e-10 ***
Residuals 42 10.214   0.243                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> #

> # Now, the residuals behave nicely:

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2),oma=c(0,0,2,0)+.1)

> plot(poisons3.lm,which=1:2,add.smooth=F)

> dev.off()

null device 
          1

(poisons07.png here)

> #

> # To complete analysis, follow-up studies should be

> # conducted for the two main effects on the transformed

> # data.