R/mixed_models.R

   1 ## File was edited in ESS mode in emacs and some indents changed. sorry.
   2 ##First read in the arguments listed at the command line
   3 args=commandArgs(TRUE)
   4
   5 ##args is now a list of character vectors
   6 ## First check to see if arguments are passed.
   7 ## Then cycle through each element of the list and evaluate the expressions.
   8
   9 if(length(args)==0){
  10     print("No arguments supplied.")
  11     ##supply default values
  12     datafile=''
  13     paramfile=''
  14     fixed_factors<-c()
  15     random_factors<-c()
  16
  17 } else {
  18     for(i in 1:length(args)){
  19         print(paste("Processing arg ", args[[i]]));
  20         eval(parse(text=args[[i]]))
  21     }
  22 }
  23
  24
  25 workdir = dirname(datafile);
  26 setwd(workdir);
  27
  28 library(lme4)
  29 #library(lmerTest)
  30 library(emmeans)
  31 library(effects)
  32 library(stringr)
  33 library(dplyr)
  34
  35 pd = read.csv(datafile, sep="\t")
  36
  37 source(paramfile)  # should give us dependent_variable and the model
  38
  39 trait = dependent_variables
  40 pd$studyYear = as.factor(pd$studyYear)
  41 print(paste("MODEL :", model))
  42 print(paste("FIXED FACTORS: ", fixed_factors));
  43 print(paste("RANDOM FACTORS: ", random_factors));
  44 print(head(pd))
  45
  46 BLUE = as.data.frame(unique(pd$germplasmName))
  47 colnames(BLUE) = "germplasmName"
  48 adjustedBLUE = BLUE
  49
  50 BLUP = BLUE
  51 adjusted_means = BLUE
  52
  53 for(i in 1:length(trait)){
  54
  55     print(paste("processing trait", trait[i]))
  56     dependent_variables = trait[i]
  57     dependent_variables = gsub(" ", "\\.", dependent_variables) # replace space with "." in variable name
  58     dependent_variables = gsub("\\|", "\\.", dependent_variables) # replace | with .
  59     dependent_variables = gsub("\\:", "\\.", dependent_variables)
  60     dependent_variables = gsub("\\-", "\\.", dependent_variables)
  61     dependent_variables = gsub("\\/", "\\.", dependent_variables)
  62
  63     print(paste("Dependent variables : ", dependent_variables))
  64
  65     model_string = paste(dependent_variables, '~', model)
  66
  67     print(paste('MODEL STRING:', model_string));
  68     #mixmodel = lmer(as.formula(model_string), data=pd)
  69
  70     #model_summary = summary(allEffects(model,se=T))
  71
  72     #pdout = model
  73     #print(pdout)
  74
  75     genotypeEffectType = as.vector(str_match(model_string, '1\\|germplasmName'))
  76     genotypeEffectType = ifelse(is.na(genotypeEffectType), 'fixed', 'random')
  77     print(paste('modeling genotypes as: ', genotypeEffectType))
  78
  79     if (genotypeEffectType=="random") {
  80
  81         mixmodel = lmer(as.formula(model_string), data=pd)
  82
  83         print("---------")
  84         print(mixmodel)
  85         print("---------")
  86
  87         res <- (ranef(mixmodel)$germplasmName)
  88
  89         blup <- res%>%mutate("germplasmName" = rownames(res))
  90         names(blup)[1] = trait[i]
  91         blup <- blup[,c("germplasmName",trait[i])]
  92
  93         blup <- as.data.frame(blup)
  94
  95         BLUP <- merge(x = BLUP, y = blup, by ="germplasmName", all=TRUE)
  96
  97         adj = coef(mixmodel)$germplasmName
  98
  99         print(paste("adj", adj));
 100
 101         adj = adj[1]  # keep only one column
 102         adj_means = adj%>%mutate("germplasmName" = rownames(adj))
 103         names(adj_means)[1] = trait[i]
 104         adj_means = as.data.frame(adj_means)
 105
 106         adjusted_means =  merge(x = adjusted_means, y = adj_means, by ="germplasmName", all=TRUE)
 107
 108     } else {
 109
 110         mixmodel = lmer(as.formula(model_string), data=pd)
 111
 112         # compute adjusted blues
 113         #
 114         adj <- summary(lsmeans(mixmodel, "germplasmName"))
 115         blue <- adj[c("germplasmName", "lsmean")]
 116         colnames(blue)[2] = trait[i]
 117         blueadj = as.data.frame(blue)
 118         adjustedBLUE =  merge(x = adjustedBLUE, y = blueadj, by ="germplasmName", all=TRUE)
 119
 120         #Computing fixed effects
 121         #feff <- (fixef(mixmodel)$germplasmName)
 122         feff<-data.frame(coef(summary(mixmodel))[ , "Estimate"])
 123         rownames(feff) <- blue$germplasmName
 124         colnames(feff) <- trait[i]
 125
 126         fixedeff <- feff%>%mutate("germplasmName" = rownames(feff))
 127         fixedeff<- fixedeff[,c("germplasmName", trait[i])]
 128         fixedeff<-as.data.frame(fixedeff)
 129
 130         #file with fixed effect
 131         print(fixedeff)
 132
 133         #file with blues
 134         print(blue)
 135
 136         BLUE = merge(x = BLUE, y = fixedeff, by="germplasmName", all=TRUE)
 137
 138
 139     }
 140 }
 141
 142 # save result files
 143 #
 144 # for random effects: file.BLUPs and file.adjustedBLUPs
 145 # for fixed effects: file.BLUEs and file.adjustedBLUEs
 146 #
 147 if (genotypeEffectType=="random") {
 148     outfile_blup = paste(datafile, ".BLUPs", sep="");
 149     sink(outfile_blup)
 150     write.table(BLUP, quote=F, sep='\t', row.names=FALSE)
 151     sink();
 152     outfile_adjmeans = paste(datafile, ".adjustedBLUPs", sep="")
 153     sink(outfile_adjmeans)
 154     write.table(adjusted_means, quote=F , sep='\t', row.names=FALSE)
 155     sink();
 156
 157 } else {   #fixed
 158     outfile_blue = paste(datafile, ".BLUEs", sep="")
 159     sink(outfile_blue)
 160     write.table(BLUE, quote=F , sep='\t', row.names=FALSE)
 161     sink();
 162     outfile_adjblues = paste(datafile, ".adjustedBLUEs", sep="")
 163     sink(outfile_adjblues)
 164     write.table(adjustedBLUE, quote=F , sep='\t', row.names=FALSE)
 165     sink();
 166 }