<commit_before>
mod.evaluation <- function(yname,D = "yearly.Te.10C.q",R = "NDVI",N = "n.birds",
                           ordered.fac.treatment = c("as.num","as.binomial"),
                           centering = c('none','CGM','CWC'),log.D = F, log.R = F,
                           standardisation = F,year2excl = NA){
  #aim: evaluation of mixed models and quantile regression for Extreme limits" Hummer project. 
  #The model to be calibrateed has the response variable (yname) and 2 predictors (D) and (R) and its interactions.
  #This model is then compared to all its sub-models.
  #All models contain random (intercept-)effect associated with site/location 'LOC'
  #
  #Arguments
  #yname: name of the response variable (e.g. as named in an attached data frame)
  #D:     name of the predictor variable representing physiological demand (e.g. minimum annal temperature)
  #R:     name of the predictor variable representing resource availability (e.g. NDVI)
  #N:     a variable by which observations are to be weighted in the model (e.g. the nr of birds sampled)
  #ordered.fac.treatment:
  #       if the response variable is an ordered factor, how should it be treated? 
  #       Options are "as.num" which treats the response variable as if it were numeric, and
  #       "as.binomial" which treat it as a bernouilli experiment.
  #centering: Should the predictor variable be centered. Options are no centering "none", "GCM" which centers
  #       the data around the grand mean, and "CWC" which centers the data around group means (both sensu EndersTofighi2007PhsychologicalMethods)
  #       Also see: http://mlrv.ua.edu/2009/vol35_1/Robinson_Schumacke_rproof.pdf
  #log.D: Should the 'D' predictor variable be log-transformed?
  #log.R: Should the 'D' predictor variable be log-transformed?
  #standardisation: should predictor variables be standardized (i.e. divided by ther standard deviation)
  #year2excl: specify any years that need to be excluded from the model
  
  
  require(lme4)
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # get the variables
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  nbirds <- NULL
  if(!is.null(N)){nbirds <- get(N)
                  cat("samples will be weighted by 'n.birds' in model calibration\n")}else{
                    cat("all samples will be equally weighted calibration\n")} 
  
  D <- get(D)
  R <- get(R)
  mod.name <- c('mod.DR','mod.D','mod.R','mod.0','mod.D_R')
  y <- get(yname)  
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # do you want a particular year excluded?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
  if (!is.na(year2excl)){cat("Excluding year(s): ",year2excl,"\n")
                         y<-y[yr!=year2excl]
                         D<-D[yr!=year2excl]
                         R<-R[yr!=year2excl]
                         LOC.<-LOC[yr!=year2excl]
                         nbirds. <- nbirds[yr!=year2excl]}else{LOC. <- LOC ; nbirds. <- nbirds}
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # do you want any of the predictors to be log transformed?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  if (log.D){cat("DEMAND variable will be log-transformed\n");D<-log(D)+1}
  if (log.R){cat("RESOURCE variable will be log-transformed\n");R<-log(R)+1}
  cat("\nCorrelation between D and R (after transform and or subset):\n")
  print(cor.test(D,R))
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # is the y-variable an ordered factor?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  if ((is.factor(y[1]))&(ordered.fac.treatment[1] == c("as.num"))){
    cat(yname, " is an ordered factor. options are to treat as numeric or in a binomial fashion.\n")    
    cat(yname, " will be converted to numeric \n")
    y <- as.numeric(y)}
  
  if ((is.factor(y[1]))&(ordered.fac.treatment[1] == c("as.binomial"))){
    cat("\n",yname, " is an ordered factor. options are to treat as numeric or in a binomial fashion !!!\n")
    cat(yname, " is an ordered factor and will be converted to a count table !!!\n")
    y <- ordered.var.2.sim.count.table(y)}
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # apply centering and standardization ####
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  CGM <- function(x){(x-mean(x,na.rm=T))}#/sd(x)}
  CWC <- function(x,LOC){mns.x <- tapply(x,LOC,mean,na.rm=T)
                         mns.x <- mns.x[match(LOC,names(mns.x))]
                         mns.x <- x - mns.x
                         return(mns.x)}  
  if (centering == 'CGM'){D <- CGM(D) ; R <- CGM(R)
                          cat('\nPredictor variables were centered at the grand mean (CGM)\n') }
  if (centering == 'CWC'){
    D <- CWC(x=D,LOC=LOC.)
    R <- CWC(x=R,LOC=LOC.)
    
    cat('\nPredictor variables were centered within clusters (CWC)\n') 
  }
  
  if(standardisation){
    cat('\nPredictor variables were standardised using grand sd\n')
    D <- D/sd(D,na.rm=T)
    R <- R/sd(R,na.rm=T)
  }
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # plot the predictor variables
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
  win.graph(w=9,h=4);par(mfrow=c(1,3))
  hist(unique(D),col="light grey");hist(unique(R),col="light grey");plot(R,D,pch=16)
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # calibrate the main model and it's sub-models
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
  if ((ordered.fac.treatment[1] == c("as.binomial"))& (!is.null(dim(y)))){
    #for models with ordinal response variables treated as binomial counts  
    mod.DR <- glmer(formula = y ~ D*R + (1|LOC.),weights=nbirds.,family="binomial")
    mod.D_R <- glmer(formula = y ~ D + R + (1|LOC.),weights=nbirds.,family="binomial")
    mod.D <- glmer(formula = y ~ D + (1|LOC.),weights=nbirds.,family="binomial")
    mod.R <- glmer(formula = y ~ R + (1|LOC.),weights=nbirds.,family="binomial")
    mod.0 <- glmer(formula = y ~ 1 + (1|LOC.),weights=nbirds.,family="binomial")
    
  }else{
    #for all other models
    mod.DR <- lmer(formula = y ~ D*R + (1|LOC.),weights=nbirds.)
    mod.D_R <- lmer(formula = y ~ D + R + (1|LOC.),weights=nbirds.)
    mod.D <- lmer(formula = y ~ D + (1|LOC.),weights=nbirds.)
    mod.R <- lmer(formula = y ~ R + (1|LOC.),weights=nbirds.)
    mod.0 <- lmer(formula = y ~ 1 + (1|LOC.),weights=nbirds.)
  }  
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #in the D*R model, under what conditions are the D (R) terms negative?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  source('D:\\Hummer_code_transfer\\Plot_conditional_effects_interaction_models.r')
  cat("Conditional effects of predictors in full model, i.e. including R*D interaction\n")
  win.graph(w=8,h=8);par(mfcol=c(2,2))
  
  r1 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,plott=F)
  r2 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,switch.vars = T,plott=F)
  
  r12 <- c(-1,1)*max(abs(c(r1,r2)),na.rm=T)
    r1 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,ylim.=r12)
  r2 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,ylim.=r12,switch.vars = T)
    
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # get the AIC  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  require(bbmle)
  
  #since 'n' is not defined in a mixed model, calculate a liberal AIC (mod.AIC.CON)
  #assuing n is the nr of data points (nobs.HI), and a liberal AIC (mod.AIC.LIB) assuming n is the nr of groupd (nobs.LO)
  nobs.HI <- nrow(mod.DR@frame)
  nobs.LO <- nlevels(mod.DR@frame$LOC)
  
  mod.AIC.LIB <- ICtab(mod.DR,mod.D,mod.R,mod.0,mod.D_R,sort=F,type="AICc",nobs=nobs.HI)
  mod.AIC.CON <- ICtab(mod.DR,mod.D,mod.R,mod.0,mod.D_R,sort=F,type="AICc",nobs=nobs.LO)
  
  mod.AIC <- cbind.data.frame(round(mod.AIC.LIB$dAICc,1),round(mod.AIC.CON$dAICc,1),mod.AIC.LIB$df)
  colnames(mod.AIC) <- c('dAICc.lib','dAICCc.con','df')
  rownames(mod.AIC) <- c('mod.DR','mod.D','mod.R','mod.0','mod.D_R')
  
  #calculate Akaike weights sensu http://www.planta.cn/forum/files_planta/2006_nobbs_ecology_153.pdf
  Akaike.weight <- function(AIC.vec){AIC.vec <- AIC.vec-min(AIC.vec,na.rm=T)
                                     Akai.weight.denom <- sum(exp(-2*AIC.vec))
                                     Akai.weight <- (exp(-2*AIC.vec))/Akai.weight.denom
                                     return(Akai.weight)
  }
  A.weight <- round(Akaike.weight(mod.AIC$dAICc.lib),2)
  mod.AIC <- cbind.data.frame(mod.AIC,A.weight)<commit_msg>moved import and source statments to top of code - will fail early, rather than later<commit_after># code for evaluating mixed models and quantile regresstion for the "extreme limits" hummingbird paper
