############################################################################################################# # # TITLE of publication # # Radinger J. and Wolter C. # # R-script for analysis of the single and joint multiple effects of habitat, dispersal and barriers # on the presence/absence of 17 fish species. Method used: GLMM # ############################################################################################################# #library("devtools"); install_github("lme4",user="lme4",ref="release") #for developement version of lme4 #library("devtools"); install_github("lme4",user="lme4") #Import libraries library(reshape2) library(lme4) library(scales) library(plyr) # Resultsdir results_dir <- "path/to/habitat_dispersal_barriers_raw_data.csv" # load df habitat_dispersal_barriers_raw_data <- read.csv(paste(results_dir,"habitat_dispersal_barriers_raw_data.csv",sep="")) head(habitat_dispersal_barriers_raw_data) ####################### Prepare dataframe ########################## habitat_dispersal_barriers_raw_data$SiteID <- as.factor(habitat_dispersal_barriers_raw_data$SiteID) habitat_dispersal_barriers_raw_data$variable <- gsub("point_fidimo_corrected","dispersal",habitat_dispersal_barriers_raw_data$variable) habitat_dispersal_barriers_raw_data$variable <- gsub("point_barrier_effect_corrected","barriers",habitat_dispersal_barriers_raw_data$variable) habitat_dispersal_barriers_raw_data$variable <- gsub("point_MaxEnt","habitat",habitat_dispersal_barriers_raw_data$variable) habitat_dispersal_barriers_raw_data$year <- unlist(lapply(strsplit(as.character(habitat_dispersal_barriers_raw_data$variable), "_"),"[",2)) habitat_dispersal_barriers_raw_data$variable <- unlist(lapply(strsplit(as.character(habitat_dispersal_barriers_raw_data$variable), "_"),"[",1)) habitat_dispersal_barriers_raw_data$variable <- factor(habitat_dispersal_barriers_raw_data$variable, levels=c("habitat","dispersal","barriers")) df_habitat_dispersal_barriers <- dcast(habitat_dispersal_barriers_raw_data, species+SiteID+presabs~variable,mean) df_habitat_dispersal_barriers$habitat_centered <- scale(df_habitat_dispersal_barriers$habitat,scale=FALSE,center=median(df_habitat_dispersal_barriers$habitat)) df_habitat_dispersal_barriers$dispersal_centered <- scale(df_habitat_dispersal_barriers$dispersal,scale=FALSE,center=median(df_habitat_dispersal_barriers$dispersal)) df_habitat_dispersal_barriers$barriers_centered <- scale(df_habitat_dispersal_barriers$barriers,scale=FALSE,center=median(df_habitat_dispersal_barriers$barriers)) ################################################ Descriptive comparisons ############################################ # Test for Normality shapiro.test(df_habitat_dispersal_barriers$dispersal[df_habitat_dispersal_barriers$presabs==0]) summary(factor(df_habitat_dispersal_barriers$presabs)) # Habitat suitability # H0: Absence records show lower habitat suitability wtest_habitat <- wilcox.test(df_habitat_dispersal_barriers$habitat[df_habitat_dispersal_barriers$presabs==0], df_habitat_dispersal_barriers$habitat[df_habitat_dispersal_barriers$presabs==1], paired=FALSE, alternative="l") # Z-statistic qnorm(wtest_habitat$p.value) # Species dispersal # H0: Absence records show lower dispersal probability wtest_dispersal <- wilcox.test(df_habitat_dispersal_barriers$dispersal[df_habitat_dispersal_barriers$presabs==0], df_habitat_dispersal_barriers$dispersal[df_habitat_dispersal_barriers$presabs==1], paired=FALSE, alternative="l") # Z-statistic qnorm(wtest_dispersal$p.value) # Barrier effect # H0: Absence records show greater effect of barriers (lower (more negative) barrier values) wtest_barriers <- wilcox.test(df_habitat_dispersal_barriers$barriers[df_habitat_dispersal_barriers$presabs==0], df_habitat_dispersal_barriers$barriers[df_habitat_dispersal_barriers$presabs==1], paired=FALSE, alternative="l") # Z-statistic qnorm(wtest_barriers$p.value) ############################################################################################# ####################################### year-mean GLMM Models ################################ #Goodness of fit (Pseudo R-squared) / correlation between: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2008q2/000713.html # http://glmm.wikidot.com/faq r2.corr.mer <- function(m) { lmfit <- lm(model.response(model.frame(m)) ~ fitted(m)) summary(lmfit)$r.squared } #### Single GLMMs #### mod_s1 <- glmer(presabs~habitat_centered+(habitat_centered|species), family=binomial,data=df_habitat_dispersal_barriers) ci_mod_s1 <- confint.merMod(mod_s1,method="boot") mod_s2 <- glmer(presabs~dispersal_centered+(dispersal_centered|species), family=binomial,data=df_habitat_dispersal_barriers) ci_mod_s2 <- confint.merMod(mod_s2,method="boot") mod_s3 <- glmer(presabs~barriers_centered+(barriers_centered|species), family=binomial,data=df_habitat_dispersal_barriers) ci_mod_s3 <- confint.merMod(mod_s3,method="boot") mod_s4 <- glmer(presabs~habitat_centered+(1|species), family=binomial,data=df_habitat_dispersal_barriers) ci_mod_s4 <- confint.merMod(mod_s4,method="boot") mod_s5 <- glmer(presabs~dispersal_centered+(1|species), family=binomial,data=df_habitat_dispersal_barriers) ci_mod_s5 <- confint.merMod(mod_s5,method="boot") mod_s6 <- glmer(presabs~barriers_centered+(1|species), family=binomial,data=df_habitat_dispersal_barriers) ci_mod_s6 <- confint.merMod(mod_s6,method="boot") # Multiple GLMMs #Random slopes and random intercept for all species and fixed additive and interaction effects (beyond optimal model) mod_m0 <- glmer(presabs~1+(1|species), family=binomial,data=df_habitat_dispersal_barriers) mod_m1 <- glmer(presabs~barriers_centered+dispersal_centered+habitat_centered+barriers_centered*habitat_centered+dispersal_centered*habitat_centered+dispersal_centered*barriers_centered+(barriers_centered+habitat_centered+dispersal_centered|species), family=binomial,data=df_habitat_dispersal_barriers) (deviance(mod_m0)-deviance(mod_m1))/deviance(mod_m0) ci_mod_m1 <- confint.merMod(mod_m1,method="boot") mod_m2 <- glmer(presabs~barriers_centered+dispersal_centered+habitat_centered+(barriers_centered+habitat_centered+dispersal_centered|species), family=binomial,data=df_habitat_dispersal_barriers) (deviance(mod_m0)-deviance(mod_m2))/deviance(mod_m0) ci_mod_m2 <- confint.merMod(mod_m2,method="boot") mod_m3 <- glmer(presabs~barriers_centered+dispersal_centered+habitat_centered+barriers_centered*habitat_centered+dispersal_centered*habitat_centered+dispersal_centered*barriers_centered+(1|species), family=binomial,data=df_habitat_dispersal_barriers) (deviance(mod_m0)-deviance(mod_m3))/deviance(mod_m0) ci_mod_m3 <- confint.merMod(mod_m3,method="boot") mod_m4 <- glmer(presabs~barriers_centered+dispersal_centered+habitat_centered+(1|species), family=binomial,data=df_habitat_dispersal_barriers) (deviance(mod_m0)-deviance(mod_m4))/deviance(mod_m0) ci_mod_m4 <- confint.merMod(mod_m4,method="boot") ############# Correlation of Random slopes of dispersal ################### ranef_df <- merge((ranef(mod_m2)$species),(plyr::ddply(habitat_dispersal_barriers_raw_data[,c("species","L","AR")],~species,summarize,L=mean(L),AR=mean(AR))),by.x="row.names",by.y="species") cor.test(ranef_df$dispersal,ranef_df$L,method="spearman") cor.test(ranef_df$habitat,ranef_df$L,method="spearman") ##################################################################################################### ############################### Development over years - model ##################################### # load df again and cast habitat_dispersal_barriers_raw_data <- read.csv(paste(results_dir,"habitat_dispersal_barriers_raw_data.csv",sep="")) df_habitat_dispersal_barriers_years <- dcast(habitat_dispersal_barriers_raw_data, species+SiteID+presabs~variable,mean) # Function to extract parameter estimates coef_lmer_comp <- function(.){ #Bootstrapped CIs bootmod <- confint.merMod(.,method="boot") data.frame(Estimates=c(unname(sigma(.) * getME(., "theta")),getME(., "beta")), BootCI_lwr=bootmod[,1], BootCI_upr=bootmod[,2], AIC=AIC(.), deviance=deviance(.)) } y <- paste(c(1:9),"y",sep="") # Scaling factor for dispersal and barrier effect (see Schielzeth 2010, Simple means to improve the interpretability of regression coefficients. Methods in Ecology and Evolution. 17, 3428-3447) scale_fidimo <- sd(c(t(df_habitat_dispersal_barriers_years[,paste("point_fidimo_corrected_",y,sep="")]))) scale_barrier <- sd(c(t(df_habitat_dispersal_barriers_years[,paste("point_barrier_effect_corrected_",y,sep="")]))) scale_habitat <- sd(rep(as.vector(df_habitat_dispersal_barriers_years[,"point_MaxEnt"]),length(y))) center_fidimo <- median(c(t(df_habitat_dispersal_barriers_years[,paste("point_fidimo_corrected_",y,sep="")]))) center_barrier <- median(c(t(df_habitat_dispersal_barriers_years[,paste("point_barrier_effect_corrected_",y,sep="")]))) center_habitat <- median(rep(as.vector(df_habitat_dispersal_barriers_years[,"point_MaxEnt"]),length(y))) mod_results_complete <- list() # Following loop may take very long time due to bootstrapping for(i in y){ print(paste("This is run ", i, sep="")) HS <- as.vector(eval(parse(text=paste("scale(df_habitat_dispersal_barriers_years$point_MaxEnt,scale=scale_habitat,center=center_habitat)",sep="")))) DI <- as.vector(eval(parse(text=paste("scale(df_habitat_dispersal_barriers_years$point_fidimo_corrected_",i,",scale=scale_fidimo,center=center_fidimo)",sep="")))) BE <- as.vector(eval(parse(text=paste("scale(df_habitat_dispersal_barriers_years$point_barrier_effect_corrected_",i,",scale=scale_barrier,center=center_barrier)",sep="")))) df_mod <- data.frame(presabs=df_habitat_dispersal_barriers_years$presabs,HS,DI,BE,species=df_habitat_dispersal_barriers_years$species) # GLMM model mod_comp <- glmer(presabs~HS+DI+BE+(1|species),family=binomial,data=df_mod) # Collect Model results (inkl. bootstrapped CI) mod_results_complete[[i]] <- cbind(coef_lmer_comp(mod_comp),Time=as.numeric(substring(i,1,1))) } saveRDS(mod_results_complete,paste(results_dir,"glmm_over_years.rds",sep=""))