Mercurial > repos > workflow4metabolomics > mixmodel4repeated_measures

comparison mixmodel_script.R @ 0:a4d89d47646f draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/tools-metabolomics commit 8d2ca678d973501b60479a8dc3f212eecd56eab8
author workflow4metabolomics
date Mon, 16 May 2022 09:25:01 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a4d89d47646f
1 ####### R functions to perform linear mixed model for repeated measures
2 ####### on a multi var dataset using 3 files as used in W4M
3 ##############################################################################################################
4 lmRepeated2FF <- function(ids, ifixfact, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]],
5 pvalCutof = 0.05, dffOption, visu = visu, tit = "", least.confounded = FALSE, outlier.limit = 3) {
6 ### function to perform linear mixed model with 1 Fixed factor + Time + random factor subject
7 ### based on lmerTest package providing functions giving the same results as SAS proc mixed
8 options(scipen = 50, digits = 5)
9
10 if (!is.numeric(ids[[ivd]])) stop("Dependant variable is not numeric")
11 if (!is.factor(ids[[ifixfact]])) stop("fixed factor is not a factor")
12 if (!is.factor(ids[[itime]])) stop("Repeated factor is not a factor")
13 if (!is.factor(ids[[isubject]])) stop("Random factor is not a factor")
14
15 ## factors
16 time <- ids[[itime]]
17 fixfact <- ids[[ifixfact]]
18 subject <- ids[[isubject]]
19 # dependant variables
20 vd <- ids[[ivd]]
21
22 # argument of the function instead of re re-running ndim <- defColRes(ids, ifixfact, itime)
23 # nfp : number of main factors + model infos (REML, varSubject) + normality test
24 nfp <- ndim[1];
25 # ncff number of comparison of the fixed factor
26 nlff <- ndim[2]; ncff <- ndim[3]
27 # nct number of comparison of the time factor
28 nlt <- ndim[4]; nct <- ndim[5]
29 # nci number of comparison of the interaction
30 nli <- ndim[6]; nci <- ndim[7]
31 # number of all lmer results
32 nresT <- ncff + nct + nci
33 ## initialization of the result vector (1 line)
34 ## 4 * because nresf for : pvalues + Etimates + lower CI + Upper CI
35 res <- data.frame(array(rep(NA, (nfp + 4 * nresT))))
36 colnames(res)[1] <- "resultLM"
37
38 ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip
39 ### after excluding NA, table function is used to seek subjects with only 1 data
40 ids <- ids[!is.na(ids[[ivd]]), ]
41 skip <- length(which(table(ids[[isubject]]) == 1))
42
43 if (skip == 0) {
44
45 mfl <- lmer(vd ~ time + fixfact + time:fixfact + (1 | subject), ids) # lmer remix
46
47 rsum <- summary(mfl, ddf = dffOption)
48 ## test Shapiro Wilks on the residus of the model
49 rShapiro <- shapiro.test(rsum$residuals)
50 raov <- anova(mfl, ddf = dffOption)
51 dlsm1 <- data.frame(difflsmeans(mfl, test.effs = NULL))
52 ddlsm1 <- dlsm1
53 ## save rownames and factor names
54 rn <- rownames(ddlsm1)
55 fn <- ddlsm1[, c(1, 2)]
56 ## writing the results on a single line
57 namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")
58 namesFactPval <- paste("pvalue ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")
59 namesInter <- rownames(ddlsm1)[-c(1:(nct + ncff))]
60 namesEstimate <- paste("estimate ", namesInter)
61 namespvalues <- paste("pvalue ", namesInter)
62 namesFactprinc <- c("pval_time", "pval_trt", "pval_inter")
63 namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")
64
65 namesFactLowerCI <- paste("lowerCI ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")
66 namesLowerCI <- paste("lowerCI ", namesInter, sep = "")
67
68 namesFactUpperCI <- paste("UpperCI ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")
69 namesUpperCI <- paste("UpperCI ", namesInter, sep = "")
70
71
72 ### lmer results on 1 vector row
73 # pvalue of shapiro Wilks test of the residuals
74 res[1, ] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals"
75 res[2, ] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"
76 res[3, ] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"
77 ### 3 principal factors pvalues results + shapiro test => nfp <- 4
78 res[c((nfp - 2):nfp), ] <- raov[, 6]
79 rownames(res)[c((nfp - 2):nfp)] <- namesFactprinc
80
81 #################### Residuals diagnostics for significants variables #########################
82 ### Il at least 1 factor is significant and visu=TRUE NL graphics add to pdf
83 ## ajout JF du passage de la valeur de p-value cutoff
84 if (length(which(raov[, 6] <= pvalCutof)) > 0 & visu == "yes") {
85 diagmflF(mfl, title = tit, pvalCutof = pvalCutof, least.confounded = least.confounded,
86 outlier.limit = outlier.limit)
87 }
88
89 # pvalue of fixed factor comparisons
90 nresf <- nresT
91 res[(nfp + 1):(nfp + nct), ] <- ddlsm1[c(1:nct), 9]
92 res[(nfp + nct + 1):(nfp + nct + ncff), ] <- ddlsm1[(nct + 1):(nct + ncff), 9]
93 rownames(res)[(nfp + 1):(nfp + nct + ncff)] <- namesFactPval
94 res[(nfp + nct + ncff + 1):(nfp + nresf), ] <- ddlsm1[(nct + ncff + 1):(nresT), 9]
95 rownames(res)[(nfp + nct + ncff + 1):(nfp + nresT)] <- namespvalues
96 # Estimate of the difference between levels of factors
97 res[(nfp + nresf + 1):(nfp + nresf + nct), ] <- ddlsm1[c(1:nct), 3]
98 res[(nfp + nresf + nct + 1):(nfp + nresf + nct + ncff), ] <- ddlsm1[(nct + 1):(nct + ncff), 3]
99 rownames(res)[(nfp + nresf + 1):(nfp + nresf + nct + ncff)] <- namesFactEstim
100 res[(nfp + nresf + nct + ncff + 1):(nfp + 2 * nresf), ] <- ddlsm1[(nct + ncff + 1):(nresT), 3]
101 rownames(res)[(nfp + nresf + nct + ncff + 1):(nfp + 2 * nresf)] <- namesEstimate
102 # lower CI of the difference between levels of factors
103 nresf <- nresf + nresT
104 res[(nfp + nresf + 1):(nfp + nresf + nct), ] <- ddlsm1[c(1:nct), 7]
105 res[(nfp + nresf + nct + 1):(nfp + nresf + nct + ncff), ] <- ddlsm1[(nct + 1):(nct + ncff), 7]
106 rownames(res)[(nfp + nresf + 1):(nfp + nresf + nct + ncff)] <- namesFactLowerCI
107 res[(nfp + nresf + nct + ncff + 1):(nfp + 2 * nresf), ] <- ddlsm1[(nct + ncff + 1):(nresf), 7]
108 rownames(res)[(nfp + nresf + nct + ncff + 1):(nfp + nresf + (nresf / 2))] <- namesLowerCI
109 # Upper CI of the difference between levels of factors
110 nresf <- nresf + nresT
111 res[(nfp + nresf + 1):(nfp + nresf + nct), ] <- ddlsm1[c(1:nct), 8]
112 res[(nfp + nresf + nct + 1):(nfp + nresf + nct + ncff), ] <- ddlsm1[(nct + 1):(nct + ncff), 8]
113 rownames(res)[(nfp + nresf + 1):(nfp + nresf + nct + ncff)] <- namesFactUpperCI
114 res[(nfp + nresf + nct + ncff + 1):(nfp + nresf + (nresT)), ] <- ddlsm1[(nct + ncff + 1):(nresT), 8]
115 rownames(res)[(nfp + nresf + nct + ncff + 1):(nfp + nresf + (nresT))] <- namesUpperCI
116
117 } else {
118 ## one of the subject has only one time, subject can't be a random variable
119 ## A repeated measure could be run instead function lme of package nlme, in next version?
120 res[1, ] <- NA
121 cat("\n Computing impossible for feature ", tit, ": at least one subject has only one time.\n")
122 }
123 tres <- data.frame(t(res))
124 rownames(tres)[1] <- nameVar
125 cres <- list(tres, rn, fn)
126 return(cres)
127 }
128
129 ##############################################################################################################
130 lmRepeated1FF <- function(ids, ifixfact = 0, itime, isubject, ivd, ndim, nameVar = colnames(ids)[[ivd]],
131 dffOption, pvalCutof = 0.05) {
132 ### function to perform linear mixed model with factor Time + random factor subject
133 ### based on lmerTest package providing functions giving the same results as SAS proc mixed
134
135
136 if (!is.numeric(ids[[ivd]])) stop("Dependant variable is not numeric")
137 if (!is.factor(ids[[itime]])) stop("Repeated factor is not a factor")
138 if (!is.factor(ids[[isubject]])) stop("Random factor is not a factor")
139 # Could be interesting here to add an experience plan check to give back a specific error message
140 # in case time points are missing for some individuals
141
142 time <- ids[[itime]]
143 subject <- ids[[isubject]]
144 vd <- ids[[ivd]] ## dependant variables (quatitative)
145
146 # nfp : nombre de facteurs principaux + model infos + normality test
147 nfp <- ndim[1]
148 # nlt number of time levels; nct number of comparisons of the time factor
149 nlt <- ndim[4]
150 nct <- ndim[5]
151 # number of all lmer results
152 nresT <- nct
153 ## initialization of the result vector (1 line)
154 res <- data.frame(array(rep(NA, (nfp + 4 * nresT))))
155 colnames(res)[1] <- "resultLM"
156
157 ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip
158 ### after excluding NA, table function is used to seek subjects with only 1 data
159 ids <- ids[!is.na(ids[[ivd]]), ]
160 skip <- length(which(table(ids[[isubject]]) == 1))
161
162 if (skip == 0) {
163
164 mfl <- lmer(vd ~ time + (1 | subject), ids) # lmer remix
165 rsum <- summary(mfl, ddf = dffOption)
166 ## test Shapiro Wilks on the residus of the model
167 rShapiro <- shapiro.test(rsum$residuals)
168 raov <- anova(mfl, ddf = dffOption)
169 ## Sum of square : aov$'Sum Sq', Mean square : aov$`Mean Sq`, proba : aov$`Pr(>F)`
170
171 ## Test of all differences estimates between levels as SAS proc mixed.
172 ## results are in diffs.lsmeans.table dataframe
173 ## test.effs=NULL perform all pairs comparisons including interaction effect
174 dlsm1 <- data.frame(difflsmeans(mfl, test.effs = NULL))
175 ddlsm1 <- dlsm1
176 ## save rownames and factor names
177 rn <- rownames(ddlsm1)
178 fn <- ddlsm1[, c(1, 2)]
179 ## writing the results on a single line
180 namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct))], sep = "")
181 namesFactPval <- paste("pvalue ", rownames(ddlsm1)[c(1:(nct))], sep = "")
182 namesFactprinc <- "pval_time"
183 namesLowerCI <- paste("lowerCI ", rownames(ddlsm1)[c(1:(nct))], sep = "")
184 namesUpperCI <- paste("upperCI ", rownames(ddlsm1)[c(1:(nct))], sep = "")
185
186 ### lmer results on 1 vector
187 # pvalue of shapiro Wilks test of the residuals
188 res[1, ] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals"
189 res[2, ] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"
190 res[3, ] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"
191
192 ### factor time pvalue results + shapiro test
193 res[nfp, ] <- raov[, 6]
194 rownames(res)[nfp] <- namesFactprinc
195
196 # pvalues of time factor levels comparisons
197 res[(nfp + 1):(nfp + nct), ] <- ddlsm1[c(1:nct), 9]
198 rownames(res)[(nfp + 1):(nfp + nct)] <- namesFactPval
199
200 # Estimates of time factor levels
201 nresf <- nresT
202 res[(nfp + nresf + 1):(nfp + nresf + nct), ] <- ddlsm1[c(1:nct), 3]
203 rownames(res)[(nfp + nresf + 1):(nfp + nresf + nct)] <- namesFactEstim
204
205 # Lower CI of the difference between levels of factors
206 # nresf is incremeted
207 nresf <- nresf + nresT
208 res[(nfp + nresf + 1):(nfp + nresf + nct), ] <- ddlsm1[c(1:nct), 7]
209 rownames(res)[(nfp + nresf + 1):(nfp + nresf + nct)] <- namesLowerCI
210 # Lower CI of the difference between levels of factors
211 # nresf is incremeted
212 nresf <- nresf + nresT
213 res[(nfp + nresf + 1):(nfp + nresf + nct), ] <- ddlsm1[c(1:nct), 8]
214 rownames(res)[(nfp + nresf + 1):(nfp + nresf + nct)] <- namesUpperCI
215
216
217 } else {
218 ## one of the subject has only one time, subject can't be a random variable
219 ## A repeated measure could be run instead function lme of package nlme, next version
220 res[1, ] <- NA
221 cat("\n Computing impossible for feature ", colnames(ids)[4], ": at least one subject has only one time.\n")
222 }
223 tres <- data.frame(t(res))
224 rownames(tres)[1] <- nameVar
225 cres <- list(tres, rn, fn)
226 return(cres)
227
228 }
229
230 ##############################################################################################################
231 defColRes <- function(ids, ifixfact, itime) {
232 ## define the size of the result file depending on the numbers of levels of the fixed and time factor.
233 ## Numbers of levels define the numbers of comparisons with pvalue and estimate of the difference.
234 ## The result file also contains the pvalue of the fixed factor, time factor and interaction
235 ## plus Shapiro normality test. This is define by nfp
236 ## subscript of fixed factor=0 means no other fixed factor than "time"
237 if (ifixfact > 0) {
238 nfp <- 6 # shapiro + subject variance + REML + time + fixfact + interaction
239 time <- ids[[itime]]
240 fixfact <- ids[[ifixfact]]
241
242 cat("\n Levels fixfact: ", levels(fixfact))
243 cat("\n Levels time: ", levels(time))
244
245 # ncff number of comparisons of the fixed factor (nlff number of levels of fixed factor)
246 nlff <- length(levels(fixfact))
247 ncff <- (nlff * (nlff - 1)) / 2
248 # nct number of comparison of the time factor (nlt number of levels of time factor)
249 nlt <- length(levels(time))
250 nct <- (nlt * (nlt - 1)) / 2
251 # nci number of comparison of the interaction
252 nli <- nlff * nlt
253 nci <- (nli * (nli - 1)) / 2
254 ndim <- c(NA, NA, NA, NA, NA, NA, NA)
255
256 ndim[1] <- nfp # pvalues of fixed factor, time factor and interaction (3columns) and shapiro test pvalue
257 ndim[2] <- nlff # number of levels of fixed factor
258 ndim[3] <- ncff # number of comparisons (2by2) of the fixed factor
259 ndim[4] <- nlt # number of levels of time factor
260 ndim[5] <- nct # number of comparisons (2by2) of the time factor
261 ndim[6] <- nli # number of levels of interaction
262 ndim[7] <- nci # number of comparisons (2by2) of the interaction
263
264 }
265 else {
266 nfp <- 4 # Mandatory columns: shapiro + subject variance + REML + time
267 time <- ids[[itime]]
268 # nct number of comparison of the time factor
269 nlt <- length(levels(time))
270 nct <- (nlt * (nlt - 1)) / 2
271 ndim <- c(NA, NA, NA, NA, NA, NA, NA)
272
273 ndim[1] <- nfp # pvalues of shapiro test + subject variance + REML + time factor
274 ndim[4] <- nlt # number of levels of time factor
275 ndim[5] <- nct # number of comparisons (2by2) of the time factor
276 }
277 return(ndim)
278 }
279
280 ##############################################################################################################
281 lmixedm <- function(datMN,
282 samDF,
283 varDF,
284 fixfact, time, subject,
285 logtr = "none",
286 pvalCutof = 0.05,
287 pvalcorMeth = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")[7],
288 dffOption,
289 visu = "no",
290 least.confounded = FALSE,
291 outlier.limit = 3,
292 pdfC,
293 pdfE
294 ) {
295 sampids <- samDF
296 dataMatrix <- datMN
297 varids <- varDF
298
299 options("scipen" = 50, "digits" = 5)
300 pvalCutof <- as.numeric(pvalCutof)
301
302 cat("\n dff computation method=", dffOption)
303 ### Function running lmer function on a set of variables described in
304 ### 3 different dataframes as used by W4M
305 ### results are merge with the metadata variables varids
306 ### ifixfact, itime, isubject are subscripts of the dependant variables. if only time factor the ifixfat is set to 0 and no diag is performed (visu="no")
307 if (fixfact == "none") {
308 ifixfact <- 0 ; visu <- "no"
309 } else ifixfact <- which(colnames(sampids) == fixfact)
310 itime <- which(colnames(sampids) == time)
311 isubject <- which(colnames(sampids) == subject)
312
313 lmmds <- dataMatrix
314 if (logtr != "log10" & logtr != "log2") logtr <- "none"
315 if (logtr == "log10") lmmds <- log10(lmmds + 1)
316 if (logtr == "log2") lmmds <- log2(lmmds + 1)
317
318 dslm <- cbind(sampids, lmmds)
319
320 nvar <- ncol(lmmds)
321 firstvar <- ncol(sampids) + 1
322 lastvar <- firstvar + ncol(lmmds) - 1
323
324 if (ifixfact > 0) dslm[[ifixfact]] <- factor(dslm[[ifixfact]])
325 dslm[[itime]] <- factor(dslm[[itime]])
326 dslm[[isubject]] <- factor(dslm[[isubject]])
327 ## call defColres to define the numbers of test and so the number of columns of results
328 ## depends on whether or not there is a fixed factor with time. If only time factor ifixfact=0
329 if (ifixfact > 0) {
330 ndim <- defColRes(dslm[, c(ifixfact, itime)], ifixfact = 1, itime = 2)
331 nColRes <- ndim[1] + (4 * (ndim[3] + ndim[5] + ndim[7]))
332 firstpval <- ndim[1] - 2
333 lastpval <- ndim[1] + ndim[3] + ndim[5] + ndim[7]
334 } else {
335 ndim <- defColRes(dslm[, itime], ifixfact = 0, itime = 1)
336 nColRes <- ndim[1] + (4 * (ndim[5]))
337 firstpval <- ndim[1]
338 lastpval <- ndim[1] + ndim[5]
339 }
340 ## initialisation of the result file
341 resLM <- data.frame(array(rep(NA, nvar * nColRes), dim = c(nvar, nColRes)))
342 rownames(resLM) <- rownames(varids)
343
344 ## PDF initialisation
345 if (visu == "yes") {
346 pdf(pdfC, onefile = TRUE, height = 15, width = 30)
347 par(mfrow = c(1, 3))
348 }
349
350
351 for (i in firstvar:lastvar) {
352
353 subds <- dslm[, c(ifixfact, itime, isubject, i)]
354
355 tryCatch({
356 if (ifixfact > 0)
357 reslmer <- lmRepeated2FF(subds, ifixfact = 1, itime = 2, isubject = 3, ivd = 4, ndim = ndim, visu = visu,
358 tit = colnames(dslm)[i], pvalCutof = pvalCutof,
359 dffOption = dffOption, least.confounded = least.confounded,
360 outlier.limit = outlier.limit)
361 else
362 reslmer <- lmRepeated1FF(subds, ifixfact = 0, itime = 1, isubject = 2, ivd = 3, ndim = ndim,
363 pvalCutof = pvalCutof, dffOption = dffOption)
364
365 resLM[i - firstvar + 1, ] <- reslmer[[1]]
366 }, error = function(e) {
367 cat("\nERROR with ", rownames(resLM)[i - firstvar + 1], ": ", conditionMessage(e), "\n")
368 }
369 )
370
371 }
372
373 if (exists("reslmer")) {
374 colnames(resLM) <- colnames(reslmer[[1]])
375 labelRow <- reslmer[[2]]
376 factorRow <- reslmer[[3]]
377 } else {
378 stop("\n- - - - -\nModel computation impossible for every single variables in the dataset: no result returned.\n- - - - -\n")
379 }
380
381 if (visu == "yes") dev.off()
382
383
384 ## pvalue correction with p.adjust library multtest
385 ## Possible methods of pvalue correction
386 AdjustMeth <- c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
387 if (length(which(pvalcorMeth == AdjustMeth)) == 0) pvalcorMeth <- "none"
388
389 if (pvalcorMeth != "none") {
390 for (k in firstpval:lastpval) {
391 resLM[[k]] <- p.adjust(resLM[[k]], method = pvalcorMeth, n = dim(resLM[k])[[1]])
392
393 }
394 }
395
396 ## for each variables, set pvalues to NA and estimates = 0 when pvalue of factor > pvalCutof value define by user
397 if (ifixfact > 0) {
398 ## time effect
399 resLM[which(resLM[, firstpval] > pvalCutof), c((lastpval + 1):(lastpval + ndim[5]))] <- 0
400 resLM[which(resLM[, firstpval] > pvalCutof), c((ndim[1] + 1):(ndim[1] + ndim[5]))] <- NA
401 ## treatment effect
402 resLM[which(resLM[, firstpval + 1] > pvalCutof), c((lastpval + ndim[5] + 1):(lastpval + ndim[5] + ndim[3]))] <- 0
403 resLM[which(resLM[, firstpval + 1] > pvalCutof), c((ndim[1] + ndim[5] + 1):(ndim[1] + ndim[5] + ndim[3]))] <- NA
404 ## interaction effect
405 resLM[which(resLM[, firstpval + 2] > pvalCutof), c((lastpval + ndim[5] + ndim[3] + 1):(lastpval + ndim[5] + ndim[3] + ndim[7]))] <- 0
406 resLM[which(resLM[, firstpval + 2] > pvalCutof), c((ndim[1] + ndim[5] + ndim[3] + 1):(ndim[1] + ndim[5] + ndim[3] + ndim[7]))] <- NA
407 } else {
408 ## time effect only
409 resLM[which(resLM[, firstpval] > pvalCutof), c((lastpval + 1):(lastpval + ndim[5]))] <- 0
410 resLM[which(resLM[, firstpval] > pvalCutof), c((firstpval + 1):(firstpval + ndim[5]))] <- NA
411 }
412
413 ## for each variable, estimates plots are performed if at least one factor is significant after p-value correction
414 pdf(pdfE, onefile = TRUE, height = 15, width = 30)
415
416 ## for each variable (in row)
417 for (i in seq_len(nrow(resLM))) {
418
419 ## if any fixed factor + time factor
420 if (ifixfact > 0)
421
422 ## if any main factor after p-value correction is significant -> plot estimates and time course
423 if (length(which(resLM[i, c(4:6)] < pvalCutof)) > 0) {
424
425 ## Plot of time course by fixfact : data prep with factors and quantitative var to be plot
426 subv <- dslm[, colnames(dslm) == rownames(resLM)[i]]
427 subds <- data.frame(dslm[[ifixfact]], dslm[[itime]], dslm[[isubject]], subv)
428 libvar <- c(fixfact, time, subject)
429 colnames(subds) <- c(libvar, rownames(resLM)[i])
430
431 ## Plot of estimates with error bars for all fixed factors and interaction
432 rddlsm1 <- t(resLM[i, ])
433 pval <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "pvalue"]
434 esti <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "estima"]
435 loci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "lowerC"]
436 upci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "UpperC"]
437 rddlsm1 <- data.frame(pval, esti, loci, upci, factorRow)
438 colnames(rddlsm1) <- c("p.value", "Estimate", "Lower.CI", "Upper.CI", colnames(factorRow))
439 rownames(rddlsm1) <- labelRow
440
441 ## function for plotting these 2 graphs
442 plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)
443
444 }
445
446 ## if only a time factor
447 if (ifixfact == 0)
448
449 ## if time factor after p-value correction is significant -> plot time course
450 if (length(which(resLM[i, 4] < pvalCutof)) > 0) {
451
452 ## Plot of time course : data prep with factors and quantitative var to be plot
453 subv <- dslm[, colnames(dslm) == rownames(resLM)[i]]
454 subds <- data.frame(dslm[[itime]], dslm[[isubject]], subv)
455 libvar <- c(time, subject)
456 colnames(subds) <- c(libvar, rownames(resLM)[i])
457
458 ## Plot of estimates with error bars for all fixed factors and interaction
459 rddlsm1 <- t(resLM[i, ])
460 pval <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "pvalue"]
461 esti <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "estima"]
462 loci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "lowerC"]
463 upci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "upperC"]
464 rddlsm1 <- data.frame(pval, esti, loci, upci, factorRow)
465 colnames(rddlsm1) <- c("p.value", "Estimate", "Lower.CI", "Upper.CI", colnames(factorRow))
466 rownames(rddlsm1) <- labelRow
467
468 ## function for plotting these 2 graphs
469 plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)
470
471 }
472
473 }
474 dev.off()
475
476 ## return result file with pvalues and estimates (exclude confidence interval used for plotting)
477 iCI <- which(substr(colnames(resLM), 4, 7) == "erCI")
478 resLM <- resLM[, - iCI]
479 resLM <- cbind(varids, resLM)
480 return(resLM)
481 }