1+ # calculate area under the curve (AUC) for multiple groups
2+ multi.group.AUC <- function (data ,subject.id ,sample.type , time ){
3+ library(pracma )
4+ # too lazy to rename objects from older fxn
5+ subject.id <- as.factor(subject.id )
6+ fact <- as.factor(sample.type ) # sample type factor
7+ tme <- as.factor(time ) # time
8+
9+ # split objects
10+ tmp.data <- split(data.frame (data ),fact )
11+ tmp.time <- split(tme ,fact )
12+ tmp.subs <- split(as.character(subject.id ),fact )
13+
14+ group.AUC <- lapply(1 : nlevels(fact ),function (i ){
15+ ddata <- tmp.data [[i ]]
16+ ttime <- tmp.time [[i ]]
17+ subs <- tmp.subs [[i ]]
18+
19+ # calculate AUC
20+ AUC <- sapply(1 : length(ddata ),function (i )
21+ {
22+
23+ obj <- split(as.data.frame(ddata [[i ]]),subs )
24+ # subtract baseline (first level of time) to correct negative AUC
25+ tmp.time <- split(ttime ,subs )
26+ base.time <- levels(ttime )[1 ]
27+ base.obj <- lapply(1 : length(obj ),function (j )
28+ {
29+ tmp <- as.numeric(as.matrix(unlist(obj [[j ]])))
30+ tmp - tmp [tmp.time [[j ]]== base.time ]
31+ })
32+ tmp <- split(as.data.frame(ttime ),subs )
33+ # x11()
34+ # plot(as.numeric(as.matrix(do.call("cbind",tmp))),as.numeric(as.matrix(do.call("cbind",base.obj))))
35+ out <- as.data.frame(sapply(1 : length(obj ),function (j )
36+ {
37+ x <- as.numeric(as.matrix(unlist(tmp [[j ]])))
38+ o <- order(x ) # need to be in order else AUC will be wrong!
39+ y <- as.numeric(as.matrix(unlist(base.obj [[j ]])))
40+ trapz(x [o ],y [o ])
41+ }))
42+ colnames(out )<- colnames(data [i ])
43+ out
44+ })
45+ tmp <- do.call(" cbind" ,AUC )
46+ rownames(tmp )<- paste(levels(fact )[i ],names(split(as.data.frame(ttime ),subs )),sep = " _" )
47+ tmp
48+ })
49+ res <- do.call(" rbind" ,group.AUC )
50+ colnames(res )<- colnames(data )
51+ return (res )
52+ }
53+
154# function to convert pattern to a single char objects name
255rename <- function (x , pattern , replace = " _" ){
356 # strangely sapply will not work without effort here
@@ -9,7 +62,6 @@ rename <- function(x, pattern, replace="_"){
962 return (x )
1063}
1164
12-
1365# relative standard deviation
1466# redo calc.stat using dplyr
1567calc.rsd <- function (data ,factor ,sig.figs = 2 ){
@@ -134,7 +186,7 @@ anova.formula.list<-function(data,formula,meta.data){
134186}
135187
136188# ANOVA with repeated measures and post hoc
137- aov.formula.list <- function (data ,formula ,meta.data = factor ,post.hoc = TRUE ,repeated = NULL ,p.adjust = " BH" ){
189+ aov.formula.list <- function (data ,formula ,meta.data = factor ,post.hoc = TRUE ,repeated = NULL ,FDR = " BH" ){
138190 # formula = formula excluding repeated terms
139191 # meta data = all factors
140192 # repeated name of factor for repeated measures
@@ -147,6 +199,7 @@ aov.formula.list<-function(data,formula,meta.data=factor,post.hoc=TRUE,repeated=
147199 results <- list (p.value = vector(" list" ,ncol(data )),post.hoc = vector(" list" ,ncol(data )))
148200 for (i in 1 : ncol(data )){
149201 model <- tryCatch(aov(as.formula(paste(" data[," ,i ," ]~" ,formula ,sep = " " )),data = tmp.data ), error = function (e ){NULL })
202+
150203 # get p-values
151204 if (! is.null(repeated )){
152205 names <- attr(model $ Within $ terms ,' term.labels' )
@@ -157,6 +210,7 @@ aov.formula.list<-function(data,formula,meta.data=factor,post.hoc=TRUE,repeated=
157210 p.values <- data.frame (t(summary(model )[[1 ]][1 : length(names ),5 ,drop = FALSE ]))
158211 dimnames(p.values )<- list (colnames(data )[i ],names )
159212 }
213+
160214 # pairwise t-tests (repeated) or TukeyHSD
161215 if (post.hoc ){
162216 if (! is.null(repeated )){
@@ -184,11 +238,16 @@ aov.formula.list<-function(data,formula,meta.data=factor,post.hoc=TRUE,repeated=
184238 post.h <- NULL
185239 }
186240
241+
187242 results $ p.value [[i ]]<- p.values
188243 results $ post.hoc [[i ]]<- post.h
189244
190245 }
191- return (list (p.values = do.call(" rbind" ,results $ p.value ),post.hoc = do.call(" rbind" ,results $ post.hoc )))
246+ tmp <- do.call(" rbind" ,results $ p.value )
247+ FDR.p <- sapply(1 : ncol(tmp ), function (i ) p.adjust(tmp [,i ],method = FDR ))
248+ colnames(FDR.p )<- colnames(tmp )
249+ p.vals <- data.frame (p.values = tmp ,FDR.p.values = FDR.p )
250+ return (list (p.values = p.vals ,post.hoc = do.call(" rbind" ,results $ post.hoc )))
192251
193252}
194253
@@ -258,35 +317,40 @@ stats.summary <- function(data,comp.obj,formula,sigfigs=3,log=FALSE,rel=1,do.sta
258317 }
259318
260319# function to carry out covariate adjustments
261- # -------------------------
262320covar.adjustment <- function (data ,formula ){
263- # set up that formula objects need to exists in the global environment --- fix this
264- # data--> subjects as rows, measurements as columns
265- # formula <- ~ character vector
266- # lm will be iteratively fit on each variable
267- # model residuals + preadjusted column median will be returned
268- data <- as.data.frame(data )
269- names(data )<- colnames(data )
270- output <- list ()
271- n <- ncol(data )
272- output <- lapply(1 : n ,function (i )
273- {
274- tryCatch(tmp <- as.formula(c(paste(paste(" data$'" ,colnames(data )[i ]," '~" ,sep = " " ),paste(formula ,sep = " +" ),sep = " " ))),
275- error = function (e ){tmp <- as.formula(c(paste(paste(" data[,i]" ," ~" ,sep = " " ),paste(formula ,sep = " +" ),sep = " " )))})
276- fit <- lm(tmp ,data = data )$ residuals
277- matrix (fit ,,1 )
278- })
279- out <- data.frame (do.call(" cbind" ,output ))
280- tryCatch(dimnames(out )<- dimnames(data ), error = function (e ){NULL }) # no clue why this throws errors randomly
281- # add back pre-adjustment column min to all
282- min <- apply(out ,2 ,min , na.rm = T )
283- adj.out <- do.call(" cbind" ,sapply(1 : ncol(out ),function (i )
284- {
285- out [,i ,drop = F ] + abs(min [i ])
286- }))
287- return (adj.out )
288- }
289-
321+ # set up that formula objects need to exists in the global environment --- fix this
322+ # data--> subjects as rows, measurements as columns
323+ # formula <- ~ character vector
324+ # lm will be iteratively fit on each variable
325+ # model residuals + preadjusted column median will be returned
326+
327+ # convert all to numeric
328+ tmp <- dimnames(data )
329+ data <- data.frame (do.call(" cbind" ,lapply(data ,as.numeric )))
330+ dimnames(data )<- tmp
331+ output <- list ()
332+ n <- ncol(data )
333+ output <- lapply(1 : n ,function (i )
334+ {
335+ tryCatch(tmp <- as.formula(c(paste(paste(" data$'" ,colnames(data )[i ]," '~" ,sep = " " ),paste(formula ,sep = " +" ),sep = " " ))),
336+ error = function (e ){tmp <- as.formula(c(paste(paste(" data[,i]" ," ~" ,sep = " " ),paste(formula ,sep = " +" ),sep = " " )))})
337+ fit <- lm(tmp ,data = data )$ residuals
338+ # need to account for missing values
339+ tmp <- data [,i ]
340+ tmp [! is.na(tmp )]<- fit
341+ matrix (tmp ,,1 )
342+ })
343+ out <- data.frame (do.call(" cbind" ,output ))
344+ tryCatch(dimnames(out )<- dimnames(data ), error = function (e ){NULL }) # no clue why this throws errors randomly
345+ # add back pre-adjustment column min to all
346+ min <- apply(out ,2 ,min , na.rm = T )
347+ adj.out <- do.call(" cbind" ,sapply(1 : ncol(out ),function (i )
348+ {
349+ out [,i ,drop = F ] + abs(min [i ])
350+ }))
351+ return (adj.out )
352+ }
353+
290354# helper function for getting statistics for making box plots
291355summarySE <- function (data = NULL , measurevar , groupvars = NULL , na.rm = FALSE ,conf.interval = .95 , .drop = TRUE ) {
292356 require(plyr )
@@ -462,14 +526,15 @@ simple.lme<-function(data,factor,subject,FDR="BH", progress=TRUE,interaction=FAL
462526# multi-LME with formula interface (also see simple.lme) # note random term should be +(1|random.term)
463527formula.lme <- function (data ,formula ,FDR = " BH" , progress = TRUE ){
464528 # data = data.frame of values to test and test factors
465- # factor = object to be tested
466- # subject = identifier for repeated measures
529+
467530 check.get.packages(c(" lme4" ," car" ))
468531 # not sure how to ignore test factors in data, cause error in loop
469532
470533 if (progress == TRUE ){ pb <- txtProgressBar(min = 0 , max = ncol(data ), style = 3 )} else {pb <- NULL }
471534 lmer.p.values <- do.call(" rbind" ,lapply(1 : ncol(data ), function (i ){
472535 if (progress == TRUE ){setTxtProgressBar(pb , i )}
536+ # avoid factor error, should avoid anova error below
537+ if (is.factor(data [,i ])| is.character(data [,i ])) data [,i ]<- fixlf(data [,i ])
473538 mod <- tryCatch(lmer(as.formula(paste0(" data[," ,i ," ]~" , formula )), data = data ),error = function (e ){NULL })
474539 if (is.null(mod )){1 } else {
475540 res <- Anova(mod )
@@ -608,6 +673,7 @@ multi.pairwise.mann.whitney<-function(data,factor,progress=TRUE,FDR="BH",qvalue=
608673 res
609674 }))
610675}
676+
611677# Tests
612678test <- function (){
613679
0 commit comments