CXE1Anovas.R

Code snipits and programs written in R, S or S-PLUS

1. ==============================================================================
2. Description: Insects from one line are bred at the same time.
3. Catapillars are fed with two different treatments Esterase+Suc, or Suc
4. and then at time = 0,1,3,5,7 they are killed and their RNA frozen at -80
5. Four real time runs are run over different days
6. Three technical replicates (3 Wells) on each 384 well plate for each insect
7. are measured using real time RT-PCR
8. Main variability stratums of variation of concern are;
9. Differences between runs, and differences between insects
10. No effect of well position across the heat block over runs should exist
11. ==============================================================================

1. =============================== Setup ========================================

folder <- "/Volumes/HD2/Clinton/Data" filename <- "esteraseRTcompare4runs2.txt" dset <- read.table(file.path(folder,filename), header=T, sep="\t", as.is=T) dset$Time <- factor(dset$Time / 24)

dset$Treatment <- factor(dset$Treatment, levels=c("Sucrose","Esterase"),

                                             labels=c("Sucrose","CXE1 dsRNA"))

dset$Run <- factor(dset$Run, levels=levels(factor(dset$Run)),

                                   labels=c("Run 4","Run 3","Run 2","Run 1"))

dset$Insect <- factor(rep(1:(nrow(dset)/3), each=3)) dset$Well <- factor(dset$Well)

1. ============================ Normalization ===================================

1. Normalization of realtime RT-PCR

dset$RE <- dset$Quantity/dset$NormFactor

1. ==============================================================================

diag <- function(x, label) {

 par(mfrow=c(2,2))

 diag.plots <- function(x) {
   a <- plot(x$fit, x$res, xlab="fitted", ylab="residuals", main="residuals vs fitted")
   print(a)
   abline(h=0)
   qqnorm(x$res)
   
   jack <- rstudent(x)
   plot(jack, ylab="Jacknife residuals", main="Jackknife")
   
   cooks <- cooks.distance(x)
   if(!all(is.na(cooks))) {       
     plot(cooks, ylab="Jacknife residuals")
   }

 }
 if(missing(label)) {
   diag.plots(obj)
 } else {
   diag.plots(objlabel)
 }
   invisible()

}

1. ==============================================================================

1. ============================ Diagnostics =====================================

1. Details

dim(dset) summary(dset) table(dset$Time, dset$Treatment) table(dset$Run, dset$Time, dset$Treatment)

quartz() par(mfrow=c(1,3)) hist(dset$RE) plot(density(dset$RE)) plot(sort(dset$RE), pch=".")

1. =============================== ANOVA ========================================
2. Response: Quantity, RE
3. Biological run: Run
4. Explanatory variables; Treatment, Time
5. Run:Treatment:Time = Each Insect technical replicate

dset$Time <- factor(dset$Time)

1. ===================== Naive ANOVA (incorrect df) =============================

obj <- aov(RE ~ Time * Treatment, data=dset) summary(obj) quartz() par(mfrow=c(2,2))

1. plot(obj)

diag(obj)

1. diagnostics

quartz() par(mfrow=c(1,1))

1. plot(obj)

dset[c(9,19,87),] diag(obj)

1. ============== Testing Run as a fixed and random effect ======================
2. Estimating Run as a fixed effect or Random effect makes no difference

obj <- aov(RE ~ Time * Treatment + Run, data=dset) summary(obj) obj <- aov(RE ~ Time * Treatment + Error(Run), data=dset) summary(obj) diag(obj, "Run")

1. ==============================================================================

1. ================ Testing Run and Run:Time:Treatment =========================
2. Estimating Run as a fixed effect or Random effect makes no difference

obj <- aov(RE ~ Time * Treatment + Run + Error(Run:Time:Treatment), data=dset) summary(obj) obj <- aov(RE ~ Time * Treatment + Error(Run + Run:Time:Treatment), data=dset) summary(obj) diag(obj, "Run") diag(obj, "Run:Time:Treatment")

1. ==============================================================================

1. =================== Run + Run:Treatment + Run:Treatment:Time =================
2. This model does not need Run:Treatment stratum, treatment is a fixed effect

