Mercurial > repos > jfrancoismartin > mixmodel4repeated_measures

comparison mixmodel_script.R @ 0:1422de181204 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/mixmodel4repeated_measures commit 6ea32b3182383c19e5333201d2385a61d8da3d50
author jfrancoismartin
date Wed, 10 Oct 2018 05:18:42 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1422de181204
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 , tit = "", least.confounded = FALSE, outlier.limit =3)
6 {
7 ### function to perform linear mixed model with 1 Fixed factor + Time + random factor subject
8 ### based on lmerTest package providing functions giving the same results as SAS proc mixed
9 options(scipen = 50, digits = 5)
10
11 if (!is.numeric(ids[[ivd]])) {stop("Dependant variable is not numeric")}
12 if (!is.factor(ids[[ifixfact]])) {stop("fixed factor is not a factor")}
13 if (!is.factor(ids[[itime]])) {stop("Repeated factor is not a factor")}
14 if (!is.factor(ids[[isubject]])) {stop("Random factor is not a factor")}
15 # a ce stade, il faudrait pr?voir des tests sur la validit? du plan d'exp?rience
16
17 time <- ids[[itime]]
18 fixfact <- ids[[ifixfact]]
19 subject <- ids[[isubject]]
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 # ## NL add
48 # ### DEPLACE APRES CALCUL PVALUES AJUSTEES ET NE FAIRE QUE SI AU MOINS 1 FACTEUR SIGNIFICATIF
49 # if(visu) diagmflF(mfl, title = tit, least.confounded = least.confounded, outlier.limit = outlier.limit)
50 # ## end of NL add
51
52 rsum <- summary(mfl,ddf = dffOption)
53 ## test Shapiro Wilks on the residus of the model
54 rShapiro <- shapiro.test(rsum$residuals)
55 raov <- anova(mfl,ddf = dffOption)
56 dlsm1 <- data.frame(difflsmeans(mfl,test.effs=NULL))
57 ddlsm1 <- dlsm1
58 ## save rownames and factor names
59 rn <- rownames(ddlsm1)
60 fn <- ddlsm1[,c(1,2)]
61 ## writing the results on a single line
62 namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
63 namesFactPval <- paste("pvalue ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
64 namesInter <- rownames(ddlsm1)[-c(1:(nct+ncff))]
65 #ncI <- nchar(namesInter)
66 namesEstimate <- paste("estimate ",namesInter)
67 namespvalues <- paste("pvalue ",namesInter)
68 namesFactprinc <- c("pval_time","pval_trt","pval_inter")
69 namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
70
71 namesFactLowerCI <- paste("lowerCI ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
72 namesLowerCI <- paste("lowerCI ",namesInter,sep="")
73
74 namesFactUpperCI <- paste("UpperCI ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
75 namesUpperCI <- paste("UpperCI ",namesInter,sep="")
76
77
78 ### lmer results on 1 vector row
79 # pvalue of shapiro Wilks test of the residuals
80 res[1,] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals"
81 res[2,] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"
82 res[3,] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"
83 ### 3 principal factors pvalues results + shapiro test => nfp <- 4
84 res[c((nfp-2):nfp),] <- raov[,6]; rownames(res)[c((nfp-2):nfp)] <- namesFactprinc
85
86 #################### Residuals diagnostics for significants variables #########################
87 ### Il at least 1 factor is significant and visu=TRUE NL graphics add to pdf
88 ## ajout JF du passage de la valeur de p-value cutoff
89 if (length(which(raov[,6]<=pvalCutof))>0 & visu == 'yes') {
90 diagmflF(mfl, title = tit, pvalCutof = pvalCutof, least.confounded = least.confounded,
91 outlier.limit = outlier.limit)
92
93 cat(" Signif ",pvalCutof)
94
95
96 }
97
98 # pvalue of fixed factor comparisons
99 nresf <- nresT
100 res[(nfp+1):(nfp+nct),] <- ddlsm1[c(1:nct),9]
101 res[(nfp+nct+1):(nfp+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),9]
102 rownames(res)[(nfp+1):(nfp+nct+ncff)] <- namesFactPval
103 res[(nfp+nct+ncff+1):(nfp+nresf),] <- ddlsm1[(nct+ncff+1):(nresT),9]
104 rownames(res)[(nfp+nct+ncff+1):(nfp+nresT)] <- namespvalues
105 # Estimate of the difference between levels of factors
106 res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),3]
107 res[(nfp+nresf+nct+1):(nfp+nresf+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),3]
108 rownames(res)[(nfp+nresf+1):(nfp+nresf+nct+ncff)] <- namesFactEstim
109 res[(nfp+nresf+nct+ncff+1):(nfp+2*nresf),] <- ddlsm1[(nct+ncff+1):(nresT),3]
110 rownames(res)[(nfp+nresf+nct+ncff+1):(nfp+2*nresf)] <- namesEstimate
111 # lower CI of the difference between levels of factors
112 nresf <- nresf + nresT
113 res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),7]
114 res[(nfp+nresf+nct+1):(nfp+nresf+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),7]
115 rownames(res)[(nfp+nresf+1):(nfp+nresf+nct+ncff)] <- namesFactLowerCI
116 res[(nfp+nresf+nct+ncff+1):(nfp+2*nresf),] <- ddlsm1[(nct+ncff+1):(nresf),7]
117 rownames(res)[(nfp+nresf+nct+ncff+1):(nfp+nresf+(nresf/2))] <- namesLowerCI
118 # Upper CI of the difference between levels of factors
119 nresf <- nresf + nresT
120 res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),8]
121 res[(nfp+nresf+nct+1):(nfp+nresf+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),8]
122 rownames(res)[(nfp+nresf+1):(nfp+nresf+nct+ncff)] <- namesFactUpperCI
123 res[(nfp+nresf+nct+ncff+1):(nfp+nresf+(nresT)),] <- ddlsm1[(nct+ncff+1):(nresT),8]
124 rownames(res)[(nfp+nresf+nct+ncff+1):(nfp+nresf+(nresT))] <- namesUpperCI
125
126
127 }
128 else
129 ## one of the subject has only one time, subject can't be a random variable
130 ## A repeated measure could be run instead function lme of package nlme, next version
131 { res[1,] <- NA
132 #cat("impossible computing\n")
133
134 # # ## NL add (useless)
135 # if(visu){
136 # grid.arrange(ggplot(data.frame()) + geom_point() + xlim(-1, 1) + ylim(-1, 1)+
137 # annotate("text", x = 0, y = 0, label = "impossible computing")+
138 # xlab(NULL) + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank())+
139 # ylab(NULL) + theme(axis.text.y=element_blank(),axis.ticks.y=element_blank())+
140 # theme(panel.grid.minor = element_blank() ,
141 # panel.grid.major = element_blank() ,
142 # panel.background = element_rect(fill = "white"))
143 # , top = textGrob(tit,gp=gpar(fontsize=40,font=4)))
144 #
145 # }
146 # # ## end of NL add
147
148 }
149 tres <- data.frame(t(res)); rownames(tres)[1] <- nameVar
150 cres <- list(tres,rn, fn)
151 return(cres)
152 }
153
154 ##############################################################################################################
155 lmRepeated1FF <- function(ids, ifixfact=0, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]],
156 dffOption,pvalCutof=0.05)
157 {
158 ### function to perform linear mixed model with factor Time + random factor subject
159 ### based on lmerTest package providing functions giving the same results as SAS proc mixed
160
161 if (!is.numeric(ids[[ivd]])) {stop("Dependant variable is not numeric")}
162 if (!is.factor(ids[[itime]])) {stop("Repeated factor is not a factor")}
163 if (!is.factor(ids[[isubject]])) {stop("Random factor is not a factor")}
164 # a ce stade, il faudrait pr?voir des tests sur la validit? du plan d'exp?rience
165
166 time <- ids[[itime]]
167 subject <- ids[[isubject]]
168 vd <- ids[[ivd]] ## dependant variables (quatitative)
169
170 # ndim <- defColRes(ids,0,itime)
171 # nfp : nombre de facteurs principaux + model infos + normality test
172 nfp <- ndim[1]
173 # nct number of comparison of the time factor
174 nlt <- ndim[4] ; nct <- ndim[5]
175 # number of all lmer results
176 nresf <- nct
177 ## initialization of the result vector (1 line)
178 res <- data.frame(array(rep(NA,(nfp+2*nresf))))
179 colnames(res)[1] <- "resultLM"
180
181 ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip
182 ### after excluding NA, table function is used to seek subjects with only 1 data
183 ids <- ids[!is.na(ids[[ivd]]),]
184 skip <- length(which(table(ids[[isubject]])==1))
185
186 if (skip==0) {
187
188 mfl <- lmer( vd ~ time + (1| subject), ids) # lmer remix
189 rsum <- summary(mfl,ddf = dffOption)
190 ## test Shapiro Wilks on the residus of the model
191 rShapiro <- shapiro.test(rsum$residuals)
192 raov <- anova(mfl,ddf = dffOption)
193 ## Sum of square : aov$'Sum Sq', Mean square : aov$`Mean Sq`, proba : aov$`Pr(>F)`
194
195 ## Test of all differences estimates between levels as SAS proc mixed.
196 ## results are in diffs.lsmeans.table dataframe
197 ## test.effs=NULL perform all pairs comparisons including interaction effect
198 dlsm1 <- difflsmeans(mfl,test.effs=NULL)
199 ddlsm1 <- dlsm1$diffs.lsmeans.table
200
201 ## writing the results on a single line
202 namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct))],sep="")
203 namesFactPval <- paste("pvalue ",rownames(ddlsm1)[c(1:(nct))],sep="")
204 namesFactprinc <- "pval_time"
205
206 ### lmer results on 1 vector
207 # pvalue of shapiro Wilks test of the residuals
208 res[1,] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals"
209 res[2,] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"
210 res[3,] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"
211
212 ### principal factor time pvalue results + shapiro test
213 res[nfp,] <- raov[,6]; rownames(res)[nfp] <- namesFactprinc
214 # pvalue of fixed factor comparisons
215 res[(nfp+1):(nfp+nct),] <- ddlsm1[c(1:nct),7]
216 rownames(res)[(nfp+1):(nfp+nct)] <- namesFactPval
217
218 # Estimate of the difference between levels of factors
219 res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),1]
220 rownames(res)[(nfp+nresf+1):(nfp+nresf+nct)] <- namesFactEstim
221 }
222 else
223 ## one of the subject has only one time, subject can't be a random variable
224 ## A repeated measure could be run instead function lme of package nlme, next version
225 { res[1,] <- NA
226 #cat("traitement impossible\n")
227 }
228 tres <- data.frame(t(res)); rownames(tres)[1] <- nameVar
229 return(tres)
230 }
231
232 ##############################################################################################################
233 defColRes <- function(ids, ifixfact, itime) {
234 ## define the size of the result file depending on the numbers of levels of the fixed and time factor.
235 ## Numbers of levels define the numbers of comparisons with pvalue and estimate of the difference.
236 ## The result file also contains the pvalue of the fixed factor, time factor and interaction
237 ## plus Shapiro normality test. This is define by nfp
238 ## subscript of fixed factor=0 means no other fixed factor than "time"
239 if (ifixfact>0){
240 nfp <- 6 # shapiro+time+fixfact+interaction+ others....
241 time <- ids[[itime]]
242 fixfact <- ids[[ifixfact]]
243
244 cat("\n levels fixfact",levels(fixfact))
245 cat("\n levels time",levels(time))
246
247 # ncff number of comparisons of the fixed factor (nlff number of levels of fixed factor)
248 nlff <- length(levels(fixfact)); ncff <- (nlff*(nlff-1))/2
249 # nct number of comparison of the time factor (nlt number of levels of time factor)
250 nlt <- length(levels(time)); nct <- (nlt*(nlt-1))/2
251 # nci number of comparison of the interaction
252 nli <- nlff*nlt; nci <- (nli*(nli-1))/2
253 ndim <- c(NA,NA,NA,NA,NA,NA,NA)
254
255 ndim[1] <- nfp # pvalues of fixed factor, time factor and interaction (3columns) and shapiro test pvalue
256 ndim[2] <- nlff # number of levels of fixed factor
257 ndim[3] <- ncff # number of comparisons (2by2) of the fixed factor
258 ndim[4] <- nlt # number of levels of time factor
259 ndim[5] <- nct # number of comparisons (2by2) of the time factor
260 ndim[6] <- nli # number of levels of interaction
261 ndim[7] <- nci # number of comparisons (2by2) of the interaction
262
263 }
264 else {
265 nfp <- 4 # shapiro+time
266 time <- ids[[itime]]
267 # nct number of comparison of the time factor
268 nlt <- length(levels(time)); nct <- (nlt*(nlt-1))/2
269 ndim <- c(NA,NA,NA,NA,NA,NA,NA)
270
271 ndim[1] <- nfp # pvalues of time factor and shapiro test pvalue
272 ndim[4] <- nlt # number of levels of time factor
273 ndim[5] <- nct # number of comparisons (2by2) of the time factor
274 }
275 return(ndim)
276 }
277
278 ##############################################################################################################
279 lmixedm <- function(datMN,
280 samDF,
281 varDF,
282 fixfact, time, subject,
283 logtr = "none",
284 pvalCutof = 0.05,
285 pvalcorMeth = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")[7],
286 dffOption,
287 visu = "no",
288 least.confounded = FALSE,
289 outlier.limit = 3,
290 pdfC,
291 pdfE
292 )
293 {
294 sampids <- samDF
295 dataMatrix <- datMN
296 varids <- varDF
297
298 options("scipen" = 50, "digits" = 5)
299 pvalCutof <- as.numeric(pvalCutof)
300
301 cat("\n dff computation method=",dffOption)
302 ### Function running lmer function on a set of variables described in
303 ### 3 different dataframes as used by W4M
304 ### results are merge with the metadata variables varids
305 ### ifixfact, itime, isubject are subscripts of the dependant variables
306 if (fixfact=="none") ifixfact <-0 else ifixfact <- which(colnames(sampids)==fixfact)
307 itime <- which(colnames(sampids)==time)
308 isubject <- which(colnames(sampids)==subject)
309
310 #lmmds <- dataMatrix[,-1]
311
312 lmmds <- dataMatrix
313 if (logtr!="log10" & logtr!="log2") logtr <- "none"
314 if (logtr=="log10") lmmds <- log10(lmmds+1)
315 if (logtr== "log2") lmmds <- log2(lmmds+1)
316
317 #idsamp <- dataMatrix[,1]
318 #lmmds <- t(lmmds)
319 dslm <- cbind(sampids,lmmds)
320
321 nvar <- ncol(lmmds); firstvar <- ncol(sampids)+1; lastvar <- firstvar+ncol(lmmds)-1
322
323 dslm[[ifixfact]] <- factor(dslm[[ifixfact]])
324 dslm[[itime]] <- factor(dslm[[itime]])
325 dslm[[isubject]] <- factor(dslm[[isubject]])
326 ## call defColres to define the numbers of test and so the number of columns of results
327 ## depends on whether or not there is a fixed factor with time. If only time factor ifixfact=0
328 if (ifixfact>0) {
329 ndim <- defColRes(dslm[,c(ifixfact,itime)],ifixfact=1,itime=2)
330 nColRes <- ndim[1]+(4*(ndim[3]+ndim[5]+ndim[7]))
331 firstpval <- ndim[1]-2
332 lastpval <- ndim[1]+ndim[3]+ndim[5]+ndim[7]
333 } else
334 {
335 ndim <- defColRes(dslm[,itime],ifixfact=0,itime=1)
336 nColRes <- ndim[1]+(2*(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 ############### test ecriture dans pdf
345 if(visu == "yes") {
346 pdf(pdfC, onefile=TRUE, height = 15, width = 30)
347 par(mfrow=c(1,3))
348 }
349 ############### fin test ecriture dans pdf
350 ## pour test : lastvar <- 15
351 cat("\n pvalCutof ", pvalCutof)
352
353 for (i in firstvar:lastvar) {
354
355 ## NL modif
356 cat("\n[",colnames(dslm)[i],"] ")
357 ## end of NL modif
358
359 subds <- dslm[,c(ifixfact,itime,isubject,i)]
360
361 ## NL modif
362 tryCatch({
363 if (ifixfact>0)
364 reslmer <- lmRepeated2FF(subds,ifixfact=1,itime=2,isubject=3, ivd=4, ndim=ndim, visu = visu,
365 tit = colnames(dslm)[i], pvalCutof=pvalCutof,
366 dffOption=dffOption,least.confounded = least.confounded,
367 outlier.limit = outlier.limit)
368 else
369 reslmer <- lmRepeated1FF(subds,ifixfact=0,1,2, ivd=3, ndim=ndim, pvalCutof=pvalCutof,dffOption)
370 ## end of NL modif
371 resLM[i-firstvar+1,] <- reslmer[[1]]
372 }, error=function(e){cat("ERROR : ",conditionMessage(e), "\n");})
373 if (i==firstvar) {
374 colnames(resLM) <- colnames(reslmer[[1]])
375 labelRow <- reslmer[[2]]
376 factorRow <- reslmer[[3]]
377 }
378 }
379 ## for debug : ifixfact=1;itime=2;isubject=3; ivd=4;tit = colnames(dslm)[i]; ids <- subds
380
381
382 ## NL add
383 if(visu == "yes") dev.off()
384 ## end of NL add
385
386 ## pvalue correction with p.adjust library multtest
387 ## Possible methods of pvalue correction
388 AdjustMeth <- c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr","none")
389 if (length(which(pvalcorMeth == AdjustMeth))==0) pvalcorMeth <- "none"
390
391 if (pvalcorMeth !="none") {
392 for (k in firstpval:lastpval){
393 resLM[[k]]=p.adjust(resLM[[k]], method=pvalcorMeth, n=dim(resLM[k])[[1]])
394
395 }
396 }
397
398 ## for each variables, set pvalues to NA and estimates = 0 when pvalue of factor > pvalCutof value define by user
399 if (ifixfact>0) {
400 ## time effect
401 resLM[which(resLM[,firstpval]> pvalCutof),c((lastpval+1):(lastpval+ndim[5]))] <- 0
402 resLM[which(resLM[,firstpval]> pvalCutof),c((ndim[1]+1):(ndim[1]+ndim[5]))] <- NA
403 ## treatment effect
404 resLM[which(resLM[,firstpval+1]> pvalCutof),c((lastpval+ndim[5]+1):(lastpval+ndim[5]+ndim[3]))] <- 0
405 resLM[which(resLM[,firstpval+1]> pvalCutof),c((ndim[1]+ndim[5]+1):(ndim[1]+ndim[5]+ndim[3]))] <- NA
406 ## interaction effect
407 resLM[which(resLM[,firstpval+2]> pvalCutof),c((lastpval+ndim[5]+ndim[3]+1):(lastpval+ndim[5]+ndim[3]+ndim[7]))] <- 0
408 resLM[which(resLM[,firstpval+2]> pvalCutof),c((ndim[1]+ndim[5]+ndim[3]+1):(ndim[1]+ndim[5]+ndim[3]+ndim[7]))] <- NA
409 } else {
410 ## time effect only
411 resLM[which(resLM[,firstpval]> pvalCutof),c((lastpval+1):(lastpval+ndim[5]))] <- 0
412 resLM[which(resLM[,firstpval]> pvalCutof),c((firstpval+1):(firstpval+ndim[5]))] <- NA
413 }
414
415 ## for each variable, estimates plots are performed if at least one factor is significant after p-value correction
416 pdf(pdfE, onefile=TRUE, height = 15, width = 30)
417 #par(mfrow=c(2,2))
418
419 ## for each variable (in row)
420 for (i in 1:nrow(resLM)) {
421 #cat("\n",rownames(resLM)[i])
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 #colnames(subds) <- c(colnames(dslm)[ifixfact],colnames(dslm)[itime],colnames(dslm)[isubject],rownames(resLM)[i] <- rownames(resLM)[i] )
429 libvar <- c(fixfact,time,subject)
430 colnames(subds) <- c(libvar,rownames(resLM)[i])
431
432 ## Plot of estimates with error bars for all fixed factors and interaction
433 rddlsm1 <- t(resLM[i,])
434 pval <- rddlsm1[substr(rownames(rddlsm1),1,6)=="pvalue"]
435 esti <- rddlsm1[substr(rownames(rddlsm1),1,6)=="estima"]
436 loci <- rddlsm1[substr(rownames(rddlsm1),1,6)=="lowerC"]
437 upci <- rddlsm1[substr(rownames(rddlsm1),1,6)=="UpperC"]
438 rddlsm1 <- data.frame(pval,esti,loci,upci,factorRow)
439 colnames(rddlsm1) <- c("p.value","Estimate","Lower.CI","Upper.CI",colnames(factorRow))
440 rownames(rddlsm1) <- labelRow
441
442 ## function for plotting these 2 graphs
443 plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)
444
445 }
446 }
447 dev.off()
448
449 ## return result file with pvalues and estimates (exclude confidence interval used for plotting)
450 iCI <- which(substr(colnames(resLM),4,7)=="erCI")
451 resLM <- resLM[,-iCI]
452 resLM <- cbind(varids,resLM)
453 return(resLM)
454 }
455
456