help file2 for R-code

get_classical_sd_v3.r

+#
+# get_classical_sd_v1.r
+#
+
+# Can get sds for nested random effects using
+get.pairwise.sd <- function(an.df, balance=FALSE) {
+  
+  if(balance) {
+  	an.df$obs <- with(an.df, ave(Concentration, Subject, FUN=function(x) sum(!is.na(x))) )
+  	an.df <- droplevels( subset(an.df, obs == 4 ) )
+  }
+  session.ag <- 
+    aggregate(an.df$Concentration, by=an.df[c("Subject", "Session")],
+    	FUN=function(x) if(length(x)==2) diff(x) else NA)
+  
+  # This includes subject-specific session IDs, so we need a session 1/2 ID
+  session.ag$Session12 <- factor( gsub( ".:(.)", "\\1", session.ag$Session ) )
+  
+  session.vars <- with(session.ag, tapply(x, Session12, var, na.rm=TRUE))
+
+  session.mean.ag <- 
+    aggregate(an.df$Concentration, by=an.df[c("Subject", "Session")],FUN=mean, na.rm=TRUE)
+  smean.diffs <- with(session.mean.ag, tapply(x, Subject, diff, na.rm=TRUE))
+  smean.diffs <- as.vector(smean.diffs) # Array causes trouble in pipes (?!).
+  var.Sessdiff <- var(smean.diffs, na.rm=TRUE)
+
+ 
+  # Session 1 changes position; session 2 does not
+  var.w <- session.vars["2"]/2
+  var.Pos <- (session.vars["1"] - session.vars["2"])/2
+  var.Sess <- (var.Sessdiff - 3*var.Pos/2 - var.w) /2
+
+  non.neg.vars <- pmax(0, c(var.w, var.Pos, var.Sess))
+  S.R <- sqrt(sum( non.neg.vars ))
+  
+  DF.w=length(na.omit(subset(session.ag, Session12=="2")$x)) -1
+  DF.Pos=length(na.omit(subset(session.ag, Session12=="1")$x)) -1 
+  DF.Sess=length(na.omit(smean.diffs))-1 
+  
+  mean.C <- mean(an.df$Concentration, na.rm=T)
+  
+  data.frame(Mean=mean.C, 
+                S.w=sqrt(var.w),                DF.w=DF.w, 
+                S.Pos = sqrt(max(0,var.Pos)),   DF.Pos=DF.Pos,
+                S.Sess = sqrt(max(0,var.Sess)), DF.Sess=DF.Sess,
+                S.R = S.R, RSD.R = S.R / mean.C)
+}
+
+
+#Also possible from anova using expected MS (see text)
+get.anova.sd <- function(an.df, raw.vars=FALSE) {
+ 
+  # First remove Subjects with missing or partial data
+  # to provide balanced experiment
+  
+  an.df$obs <- with(an.df, ave(Concentration, Subject, FUN=function(x) sum(!is.na(x))) )
+  
+  an.df.x <- subset(an.df, obs == 4 )
+  
+  a <- anova(lm(Concentration~Subject/Session/Position, an.df.x))
+
+  # variances derived from expected mean squares:
+  var.w <- a[4,3]
+  var.Pos <- a[3,3] - a[4,3]
+  var.Sess <- a[2,3]/2 - 3*a[3,3]/4 +  a[4,3] / 4
+  # var.Sess <- a[2,3]/2 - 3*var.Pos / 2 -  var.w
+  
+  non.neg.vars <- pmax(0, c(var.w, var.Pos, var.Sess))
+  S.R <- sqrt(sum( non.neg.vars ))
+  
+  DF.w <- a[4,1] 
+  DF.Pos= a[3,1] 
+  DF.Sess=a[2,1] 
+  
+  mean.C <- mean(an.df$Concentration, na.rm=T)
+  
+  rv <- data.frame(Mean=mean.C, 
+             S.w=sqrt(var.w),                DF.w=DF.w, 
+             S.Pos = sqrt(max(0,var.Pos)),   DF.Pos=DF.Pos,
+             S.Sess = sqrt(max(0,var.Sess)), DF.Sess=DF.Sess,
+             S.R = S.R, RSD.R = S.R / mean.C)
+  if(!raw.vars) {
+  	rv 
+  } else  {
+  	RawVars<-data.frame(Within=var.w, Position=var.Pos, Session=var.Sess)
+  	RawVars <- rbind(RawVars, sqrt(pmax(0, c(var.w, var.Pos, var.Sess))))
+  	list(SD=rv, RawVars=RawVars )
+  }
+}
+
+