# modified from <NAME>

# import packages
require(lme4)
require(bbmle)

# source data and r files
source('/Users/sarah/Documents/GitHub/extreme_limits/Plot_conditional_effects_interaction_models.r')


mod.evaluation <- function(yname,D = "yearly.Te.10C.q",R = "NDVI",N = "n.birds",
                           ordered.fac.treatment = c("as.num","as.binomial"),
                           centering = c('none','CGM','CWC'),log.D = F, log.R = F,
                           standardisation = F,year2excl = NA){
  #aim: evaluation of mixed models and quantile regression for Extreme limits" Hummer project. 
  #The model to be calibrateed has the response variable (yname) and 2 predictors (D) and (R) and its interactions.
  #This model is then compared to all its sub-models.
  #All models contain random (intercept-)effect associated with site/location 'LOC'
  #
  #Arguments
  #yname: name of the response variable (e.g. as named in an attached data frame)
  #D:     name of the predictor variable representing physiological demand (e.g. minimum annal temperature)
  #R:     name of the predictor variable representing resource availability (e.g. NDVI)
  #N:     a variable by which observations are to be weighted in the model (e.g. the nr of birds sampled)
  #ordered.fac.treatment:
  #       if the response variable is an ordered factor, how should it be treated? 
  #       Options are "as.num" which treats the response variable as if it were numeric, and
  #       "as.binomial" which treat it as a bernouilli experiment.
  #centering: Should the predictor variable be centered. Options are no centering "none", "GCM" which centers
  #       the data around the grand mean, and "CWC" which centers the data around group means (both sensu EndersTofighi2007PhsychologicalMethods)
  #       Also see: http://mlrv.ua.edu/2009/vol35_1/Robinson_Schumacke_rproof.pdf
  #log.D: Should the 'D' predictor variable be log-transformed?
  #log.R: Should the 'D' predictor variable be log-transformed?
  #standardisation: should predictor variables be standardized (i.e. divided by ther standard deviation)
  #year2excl: specify any years that need to be excluded from the model
  
  

  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # get the variables
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  nbirds <- NULL
  if(!is.null(N)){nbirds <- get(N)
                  cat("samples will be weighted by 'n.birds' in model calibration\n")}else{
                    cat("all samples will be equally weighted calibration\n")} 
  
  D <- get(D)
  R <- get(R)
  mod.name <- c('mod.DR','mod.D','mod.R','mod.0','mod.D_R')
  y <- get(yname)  
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # do you want a particular year excluded?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
  if (!is.na(year2excl)){cat("Excluding year(s): ",year2excl,"\n")
                         y<-y[yr!=year2excl]
                         D<-D[yr!=year2excl]
                         R<-R[yr!=year2excl]
                         LOC.<-LOC[yr!=year2excl]
                         nbirds. <- nbirds[yr!=year2excl]}else{LOC. <- LOC ; nbirds. <- nbirds}
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # do you want any of the predictors to be log transformed?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  if (log.D){cat("DEMAND variable will be log-transformed\n");D<-log(D)+1}
  if (log.R){cat("RESOURCE variable will be log-transformed\n");R<-log(R)+1}
  cat("\nCorrelation between D and R (after transform and or subset):\n")
  print(cor.test(D,R))
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # is the y-variable an ordered factor?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  if ((is.factor(y[1]))&(ordered.fac.treatment[1] == c("as.num"))){
    cat(yname, " is an ordered factor. options are to treat as numeric or in a binomial fashion.\n")    
    cat(yname, " will be converted to numeric \n")
    y <- as.numeric(y)}
  
  if ((is.factor(y[1]))&(ordered.fac.treatment[1] == c("as.binomial"))){
    cat("\n",yname, " is an ordered factor. options are to treat as numeric or in a binomial fashion !!!\n")
    cat(yname, " is an ordered factor and will be converted to a count table !!!\n")
    y <- ordered.var.2.sim.count.table(y)}
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # apply centering and standardization ####
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  CGM <- function(x){(x-mean(x,na.rm=T))}#/sd(x)}
  CWC <- function(x,LOC){mns.x <- tapply(x,LOC,mean,na.rm=T)
                         mns.x <- mns.x[match(LOC,names(mns.x))]
                         mns.x <- x - mns.x
                         return(mns.x)}  
  if (centering == 'CGM'){D <- CGM(D) ; R <- CGM(R)
                          cat('\nPredictor variables were centered at the grand mean (CGM)\n') }
  if (centering == 'CWC'){
    D <- CWC(x=D,LOC=LOC.)
    R <- CWC(x=R,LOC=LOC.)
    
    cat('\nPredictor variables were centered within clusters (CWC)\n') 
  }
  
  if(standardisation){
    cat('\nPredictor variables were standardised using grand sd\n')
    D <- D/sd(D,na.rm=T)
    R <- R/sd(R,na.rm=T)
  }
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # plot the predictor variables
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
  win.graph(w=9,h=4);par(mfrow=c(1,3))
  hist(unique(D),col="light grey");hist(unique(R),col="light grey");plot(R,D,pch=16)
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # calibrate the main model and it's sub-models
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
  if ((ordered.fac.treatment[1] == c("as.binomial"))& (!is.null(dim(y)))){
    #for models with ordinal response variables treated as binomial counts  
    mod.DR <- glmer(formula = y ~ D*R + (1|LOC.),weights=nbirds.,family="binomial")
    mod.D_R <- glmer(formula = y ~ D + R + (1|LOC.),weights=nbirds.,family="binomial")
    mod.D <- glmer(formula = y ~ D + (1|LOC.),weights=nbirds.,family="binomial")
    mod.R <- glmer(formula = y ~ R + (1|LOC.),weights=nbirds.,family="binomial")
    mod.0 <- glmer(formula = y ~ 1 + (1|LOC.),weights=nbirds.,family="binomial")
    }
  else{
    #for all other models
    mod.DR <- lmer(formula = y ~ D*R + (1|LOC.),weights=nbirds.)
    mod.D_R <- lmer(formula = y ~ D + R + (1|LOC.),weights=nbirds.)
    mod.D <- lmer(formula = y ~ D + (1|LOC.),weights=nbirds.)
    mod.R <- lmer(formula = y ~ R + (1|LOC.),weights=nbirds.)
    mod.0 <- lmer(formula = y ~ 1 + (1|LOC.),weights=nbirds.)
  }  
  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  #in the D*R model, under what conditions are the D (R) terms negative?
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  cat("Conditional effects of predictors in full model, i.e. including R*D interaction\n")
  win.graph(w=8,h=8);par(mfcol=c(2,2))
  
  r1 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,plott=F)
  r2 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,switch.vars = T,plott=F)
  
  r12 <- c(-1,1)*max(abs(c(r1,r2)),na.rm=T)
    r1 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,ylim.=r12)
  r2 <- Plot_conditional_effects_interaction_models(mymod=mod.DR,ylim.=r12,switch.vars = T)
    
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  # get the AIC  
  #++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  #since 'n' is not defined in a mixed model, calculate a liberal AIC (mod.AIC.CON)
  #assuing n is the nr of data points (nobs.HI), and a liberal AIC (mod.AIC.LIB) assuming n is the nr of groupd (nobs.LO)
  nobs.HI <- nrow(mod.DR@frame)
  nobs.LO <- nlevels(mod.DR@frame$LOC)
  
  mod.AIC.LIB <- ICtab(mod.DR,mod.D,mod.R,mod.0,mod.D_R,sort=F,type="AICc",nobs=nobs.HI)
  mod.AIC.CON <- ICtab(mod.DR,mod.D,mod.R,mod.0,mod.D_R,sort=F,type="AICc",nobs=nobs.LO)
  
  mod.AIC <- cbind.data.frame(round(mod.AIC.LIB$dAICc,1),round(mod.AIC.CON$dAICc,1),mod.AIC.LIB$df)
  colnames(mod.AIC) <- c('dAICc.lib','dAICCc.con','df')
  rownames(mod.AIC) <- c('mod.DR','mod.D','mod.R','mod.0','mod.D_R')
  
  #calculate Akaike weights sensu http://www.planta.cn/forum/files_planta/2006_nobbs_ecology_153.pdf
  Akaike.weight <- function(AIC.vec){AIC.vec <- AIC.vec-min(AIC.vec,na.rm=T)
                                     Akai.weight.denom <- sum(exp(-2*AIC.vec))
                                     Akai.weight <- (exp(-2*AIC.vec))/Akai.weight.denom
                                     return(Akai.weight)
  }
  A.weight <- round(Akaike.weight(mod.AIC$dAICc.lib),2)
  mod.AIC <- cbind.data.frame(mod.AIC,A.weight)
{'fraction_non_alphanumeric': 0.2522465800469349, 'fraction_numerical': 0.009844885810772136, 'mean_word_length': 3.3625468164794006, 'pattern_counts': {'":': 0, '<': 118, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3354049', 'n_tokens_mistral': 6440, 'n_tokens_neox': 5786, 'n_words': 1323}