obj <- aov(RE ~ Time * Treatment + Error(Run/Treatment/Time), data=dset) summary(obj) diag(obj, "Run") diag(obj, "Run:Treatment") diag(obj, "Run:Treatment:Time")

1. ==============================================================================

1. ========================= Testing Run and Insect =============================
2. Insect is the same as Run:Treatment:Time

default <- options(show.signif.stars=F) obj <- aov(RE ~ Time * Treatment + Run + Error(Insect), data=dset) summary(obj) obj <- aov(RE ~ Time * Treatment + Error(Run + Insect), data=dset) summary(obj)

summary(obj, split=list(Time=list(L=1, Q=2, C=3)))

options(default)

1. ==============================================================================

aovtable <- function(obj) {

 library(Hmisc)
 Sys.putenv("DISPLAY"=":0")
 anova.header <- t(c("$Source$","$df$","$SS$","$MS$","$F$","$p{\\;}value$"))
 
 Run <- matrix("", nc=6, nr=1)
 tmp <- summary(obj)$"Error: Run"1[1:3]
 Run[,1] <- rownames(tmp[1])
 Run[1,2] <- unlist(tmp[1])
 Run[1,3:4] <- round(unlist(tmp[2:3]),2)
 
 Insect <- matrix("", nc=6, nr=4)
 Insect[,1] <- rownames(summary(obj)$"Error: Insect"1)
 Insect[,2] <- summary(obj)$"Error: Insect"1[,1]
 Insect[,3] <- round(summary(obj)$"Error: Insect"1[,2],2)
 Insect[,4] <- round(summary(obj)$"Error: Insect"1[,3],2)
 Insect[,5] <- round(summary(obj)$"Error: Insect"1[,4],2)
 Insect[,6] <-  signif(summary(obj)$"Error: Insect"1[,5],2)
 
 Replicates <- matrix("", nc=6, nr=1)
 tmp <- summary(obj)$"Error: Within"1[1:3]
 Replicates[,1] <- rownames(tmp[1])
 Replicates[1,2] <- unlist(tmp[1])
 Replicates[1,3:4] <- round(unlist(tmp[2:3]),2)
 
 aovTable <-  rbind(
                    c("$Error(Run)$", rep("",5)), Run,
                    c("$Error(Insect)$", rep("",5)), Insect,
                    c("$Error(Reps)$", rep("",5)), Replicates
                    )
 
 aovTable[is.na(aovTable)] <- ""
 dimnames(aovTable)2 <- anova.header
 latex.default(aovTable, file="/Volumes/HD2/Clinton/Latex/CXE1table.tex")#, rowlabel.just="r")
 invisible()

}

aovtable(obj)

1. latex(aovTable, file="defaultTable")
2. latex(aovTable, title="cTable", ctable=T)
3. latex(aovTable, title="bookTable", booktabs=T)

                                       #options(digits=2)

names(summary(obj)) summary(obj)$"Error: Run"1[1:3] summary(obj)$"Error: Insect"1 summary(obj)$"Error: Within"1[1:3]

latex(anova.header)

1. ============================= LSD calculation ================================
2. Extracting the correct residual

MSE <- summary(obj$Insect)1[4,"Mean Sq"] resid.df <- summary(obj$Insect)1[4,"Df"] CXElsd <- qt(0.995, resid.df) * sqrt(2 * MSE / resid.df)

1. use tukey ??? qtukey(0.975, 5,4) * MSE

1. ==============================================================================

ls()

1. fullobj <- aov(RE ~ Time * Treatment + Run + Run:Time + Run:Treatment + Error(Run:Time:Treatment), data=dset)
2. summary(fullobj)

1. plot(cooks, ylab="Jacknife residuals")
2. identify(2:40, cooks, 2:40)

1. Within doesnt make sence, fit =0

plot(obj$Within$fit, obj$Within$res)

1. first principles pvalue calculation example

1-pf(1.47117 / 0.03194 , 1, 30)

1. There is a significant interaction effect (Time:Treatment), but the main effect of Treatment
2. has the most signal in the F values. The experimental size has probably contributed to the
3. significance in the Time:Treatment interaction.
4. Check assumptions after taking into account Interaction information
5. Residuals normally distributed (Slight departures in normality)