#this script takes models for male bee project, computes estimated median odds ratios, and then plots the MORs. #model fitting takes a really long time! here are the "allfit" model objects, which are lists of merMod objects load("models/phenology.rda") ###Now using rtime[[5]]. This wrap for neldermead gives low likelihood load("models/morecrossed_allfit.rda") # for this one I use the first one fm3.all[[2]] #loAD DATA library(merTools) library(lme4) library(parallel) library(plyr) library(tidyverse) library(ggstance) # # #get approximate MORs # MOR.fm3<-data.frame(RE=as.data.frame(VarCorr(fm3.all[[2]]))$grp, MOR=exp(sqrt(2*as.data.frame(VarCorr(fm3.all[[2]]))$vcov)*0.6745)) # MOR.fm3 #just the OR for the fixed effect with round 1 as reference level # OR.fe.fm2<-exp(as.data.frame(fixef(fm2.all[[1]]))) # write.csv(MOR.fm3, "data/no_phenology_MOR.csv") # write.csv(MOR.fm2, "data/phenology_MOR.csv") # write.csv(OR.fe.fm2, "data/phenology_fixed_OR.csv") #now for something more sophisticated,a bootstrapping for MOR with 95% CrI #Milan Seth (https://stats.stackexchange.com/users/188609/milan-seth), How to implement credible 95% interval for median odds ratio using JAGS?, URL (version: 2017-12-14): https://stats.stackexchange.com/q/318912 #fm1<-fm2.all[[1]] ############################################ ########### SET the model MORcalc<-function(fm1){ # fm1<-fm3.all[[2]] t<- lme4::ranef(fm1,condVar = TRUE) Sys.time()->start; no_cores <- detectCores() - 1 cl <- makeCluster(no_cores) clusterExport(cl=cl, varlist=ls(),envir = environment()) MOR_CrI<-lapply(t, FUN=function(x){var<- as.numeric(unlist(attr(x,"postVar"))) # MOR_CrI_phenology<-parLapply(cl, t, fun=function(x){var<- as.numeric(unlist(attr(x,"postVar"))) est<-as.numeric(unlist(x)) ### bootstrap #### iterate over replicates (this is what actually needs to be parallel to speed up) mor_boot<-parLapply(cl, 1:1000, fun=function(i){ ### draw vector of area random effects from normal drw<- rnorm(n = length(est),mean = est,sd = sqrt(var)) ### create data frame with all possible pairs s<- combn(drw,2) #### estimate MOR and save in output morboot<-median(exp(abs(s[1,]- s[2,]))) return(morboot) }) mor_boot<-unlist(mor_boot) ### bootstrap median and 95% CI q<-quantile(mor_boot, c(.025,.5,.975), na.rm=T) return(q) } ) stopCluster(cl) print(Sys.time()-start); return(MOR_CrI) } implicit<-cbind(ldply(MORcalc(fm3.all[[2]])), mod=rep("implicit",6)) names(implicit)<-c("term", "lcl", "median", "ucl", "mod") phenology<-cbind(ldply(MORcalc(rtime[[5]])), mod=rep("phenology",12)) names(phenology)<-c("term", "lcl", "median", "ucl", "mod") write.csv(implicit,"data/updated_implicit_MORs.csv", row.names = F) write.csv(phenology, "data/phenology_MOR.csv", row.names=F) implicit<-read.csv("data/updated_implicit_MORs.csv") phenfull<-read.csv("data/phenology_MOR.csv") names(phenfull)<-c("term", "lcl", "median", "ucl", "mod") names(implicit)<-c("term", "lcl", "median", "ucl", "mod") # #combine the MORs into a single df, fix idiotic name orders so that it looks right is. mors<-rbind(implicit[,-6], phenfull) # mors$term<-as.character(mors$term) mors[mors$term=="bee:site", "term"]<-"site:bee" mors[mors$term=="plant_code:site", "term"]<-"site:plant_code" mors[,"term"]<-sub("plant_code", "flower", mors[,"term"]) mors[,"term"]<-sub("bee", "bee", mors[,"term"]) #this is generating ghost levels intentionally so that it's easy to plot them with space between terms in the right spots. #create a column for the plotting symbol mors$pshape<-rep(25,18) mors$pshape[grep("bee", as.character(as.factor(mors$term)))]<-24 mors$pshape[grep("flower", as.character(as.factor(mors$term)))]<-23 mors<-rbind(mors, c("", NA, NA, NA,NA, NA), c(" ", NA, NA, NA, NA,NA), c(" ", NA, NA,NA,NA, NA),c(" ", NA, NA,NA,NA, NA),c(" ", NA, NA,NA,NA, NA),c(" ", NA, NA,NA,NA, NA)) mors$lcl<-as.numeric(mors$lcl) mors$ucl<-as.numeric(mors$ucl) mors$median<-as.numeric(mors$median) #save the data to avoid rerunning (if model fit/data haven't changed) # write.csv(mors, "data/mors.csv") mors<-read.csv("data/mors.csv") mors[,"term"]<-ordered(mors[,"term"], levels = rev(c("bee", "flower", "site", "", "flower:bee", "site:bee", "site:flower"," ", "sampling_round", " ","bee:sampling_round", "flower:sampling_round", "site:sampling_round", "site:bee:sampling_round", "site:flower:sampling_round")) ) ################################################################# ### plot model estimates cairo_ps("figures/MOR_Revised.eps",width=7.3, height=4.5) mors %>% filter(is.na(term)==F) %>% ggplot(aes(median, term, color=mod))+ geom_vline(xintercept = 1, linetype="dashed")+ geom_errorbarh(aes (xmin=lcl, xmax=ucl, color=mod), position=position_dodgev(height=0.3), width=0)+ geom_point(position=position_dodgev(height=0.3), size=2)+ labs(x="median odds ratio", y="")+ scale_x_log10(breaks=c(1,2,4,8))+ theme_classic()+ scale_y_discrete(labels=rev(c("bee", "flower", "site", "", "flower : bee", "site : bee", "site : flower","", "round","", "bee : round", "flower : round", "site : round", "site : bee : round", "site : flower : round")))+ theme(text=element_text(size=15),axis.text=element_text(color="black"),axis.ticks.y=element_blank(), legend.position = "none") dev.off() ############################## #what is correlation between plant effects in m1 and m2? cor(ranef(rtime[[5]])$plant_code, ranef(fm3.all[[2]])$plant_code, method="spearman") cor(ranef(rtime[[5]])$plant_code, ranef(fm3.all[[2]])$plant_code, method="pearson")