PBP1Anovas.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 <- "PBPRTcompare2runs.txt" dset <- read.table(file.path(folder,filename), header=T, sep="\t", as.is=T)

dset$Time <- factor(dset$Time) dset$Treatment <- factor(dset$Treatment, levels = sort(unique(dset$Treatment)),

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

dset$Run <- factor(dset[,1], levels = sort(unique(dset$Run)), labels=c("Run 2","Run 1")) dset$Insect <- factor(rep(1:(nrow(dset)/3), each=3)) dset$Well <- factor(dset$Well)

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

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. =============================== ANOVA ========================================
2. Response: Quantity, RE
3. Biological run: Run
4. Explanatory variables; Treatment, Time
5. Run:Treatment:Time = Each Insect technical replicate

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(aovTable, file="/Volumes/HD2/Clinton/Latex/PBP11table.tex")
 invisible()

}

aovtable(obj)

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

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

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

1. =========================== Tukey LSD calculation ============================
2. First principles

PBPtukeylsd <- qtukey(0.975, 4, resid.df) * sqrt(2 * MSE / resid.df)

1. Not working!!!
2. TukeyHSD(obj$Insect, "Time:Treatment", ordered=F)
3. library(multcomp)
4. simint.formula(formula = RE ~ Time * Treatment + Error(Run + Insect), data = dset,
5. whichf = "time", type = "Tukey")
6. ==============================================================================

1. diag(obj, "(Intercept)")

diag(obj, "Insect") diag(obj, "Within")