R/mixed_models.R

   1
   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 if(length(args)==0){
   9     print("No arguments supplied.")
  10     ##supply default values
  11     datafile=''
  12     paramfile=''
  13 }else{
  14     for(i in 1:length(args)){
  15          eval(parse(text=args[[i]]))
  16     }
  17 }
  18
  19 workdir = dirname(datafile);
  20 setwd(workdir);
  21
  22 library(lme4)
  23 library(lmerTest)
  24 #library(emmeans)
  25 library(effects)
  26 library(phenoAnalysis)
  27 library(stringr)
  28 library(dplyr)
  29
  30 pd = read.csv(datafile, sep="\t")
  31 source(paramfile)  # should give us dependent_variable and the model
  32
  33 pd$studyYear = as.factor(pd$studyYear)
  34 print(paste("MODEL :", model))
  35
  36 print(head(pd))
  37 dependent_variable = gsub(" ", "\\.", dependent_variable) # replace space with "." in variable name
  38 dependent_variable = gsub("\\|", "\\.", dependent_variable) # replace | with .
  39 dependent_variable = gsub("\\:", "\\.", dependent_variable)
  40 dependent_variable = gsub("\\-", "\\.", dependent_variable)
  41 dependent_variable = gsub("\\/", "\\.", dependent_variable)
  42 print(paste("Dependent variable : ", dependent_variable))
  43 model_string = paste(dependent_variable, '~', model)
  44
  45 print(paste('MODEL STRING:', model_string));
  46 model = lmer(as.formula(model_string), data=pd)
  47
  48
  49 #model_summary = summary(allEffects(model,se=T))
  50
  51 pdout = model
  52 print(pdout)
  53
  54 genotypeEffectType = as.vector(str_match(model_string, '1\\|germplasmName'))
  55 genotypeEffectType = ifelse(is.na(genotypeEffectType), 'fixed', 'random')
  56 print(paste('modeling genotypes as: ', genotypeEffectType))
  57
  58 adjusted_means = getAdjMeans(modelOut=model,
  59     traitName=dependent_variable,
  60     genotypeEffectType=genotypeEffectType,
  61     adjMeansVariable='germplasmName')
  62
  63 print(head(adjusted_means))
  64
  65 outfile = paste(datafile, ".results", sep="")
  66 print(outfile)
  67 print(model)
  68 #print(model_summary)
  69 print(colnames(model))
  70 print(ranef(model))
  71 print(adjusted_means)
  72 sink(outfile)
  73 #write.table(ranef(model)$germplasmName)
  74
  75 write.table(select(adjusted_means, 'germplasmName', dependent_variable), row.names = F)
  76 sink();