library(lme4) library(lmerTest) library(optimx) library(effects) library(Hmisc) library(cAIC4) library(car) library(MuMIn) library(MASS) library(party) library(gridExtra) library(latticeExtra) setwd(".../DATA") {dat=read.csv("rsf_table.csv", sep=";") dat$pb = factor(dat$pb) dat$anid = factor(dat$anid) } #======================================================================================= # Correlation analysis, Spearman #======================================================================================= names(dat) with(dat,cor.test(siveg,nveg,method="spearman")) # Check for collinearity between independent variables, Spearman's Rho>0.7 mat.data<-as.matrix(dat) cor(dat[,-c(1:2)], method="spearman") rcorr(mat.data[,-c(1:2)], type="spearman") cor.mat=round(rcorr(mat.data[,-c(1:2)], type="spearman")$r,2) cor.mat[abs(cor.mat)<.5]=paste(".") noquote(cor.mat) #Calculate correlations for a selection of variables t.data<-with(dat,data.frame(suf,g,mun,hc,r,a,u,cvbiot,nbiot,sibiot,nveg,siveg)) mat.data<-as.matrix(t.data) rcorr(mat.data[,1:length(t.data)], type="spearman") #======================================================================================= # Visualization of correlations #======================================================================================= #correlation tree round(cor(dat[,-c(1:2,14)],use="complete.obs", method="spearman"), 3) #var <- c(3:ncol(dat)) var = c(3:9, 11, 13, 15, 16, 22) v <- as.formula(paste("~", names(dat[,var]), collapse="+")) plot(varclus(v, data =dat[,var])) # yielding the figure below abline(h=.3) # PCA Biplot bi.dat=na.omit(dat) bi <- princomp(bi.dat[,-c(1:2,14)], cor = TRUE) biplot(bi) #============================================= # univariate GLM - stepwise model building #============================================= names(dat) # Positions of continuous Variables n=ncol(dat) varlist = names(dat)[3:n] response <- dat$pb datsq = subset(dat, select=-c(dw, dpuo, dbuh, dmun, dd, dg, da)) varsq = names(datsq)[3:ncol(datsq)] datsq = dat # null M0 <- glmer(response ~ 1 + (1|anid), na.action=na.omit, family=binomial, data=dat) summary(M0) # loop for simple univariate model building models <- lapply(varlist, function(x) { glmer(substitute(response ~ i + (1|anid), list(i = as.name(x))), data = dat, family=binomial) # binomial(link="logit") }) # loop for squared term univariate model building [x + I(x^2)] sqmodels <- lapply(varsq, function(x) { glmer(substitute(response ~ i + I(i^2) + (1|anid), list(i = as.name(x))), data = datsq, family=binomial) }) ## Output of model summaries lapply(models, summary) lapply(sqmodels, summary) lapply(models, AIC) lapply(sqmodels, AIC) for(i in 1:length(models)){ assign(paste0("M", i), models[[i]]) if(i<=length(varsq)) assign(paste0("Msq", i), sqmodels[[i]]) } anova(M0, M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, M17, M18, M19, M20, M21) anova(Msq1, Msq2, Msq3, Msq4, Msq5, Msq6, Msq7, Msq8) anova(M11, M12) anova(M13, M14, M15, M16) #============================================= # multivariate GLM - model building #============================================= # transform distance units m to km { dat3 = subset(dat, select=-c(d, a, nbiot, sibiot, cvbiot, nveg, cvveg, dbuh, dmun, da)) dat3$dw = dat3$dw/1000 dat3$dpuo = dat3$dpuo/1000 dat3$dd = dat3$dd/1000 dat3$dg = dat3$dg/1000} # full list of variables varlist = c('dw', 'buh', 'mun', 'd', 'hc', 'r', 'g', 'a', 'dpuo', 'DG', 'nveg', 'siveg', 'dbuh', 'dmun', 'dd', 'dg', 'da') varsq = c('buh', 'mun', 'd', 'hc', 'r', 'g', 'a', 'DG', 'nveg', 'siveg') # stepwise reduced list varlist = c('hc') varsq = varlist # loop for simple univariate model building models <- lapply(varlist, function(x) { glmer(substitute(response ~ dpuo + I(dpuo^2) + siveg + I(siveg^2) + dw + I(dw^2) + dd + I(dd^2) + r + I(r^2) + buh + I(buh^2) + DG + dg + I(dg^2) + mun + I(mun^2) + i + (1|anid), list(i = as.name(x))), data = dat3, family=binomial) }) # loop for squared term univariate model building [x + I(x^2)] sqmodels <- lapply(varsq, function(x) { glmer(substitute(response ~ dpuo + I(dpuo^2) + siveg + I(siveg^2) + dw + I(dw^2) + dd + I(dd^2) + r + I(r^2) + buh + I(buh^2) + DG + dg + I(dg^2) + mun + I(mun^2) + i + I(i^2) + (1|anid), list(i = as.name(x))), data = dat3, family=binomial) }) # Output summaries for model selection data.frame(i="_", frml=unlist(lapply(models,function(x){as.character(formula(x)[[3]])[2]})), i="_", M=paste0("M",1:length(models)), i="_", npar=unlist(lapply(models,function(x){ncol(x@frame)})), i="_", AIC=unlist(lapply(models, AIC))) data.frame(i="_", frml=unlist(lapply(sqmodels,function(x){as.character(formula(x)[[3]])[2]})), i="_", M=paste0("Msq",1:length(sqmodels)), i="_", npar=unlist(lapply(sqmodels,function(x){ncol(x@frame)})), i="_", AIC=unlist(lapply(sqmodels, AIC))) #============================================= # full additive multivariate Model #============================================= # Best stepwise built multivariate Model (Distances in km) Mikm = glmer(pb ~ dpuo + I(dpuo^2) + siveg + I(siveg^2) + dw + I(dw^2) + dd + I(dd^2) + r + I(r^2) + buh + I(buh^2) + DG + dg + I(dg^2) + mun + I(mun^2) + hc + (1|anid), data = dat3, family=binomial) summary(Mikm) # Effect plots vl = sub('I', NA, names(coefficients(Mikm)$anid)[-1]) vl = vl[!is.na(vl)] for(i in 1:length(vl)){ assign(paste0("Mikm.", vl[i]), plot(effect(vl[i], Mikm), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))))) if(i==1) Mikmplot = paste0("Mikm.", sub(':', '.', vl[i])) else Mikmplot = c(Mikmplot, paste0("Mikm.", sub(':', '.', vl[i]))) }; noquote(paste0(Mikmplot,",")) grid.arrange(Mikm.dpuo, Mikm.dw, Mikm.dg, Mikm.hc, Mikm.r, Mikm.mun, Mikm.buh, Mikm.dd, Mikm.DG, Mikm.siveg, nrow=3, ncol=4) #============================================= # add Interactions #============================================= # Mi4a Mi4a = glmer(pb ~ dpuo + dw + dg + hc + r + I(r^2) + mun + buh + I(buh^2) + dd + I(dd^2) + DG + siveg + I(siveg^2) + hc : (mun + buh + dd + I(dd^2) + DG) + (r + I(r^2)) : (buh + I(buh^2) + DG + siveg + I(siveg^2)) + (1|anid), data = dat3, family=binomial) summary(Mi4a) # AIC=8450.5 # unrealistic variances (variance inflation), drop squared r interaction # Mi4b: Mikm +(hc+r):(mun12+buh12+dd12+DG12+siveg12) -r:siveg12 -hc:siveg2 Mi4b = glmer(pb ~ dpuo + I(dpuo^2) + dw + I(dw^2) + dg + I(dg^2) + hc + r + I(r^2) + mun + I(mun^2) + buh + I(buh^2) + dd + I(dd^2) + DG + siveg + I(siveg^2) + hc : (mun + I(mun^2) + buh + I(buh^2) + dd + I(dd^2) + DG + siveg) + r : (mun + I(mun^2) + buh + I(buh^2) + dd + I(dd^2) + DG) + (1|anid), data = dat3, family=binomial) summary(Mi4b) # AIC=8247.2 # further variable interactions dropped because of extreme variance inflation # Mi4c: Mi4b -hc:dd2 -hc:buh2 -r:mun2 -r:dd2 -r:buh12 -hc:mun2 -siveg2 Mi4c = glmer(pb ~ dpuo + I(dpuo^2) + dw + I(dw^2) + dg + I(dg^2) + hc + r + I(r^2) + mun + I(mun^2) + buh + I(buh^2) + dd + I(dd^2) + DG + siveg + hc : (mun + buh + dd + DG + siveg) + r : (mun + dd + DG) + (1|anid), data = dat3, family=binomial) summary(Mi4c) # AIC=8265.7 ############## best model accuracy without variance inflation summary(Mi4c)$coef #Mi4c: diversity indices dropped to avoid NA Mi4cohneNA = glmer(pb ~ dpuo + I(dpuo^2) + dw + I(dw^2) + dg + I(dg^2) + hc + r + I(r^2) + mun + I(mun^2) + buh + I(buh^2) + dd + I(dd^2) + hc : (mun + buh + dd) + r : (mun + dd) + (1|anid), data = dat3, family=binomial) summary(Mi4cohneNA) summary(Mi4cohneNA)$coef #============================================= # Effect plots #============================================= ### vereinfachter Interaktionsplot Minput = Mi4c ## enter model to plot vl = sub('I', NA, names(coefficients(Minput)$anid)[-1]) vl = vl[!is.na(vl)] for(i in 1:length(vl)){ assign(paste0("Mi.", sub(':', '.', vl[i])), plot(effect(vl[i], Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))))) if(i==1) Miplot = paste0("Mi.", sub(':', '.', vl[i])) else Miplot = c(Miplot, paste0("Mi.", sub(':', '.', vl[i]))) } noquote(paste0(Miplot,",")) png("../Mi4cNA1.png", width=1200, height=800, units="px", pointsize=20) grid.arrange(Mi.dpuo, Mi.dw, Mi.dg, Mi.hc, Mi.r, Mi.mun, Mi.buh, Mi.dd, Mi.DG, Mi.siveg, nrow=3, ncol=4) dev.off() png("../Mi4cNA2.png", width=1200, height=800, units="px", pointsize=20) grid.arrange(Mi.hc.mun, Mi.hc.buh, Mi.hc.dd, Mi.hc.DG, Mi.hc.siveg, nrow=2, ncol=3) dev.off() png("../Mi4cNA3.png", width=1200, height=800, units="px", pointsize=20) grid.arrange(Mi.r.mun, Mi.r.dd, Mi.r.DG, nrow=2, ncol=3) dev.off() # result plot 1) Main Effects png("../Mi4c_maineffects.png", width=800, height=1000, units="px", pointsize=20) lattice.options( layout.heights=list(bottom.padding=list(x=-.8), top.padding=list(x=-.1)), layout.widths=list(left.padding=list(x=0), right.padding=list(x=-1))) grid.arrange(plot(effect("dpuo", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="Distance variables", xlab="Distance to infrastructure [km]"), plot(effect("hc", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="Availability variables", xlab="Sand heaths"), plot(effect("DG", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="Diversity indices", xlab="Cover rate of food plants"), plot(effect("dw", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Distance to forest [km]"), plot(effect("r", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Natural grasslands"), plot(effect("siveg", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Shannon Index of plant species"), plot(effect("dg", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Distance to pastures [km]"), plot(effect("mun", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Raised bogs, mires, fens"), plot(effect("dd", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Distance to dunes [km]"), plot(effect("buh", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1))), main="", xlab="Loose shrub formations"), nrow=4, ncol=3, layout_matrix=rbind(c(1, 2, 3), c(4, 5, 6), c(7, 8, NA), c(9, 10, NA))) dev.off() # result plot 2) Interaction Effects png("../Mi4c_interactioneffects1.png", width=800, height=1000, units="px", pointsize=20) lattice.options( layout.heights = list(bottom.padding=list(x=-.8), top.padding=list(x=-.1)), layout.widths = list(left.padding=list(x=0), right.padding=list(x=-1))) trellis.par.set(list(axis.components = list(top = list(tck=0)))) plist = list(plot(effect("hc:mun", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="Sand heath interactions", xlab=""), plot(effect("hc:buh", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="", xlab=""), plot(effect("hc:dd", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="", xlab=""), plot(effect("hc:DG", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="", xlab=""), plot(effect("hc:siveg", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="", xlab="Sand heath availability"), plot(effect("r:mun", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="Grassland interactions", xlab=""), plot(effect("r:dd", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="", xlab=""), plot(effect("r:DG", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1)), x=list(cex=.9)), main="", xlab="Grassland availability")) plist = lapply(plist, function(x){update(x, scales=list(alternating=F, x=list(rot=90, at=seq(0,1,.25), labels=seq(0,1,.25))))}) grid.arrange(grobs=plist, nrow=5, ncol=2, layout_matrix=cbind(c(1, 2, 3, 4, 5), c(6, NA, 7, 8, NA))) dev.off() str(trellis.par.get("axis.components")) mytheme <- list(axis.components = list(left = list(tck=1), right = list(tck=1))) #Bsp stackoverflow tck=4 mytheme <- list(axis.components = list(scales = list(alternating = F, x=list(rot=45, at=seq(0,1,.3), labels=seq(0,1,.3)))), axis.text=list(alpha=1, cex=1)) #, x=list(rot=45), alternating = F, fontsize=list(text=12), par.xlab.text=list(cex=1) mytheme <- list(axis.components = list(top = list(tck=0))) trellis.par.set(mytheme) plot(effect("dpuo", Minput), type="response", rescale.axis=F, axes=list(y=list(lim=c(0,1), cex=.8), x=list(cex=.8)), main="Distance variables", xlab="Distance to infrastructure [km]") ar=plot(effect("hc:mun", Minput), type="response", rescale.axis=F, layout=c(5,1), axes=list(y=list(lim=c(0,1), cex=1), x=list(cex=.9)), main="Sand heath interactions", xlab="") update(ar, scales=list(x=list(at=seq(0,1,.3), labels=seq(0,1,.3)))) update(ar, scales=list(alternating=F, x=list(rot=90, at=seq(0,1,.25), labels=seq(0,1,.25)))) update(ar, scales=list(alternating=F, x=list(rot=90)))