--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mixmodel_script.R	Wed Oct 10 05:18:42 2018 -0400
@@ -0,0 +1,456 @@
+#######  R functions to perform linear mixed model for repeated measures 
+#######  on a multi var dataset using 3 files as used in W4M
+##############################################################################################################
+lmRepeated2FF <- function(ids, ifixfact, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]], 
+                          pvalCutof=0.05,dffOption, visu , tit = "", least.confounded = FALSE, outlier.limit =3) 
+   {
+   ### function to perform linear mixed model with 1 Fixed factor + Time + random factor subject
+   ### based on lmerTest package providing functions giving the same results as SAS proc mixed
+   options(scipen = 50, digits = 5)
+   
+   if (!is.numeric(ids[[ivd]]))     {stop("Dependant variable is not numeric")}
+   if (!is.factor(ids[[ifixfact]])) {stop("fixed factor is not a factor")}
+   if (!is.factor(ids[[itime]]))    {stop("Repeated factor is not a factor")}
+   if (!is.factor(ids[[isubject]])) {stop("Random factor is not a factor")}
+   # a ce stade, il faudrait pr?voir des tests sur la validit? du plan d'exp?rience
+   
+   time <- ids[[itime]]
+   fixfact <- ids[[ifixfact]]
+   subject <- ids[[isubject]]
+   vd <- ids[[ivd]]
+   
+   # argument of the function instead of re re-running ndim <- defColRes(ids,ifixfact,itime)
+   # nfp : number of main factors + model infos (REML, varSubject) + normality test 
+   nfp <- ndim[1];
+   # ncff number of comparison of the fixed factor
+   nlff <- ndim[2];  ncff <- ndim[3]
+   # nct number of comparison of the time factor
+   nlt <- ndim[4] ; nct <- ndim[5]
+   # nci number of comparison of the interaction
+   nli <- ndim[6];  nci <- ndim[7]
+   # number of all lmer results
+   nresT <- ncff+nct+nci
+   ## initialization of the result vector (1 line)
+   ## 4 * because nresf for : pvalues + Etimates + lower CI + Upper CI
+   res <- data.frame(array(rep(NA,(nfp + 4 * nresT))))
+   colnames(res)[1] <- "resultLM"
+   
+   ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip
+   ### after excluding NA, table function is used to seek subjects with only 1 data
+   ids <- ids[!is.na(ids[[ivd]]),]
+   skip <- length(which(table(ids[[isubject]])==1))
+   
+   if (skip==0) {
+      
+      mfl <- lmer( vd ~ time + fixfact + time:fixfact + (1| subject), ids) # lmer remix
+      
+      # ## NL add   
+      # ###  DEPLACE APRES CALCUL PVALUES AJUSTEES ET NE FAIRE QUE SI AU MOINS 1 FACTEUR SIGNIFICATIF
+      #  if(visu) diagmflF(mfl, title = tit, least.confounded = least.confounded, outlier.limit = outlier.limit)
+      # ## end of NL add
+      
+      rsum <- summary(mfl,ddf = dffOption)
+      ## test Shapiro Wilks on the residus of the model 
+      rShapiro <- shapiro.test(rsum$residuals)
+      raov <- anova(mfl,ddf = dffOption)
+      dlsm1  <- data.frame(difflsmeans(mfl,test.effs=NULL))
+      ddlsm1 <- dlsm1
+      ## save rownames and factor names
+      rn <- rownames(ddlsm1)
+      fn <- ddlsm1[,c(1,2)]
+      ## writing the results on a single line
+      namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
+      namesFactPval <- paste("pvalue ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
+      namesInter <- rownames(ddlsm1)[-c(1:(nct+ncff))]
+      #ncI <- nchar(namesInter)
+      namesEstimate <- paste("estimate ",namesInter)
+      namespvalues <- paste("pvalue ",namesInter)
+      namesFactprinc <- c("pval_time","pval_trt","pval_inter")
+      namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
+      
+      namesFactLowerCI <- paste("lowerCI ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
+      namesLowerCI <- paste("lowerCI ",namesInter,sep="")
+      
+      namesFactUpperCI <- paste("UpperCI ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")
+      namesUpperCI <- paste("UpperCI ",namesInter,sep="")
+      
+      
+      ### lmer results on 1 vector row
+      # pvalue of shapiro Wilks test of the residuals
+      res[1,] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals" 
+      res[2,] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"
+      res[3,] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"
+      ### 3 principal factors pvalues results + shapiro test =>  nfp <- 4
+      res[c((nfp-2):nfp),] <- raov[,6]; rownames(res)[c((nfp-2):nfp)] <- namesFactprinc
+      
+      ####################  Residuals diagnostics for significants variables #########################
+      ### Il at least 1 factor is significant and visu=TRUE NL graphics add to pdf
+      ## ajout JF du passage de la valeur de p-value cutoff
+       if (length(which(raov[,6]<=pvalCutof))>0 & visu == 'yes')  {
+          diagmflF(mfl, title = tit, pvalCutof = pvalCutof, least.confounded = least.confounded, 
+                   outlier.limit = outlier.limit)
+          
+          cat(" Signif ",pvalCutof)
+          
+          
+       }
+      
+      # pvalue of fixed factor comparisons
+      nresf <- nresT
+      res[(nfp+1):(nfp+nct),] <- ddlsm1[c(1:nct),9]
+      res[(nfp+nct+1):(nfp+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),9]
+      rownames(res)[(nfp+1):(nfp+nct+ncff)] <- namesFactPval
+      res[(nfp+nct+ncff+1):(nfp+nresf),] <- ddlsm1[(nct+ncff+1):(nresT),9]
+      rownames(res)[(nfp+nct+ncff+1):(nfp+nresT)] <- namespvalues
+      # Estimate of the difference between levels of factors
+      res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),3]
+      res[(nfp+nresf+nct+1):(nfp+nresf+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),3]
+      rownames(res)[(nfp+nresf+1):(nfp+nresf+nct+ncff)] <- namesFactEstim
+      res[(nfp+nresf+nct+ncff+1):(nfp+2*nresf),] <- ddlsm1[(nct+ncff+1):(nresT),3]
+      rownames(res)[(nfp+nresf+nct+ncff+1):(nfp+2*nresf)] <- namesEstimate
+      # lower CI of the difference between levels of factors
+      nresf <- nresf + nresT
+      res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),7]
+      res[(nfp+nresf+nct+1):(nfp+nresf+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),7]
+      rownames(res)[(nfp+nresf+1):(nfp+nresf+nct+ncff)] <- namesFactLowerCI
+      res[(nfp+nresf+nct+ncff+1):(nfp+2*nresf),] <- ddlsm1[(nct+ncff+1):(nresf),7]
+      rownames(res)[(nfp+nresf+nct+ncff+1):(nfp+nresf+(nresf/2))] <- namesLowerCI
+      # Upper CI of the difference between levels of factors
+      nresf <- nresf + nresT
+      res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),8]
+      res[(nfp+nresf+nct+1):(nfp+nresf+nct+ncff),] <- ddlsm1[(nct+1):(nct+ncff),8]
+      rownames(res)[(nfp+nresf+1):(nfp+nresf+nct+ncff)] <- namesFactUpperCI
+      res[(nfp+nresf+nct+ncff+1):(nfp+nresf+(nresT)),] <- ddlsm1[(nct+ncff+1):(nresT),8]
+      rownames(res)[(nfp+nresf+nct+ncff+1):(nfp+nresf+(nresT))] <- namesUpperCI
+      
+
+   }
+   else
+      ## one of the subject has only one time, subject can't be a random variable
+      ## A repeated measure could be run instead function lme of package nlme, next version
+   {   res[1,] <- NA
+      #cat("impossible computing\n")
+   
+   #  # ## NL add (useless)
+   #    if(visu){
+   #       grid.arrange(ggplot(data.frame()) + geom_point() + xlim(-1, 1) + ylim(-1, 1)+
+   #                    annotate("text", x = 0, y = 0, label = "impossible computing")+
+   #                    xlab(NULL) +  theme(axis.text.x=element_blank(),axis.ticks.x=element_blank())+
+   #                    ylab(NULL) +  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank())+
+   #                    theme(panel.grid.minor = element_blank() ,
+   #                          panel.grid.major = element_blank() ,
+   #                          panel.background = element_rect(fill = "white"))
+   #                 , top = textGrob(tit,gp=gpar(fontsize=40,font=4)))
+   #    
+   #    }
+   # # ## end of NL add
+   
+   }
+   tres <- data.frame(t(res)); rownames(tres)[1] <- nameVar
+   cres <- list(tres,rn, fn)
+   return(cres)
+}
+
+##############################################################################################################
+lmRepeated1FF <- function(ids, ifixfact=0, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]], 
+                          dffOption,pvalCutof=0.05) 
+   {
+   ### function to perform linear mixed model with factor Time + random factor subject
+   ### based on lmerTest package providing functions giving the same results as SAS proc mixed
+   
+   if (!is.numeric(ids[[ivd]]))     {stop("Dependant variable is not numeric")}
+   if (!is.factor(ids[[itime]]))    {stop("Repeated factor is not a factor")}
+   if (!is.factor(ids[[isubject]])) {stop("Random factor is not a factor")}
+   # a ce stade, il faudrait pr?voir des tests sur la validit? du plan d'exp?rience
+   
+   time <- ids[[itime]]
+   subject <- ids[[isubject]]
+   vd <- ids[[ivd]] ## dependant variables (quatitative)
+   
+   # ndim <- defColRes(ids,0,itime)
+   # nfp : nombre de facteurs principaux + model infos + normality test 
+   nfp <- ndim[1]
+   # nct number of comparison of the time factor
+   nlt <- ndim[4] ; nct <- ndim[5]
+   # number of all lmer results
+   nresf <- nct
+   ## initialization of the result vector (1 line)
+   res <- data.frame(array(rep(NA,(nfp+2*nresf))))
+   colnames(res)[1] <- "resultLM"
+   
+   ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip
+   ### after excluding NA, table function is used to seek subjects with only 1 data
+   ids <- ids[!is.na(ids[[ivd]]),]
+   skip <- length(which(table(ids[[isubject]])==1))
+   
+   if (skip==0) {
+      
+      mfl <- lmer( vd ~ time + (1| subject), ids) # lmer remix
+      rsum <- summary(mfl,ddf = dffOption)
+      ## test Shapiro Wilks on the residus of the model 
+      rShapiro <- shapiro.test(rsum$residuals)
+      raov <- anova(mfl,ddf = dffOption)
+      ## Sum of square : aov$'Sum Sq', Mean square : aov$`Mean Sq`, proba : aov$`Pr(>F)`
+      
+      ## Test of all differences estimates between levels as SAS proc mixed. 
+      ## results are in diffs.lsmeans.table dataframe
+      ## test.effs=NULL perform all pairs comparisons including interaction effect
+      dlsm1  <- difflsmeans(mfl,test.effs=NULL)
+      ddlsm1 <- dlsm1$diffs.lsmeans.table
+      
+      ## writing the results on a single line
+      namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct))],sep="")
+      namesFactPval <- paste("pvalue ",rownames(ddlsm1)[c(1:(nct))],sep="")
+      namesFactprinc <- "pval_time"
+      
+      ### lmer results on 1 vector
+      # pvalue of shapiro Wilks test of the residuals
+      res[1,] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals" 
+      res[2,] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"
+      res[3,] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"
+      
+      ### principal factor time pvalue results + shapiro test 
+      res[nfp,] <- raov[,6]; rownames(res)[nfp] <- namesFactprinc
+      # pvalue of fixed factor comparisons
+      res[(nfp+1):(nfp+nct),] <- ddlsm1[c(1:nct),7]
+      rownames(res)[(nfp+1):(nfp+nct)] <- namesFactPval
+      
+      # Estimate of the difference between levels of factors
+      res[(nfp+nresf+1):(nfp+nresf+nct),] <- ddlsm1[c(1:nct),1]
+      rownames(res)[(nfp+nresf+1):(nfp+nresf+nct)] <- namesFactEstim
+   }
+   else
+      ## one of the subject has only one time, subject can't be a random variable
+      ## A repeated measure could be run instead function lme of package nlme, next version
+   {   res[1,] <- NA
+   #cat("traitement impossible\n")
+   }
+   tres <- data.frame(t(res)); rownames(tres)[1] <- nameVar
+   return(tres)
+}
+
+##############################################################################################################
+defColRes <- function(ids, ifixfact, itime) {
+   ## define the size of the result file depending on the numbers of levels of the fixed and time factor.
+   ## Numbers of levels define the numbers of comparisons with pvalue and estimate of the difference.
+   ## The result file also contains the pvalue of the fixed factor, time factor and interaction
+   ## plus Shapiro normality test. This is define by nfp 
+   ## subscript of fixed factor=0 means no other fixed factor than "time"
+   if (ifixfact>0){
+      nfp <- 6 # shapiro+time+fixfact+interaction+ others....
+      time <- ids[[itime]]
+      fixfact <- ids[[ifixfact]]
+      
+      cat("\n levels fixfact",levels(fixfact))
+      cat("\n levels time",levels(time))
+      
+      # ncff number of comparisons of the fixed factor (nlff number of levels of fixed factor)
+      nlff <- length(levels(fixfact)); ncff <- (nlff*(nlff-1))/2
+      # nct number of comparison of the time factor (nlt number of levels of time factor)
+      nlt <- length(levels(time)); nct <- (nlt*(nlt-1))/2
+      # nci number of comparison of the interaction
+      nli <- nlff*nlt;  nci <- (nli*(nli-1))/2
+      ndim <- c(NA,NA,NA,NA,NA,NA,NA)
+      
+      ndim[1] <- nfp   # pvalues of fixed factor, time factor and interaction (3columns) and shapiro test pvalue
+      ndim[2] <- nlff  # number of levels of fixed factor
+      ndim[3] <- ncff  # number of comparisons (2by2) of the fixed factor
+      ndim[4] <- nlt   # number of levels of time factor
+      ndim[5] <- nct   # number of comparisons (2by2) of the time factor
+      ndim[6] <- nli   # number of levels of interaction
+      ndim[7] <- nci   # number of comparisons (2by2) of the interaction
+      
+   } 
+   else {
+      nfp <- 4 # shapiro+time
+      time <- ids[[itime]]
+      # nct number of comparison of the time factor
+      nlt <- length(levels(time)); nct <- (nlt*(nlt-1))/2
+      ndim <- c(NA,NA,NA,NA,NA,NA,NA)
+      
+      ndim[1] <- nfp   # pvalues of time factor and shapiro test pvalue
+      ndim[4] <- nlt   # number of levels of time factor
+      ndim[5] <- nct   # number of comparisons (2by2) of the time factor
+   }
+   return(ndim)
+}
+
+##############################################################################################################
+lmixedm <- function(datMN, 
+                    samDF,
+                    varDF,
+                    fixfact, time, subject,
+                    logtr = "none", 
+                    pvalCutof = 0.05,
+                    pvalcorMeth = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")[7],
+                    dffOption,
+                    visu = "no",
+                    least.confounded = FALSE,
+                    outlier.limit = 3,
+					     pdfC,
+					     pdfE
+                    )
+   {
+   sampids <- samDF
+   dataMatrix <- datMN
+   varids <- varDF
+   
+   options("scipen" = 50, "digits" = 5)
+   pvalCutof <- as.numeric(pvalCutof)
+   
+   cat("\n dff computation method=",dffOption)
+   ### Function running lmer function on a set of variables described in 
+   ### 3 different dataframes as used by W4M
+   ### results are merge with the metadata variables varids
+   ### ifixfact, itime, isubject are subscripts of the dependant variables
+   if (fixfact=="none")  ifixfact <-0 else ifixfact <- which(colnames(sampids)==fixfact)
+   itime    <- which(colnames(sampids)==time)
+   isubject <- which(colnames(sampids)==subject)
+   
+    #lmmds <- dataMatrix[,-1]
+   
+   lmmds <- dataMatrix
+   if (logtr!="log10" & logtr!="log2") logtr <- "none"
+   if (logtr=="log10") lmmds <- log10(lmmds+1)
+   if (logtr== "log2") lmmds <- log2(lmmds+1)
+   
+   #idsamp <- dataMatrix[,1]
+   #lmmds <- t(lmmds)
+   dslm <- cbind(sampids,lmmds)
+
+   nvar <- ncol(lmmds); firstvar <- ncol(sampids)+1; lastvar <- firstvar+ncol(lmmds)-1
+
+   dslm[[ifixfact]] <- factor(dslm[[ifixfact]])
+   dslm[[itime]]    <- factor(dslm[[itime]])  
+   dslm[[isubject]] <- factor(dslm[[isubject]])
+   ## call defColres to define the numbers of test and so the number of columns of results
+   ## depends on whether or not there is a fixed factor with time. If only time factor ifixfact=0
+   if (ifixfact>0) {
+      ndim <- defColRes(dslm[,c(ifixfact,itime)],ifixfact=1,itime=2)
+      nColRes <- ndim[1]+(4*(ndim[3]+ndim[5]+ndim[7]))
+      firstpval <- ndim[1]-2
+      lastpval <- ndim[1]+ndim[3]+ndim[5]+ndim[7]
+  } else 
+   {
+      ndim <- defColRes(dslm[,itime],ifixfact=0,itime=1) 
+      nColRes <- ndim[1]+(2*(ndim[5]))
+      firstpval <- ndim[1]
+      lastpval <- ndim[1]+ndim[5]
+   }
+   ## initialisation of the  result file 
+   resLM <- data.frame(array(rep(NA,nvar*nColRes),dim=c(nvar,nColRes)))
+   rownames(resLM) <- rownames(varids)
+
+   ###############  test ecriture dans pdf
+   if(visu == "yes") {
+      pdf(pdfC, onefile=TRUE, height = 15, width = 30)
+      par(mfrow=c(1,3))
+   }
+   ###############  fin test ecriture dans pdf
+   ## pour test : lastvar <- 15
+   cat("\n pvalCutof ", pvalCutof)
+   
+   for (i in firstvar:lastvar) {
+      
+      ## NL modif
+      cat("\n[",colnames(dslm)[i],"] ")
+      ## end of NL modif
+      
+      subds <- dslm[,c(ifixfact,itime,isubject,i)]
+      
+      ## NL modif
+      tryCatch({
+         if (ifixfact>0)
+            reslmer <- lmRepeated2FF(subds,ifixfact=1,itime=2,isubject=3, ivd=4, ndim=ndim, visu = visu, 
+                                     tit = colnames(dslm)[i], pvalCutof=pvalCutof,
+                                     dffOption=dffOption,least.confounded = least.confounded, 
+                                     outlier.limit = outlier.limit) 
+         else 
+            reslmer <- lmRepeated1FF(subds,ifixfact=0,1,2, ivd=3, ndim=ndim, pvalCutof=pvalCutof,dffOption)
+         ## end of NL modif
+         resLM[i-firstvar+1,] <- reslmer[[1]]
+      }, error=function(e){cat("ERROR : ",conditionMessage(e), "\n");})
+      if (i==firstvar) {
+         colnames(resLM) <- colnames(reslmer[[1]])
+         labelRow <- reslmer[[2]]
+         factorRow <- reslmer[[3]]
+      }
+   }
+   ## for debug : ifixfact=1;itime=2;isubject=3; ivd=4;tit = colnames(dslm)[i]; ids <- subds
+
+   
+   ## NL add
+   if(visu == "yes") dev.off()
+   ## end of NL add
+   
+   ## pvalue correction with p.adjust library multtest
+   ## Possible methods of pvalue correction
+   AdjustMeth <- c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr","none")
+   if (length(which(pvalcorMeth == AdjustMeth))==0) pvalcorMeth <- "none"
+
+   if (pvalcorMeth !="none") {
+      for (k in firstpval:lastpval){
+         resLM[[k]]=p.adjust(resLM[[k]], method=pvalcorMeth, n=dim(resLM[k])[[1]])
+         
+      }
+   }
+
+   ## for each variables, set pvalues to NA and estimates = 0 when pvalue of factor > pvalCutof value define by user
+   if (ifixfact>0) {
+      ## time effect
+      resLM[which(resLM[,firstpval]> pvalCutof),c((lastpval+1):(lastpval+ndim[5]))] <- 0
+      resLM[which(resLM[,firstpval]> pvalCutof),c((ndim[1]+1):(ndim[1]+ndim[5]))] <- NA
+   ## treatment effect
+      resLM[which(resLM[,firstpval+1]> pvalCutof),c((lastpval+ndim[5]+1):(lastpval+ndim[5]+ndim[3]))] <- 0
+      resLM[which(resLM[,firstpval+1]> pvalCutof),c((ndim[1]+ndim[5]+1):(ndim[1]+ndim[5]+ndim[3]))] <- NA
+   ## interaction effect
+      resLM[which(resLM[,firstpval+2]> pvalCutof),c((lastpval+ndim[5]+ndim[3]+1):(lastpval+ndim[5]+ndim[3]+ndim[7]))] <- 0
+      resLM[which(resLM[,firstpval+2]> pvalCutof),c((ndim[1]+ndim[5]+ndim[3]+1):(ndim[1]+ndim[5]+ndim[3]+ndim[7]))] <- NA
+   } else {
+      ## time effect only
+      resLM[which(resLM[,firstpval]> pvalCutof),c((lastpval+1):(lastpval+ndim[5]))] <- 0
+      resLM[which(resLM[,firstpval]> pvalCutof),c((firstpval+1):(firstpval+ndim[5]))] <- NA
+   }
+   
+   ## for each variable, estimates plots are performed if at least one factor is significant after p-value correction
+   pdf(pdfE, onefile=TRUE, height = 15, width = 30)
+   #par(mfrow=c(2,2))
+   
+   ## for each variable (in row)   
+   for (i in 1:nrow(resLM)) {
+      #cat("\n",rownames(resLM)[i])
+      ## if any main factor after p-value correction is significant -> plot estimates and time course
+      if (length(which(resLM[i,c(4:6)]<pvalCutof))>0) {
+         
+         ## Plot of time course by fixfact : data prep with factors and quantitative var to be plot
+         subv <- dslm[,colnames(dslm)==rownames(resLM)[i]]
+         subds <- data.frame(dslm[[ifixfact]],dslm[[itime]], dslm[[isubject]],subv)
+         #colnames(subds) <- c(colnames(dslm)[ifixfact],colnames(dslm)[itime],colnames(dslm)[isubject],rownames(resLM)[i] <- rownames(resLM)[i] )
+         libvar <- c(fixfact,time,subject)
+         colnames(subds) <- c(libvar,rownames(resLM)[i])
+         
+         ## Plot of estimates with error bars for all fixed factors and interaction
+         rddlsm1 <- t(resLM[i,])
+         pval <- rddlsm1[substr(rownames(rddlsm1),1,6)=="pvalue"]
+         esti <- rddlsm1[substr(rownames(rddlsm1),1,6)=="estima"]
+         loci <- rddlsm1[substr(rownames(rddlsm1),1,6)=="lowerC"]
+         upci <- rddlsm1[substr(rownames(rddlsm1),1,6)=="UpperC"]
+         rddlsm1 <- data.frame(pval,esti,loci,upci,factorRow)
+         colnames(rddlsm1) <- c("p.value","Estimate","Lower.CI","Upper.CI",colnames(factorRow))
+         rownames(rddlsm1) <- labelRow
+
+         ## function for plotting these 2 graphs
+         plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)
+    
+      }
+   }
+   dev.off()
+   
+   ## return result file with pvalues and estimates (exclude confidence interval used for plotting)
+   iCI <- which(substr(colnames(resLM),4,7)=="erCI")
+   resLM <- resLM[,-iCI] 
+   resLM <- cbind(varids,resLM)
+   return(resLM)
+}
+
+