Last active
November 18, 2015 23:48
-
-
Save mjwestgate/a332ac264a17c60de9be to your computer and use it in GitHub Desktop.
Code for analysing the ACT Frogwatch dataset (2002-2014)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# This script analyses data from the ACT Frogwatch program, as described by Westgate et al. (2015) | |
# 'Citizen science program shows urban areas have lower occurrence of frog species, but not accelerated declines' | |
# PLOS ONE doi: 10.1371/journal.pone.0140973 | |
# http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0140973 | |
# Data for these analyses is available at Dryad: | |
# https://datadryad.org/resource/doi:10.5061/dryad.75s51 | |
# Below I assume that relevant datasheets are available as .csv files | |
# Note: Section 2 onwards requires functions written specifically for this analysis | |
# see 'Frogwatch_2015_functions.R' | |
# call necessary libraries and defaults | |
library(boot) # logit and inv.logit | |
library(lme4) # glmms | |
library(spdep) # autocovariates | |
library(AICcmodavg) # model selection by AICc | |
## 1. IMPORT AND PROCESS RAW DATA | |
# get csv data | |
data<-read.csv("observations.csv", stringsAsFactors=FALSE) | |
sites<-read.csv("sites.csv", stringsAsFactors=FALSE) | |
# sort out covariates, data structure etc. | |
data$veg.local<-apply(data[, c("veg.aqu", "veg.bank")], 1, sum) | |
data$pc.canopy<-logit(data$pc.canopy) | |
data$pc.urban<-logit(data$pc.urban) | |
# extract and organize species data | |
species<-data[, c(11:18)] | |
species.names<-colnames(species) | |
n.species<-length(species.names) | |
species.names.pretty<-c("C. parinsignifera", "C. signifera", | |
"Lim. dumerilii", "Lim. peronii", "Lim. tasmaniensis", | |
"L. peronii", "L. verreauxii", "U. laevigata") | |
species$richness<-apply(species, 1, sum) | |
# get subset of predictor variables | |
predictors<-data[, c(1, 3, 4, 2, 5, 19, 9, 10)] | |
# Use 'sites' dataset to calculate autocorrelation for each species across all years | |
# Note: this code altered from Dormann et al. 2007 Ecography 30: 609-628 | |
# get neighbour lists etc. | |
coords<-as.matrix(sites[, c("latitude", "longitude")]) | |
nb.list<-dnearneigh(coords, 0, 20000) # units are in metres, ergo kernel is 20 km | |
nb.weights<-nb2listw(nb.list) | |
# calculate autocovariate for each response variable | |
ac.fun<-function(sp.name){autocov_dist(sites[, sp.name], coords, nbs=20000, type="inverse")} | |
ac.data<-sapply(species.names, FUN= ac.fun) | |
ac.data<-as.data.frame(ac.data) | |
colnames(ac.data)<-paste("ac", colnames(ac.data), sep=".") | |
ac.data$ac.richness<-autocov_dist(sites$richness, coords, nbs=20000, type="inverse") | |
# merge with observation data | |
ac.data$site.code<-sites$site.code | |
predictors<-merge(predictors, ac.data, by="site.code") | |
# format 'predictors' for analysis | |
for(i in 1:3){predictors[, i]<-as.factor(predictors[, i])} | |
for(i in 4:ncol(predictors)){predictors[, i]<-scale(predictors[, i])} | |
# get a single dataset for analysis | |
data.final<-cbind(predictors, species) | |
# 2. MODEL USING GLMMS | |
# set out model formulae for each species | |
basic.model<-"spp ~ veg.local + log.size + water.type + ac" | |
# the above variables are included in all models. This is because we want to control for their effects. | |
# Only one of them (veg.local) is of interest to us, in terms of its' potential additive or interactive effects with urbanisation. | |
model.list<-c( | |
simple= "", | |
canopy=" + pc.canopy", | |
urban=" + pc.urban", | |
year=" + year", | |
urban.canopy=" + pc.urban + pc.canopy", | |
urban.year=" + pc.urban + year", | |
year.canopy=" + year + pc.canopy", | |
all.additive=" + pc.urban + year + pc.canopy", | |
urbanXyear=" + pc.urban * year", | |
urbanXveg=" + pc.urban + pc.urban:veg.local", | |
urbanXcanopy=" + pc.urban * pc.canopy", | |
urbanXyearXveg=" + pc.urban * year + pc.urban:veg.local + veg.local:year + pc.urban:year:veg.local", | |
urbanXyearXcanopy=" + pc.urban * year * pc.canopy", | |
full.model="+ pc.urban * year * veg.local * pc.canopy") | |
model.list.final<-paste(basic.model, model.list, " + (1 | site.code)", sep="") | |
names(model.list.final)<-names(model.list) | |
model.list.final<-as.list(model.list.final) | |
# calculate individual species models | |
# NOTE: This can take a long time to run - be patient! | |
final.models<-sapply(species.names, FUN=function(x){ | |
urban.models(x, source=data.final, model.list= model.list.final)}) | |
names(final.models)<-paste(rep(species.names, each=2), | |
rep(c("aictab", "model"), length(species.names)), sep=".") | |
# add species richness models | |
richness.models<-urban.models("richness", | |
source=data.final, model.list=model.list.final, richness=TRUE) | |
names(richness.models)<-c("richness.aictab", "richness.model") | |
final.models<-append(final.models, richness.models) | |
# get model selection results | |
n.models<-length(model.list.final) | |
aictab.final<-as.data.frame(do.call("rbind", final.models[c(1:n.models)*2-1])) | |
aictab.final$species<-rep(c(species.names, "richness"), each= n.models) | |
aictab.final<-aictab.final[, c(ncol(aictab.final), 1:(ncol(aictab.final)-1))] | |
# write.csv(aictab.final, "model_selection_results.csv", row.names=FALSE) # export | |
# get coefficients for 'final' models | |
data.names<-paste(c(species.names, "richness"), "model", sep=".") | |
coef.list<-lapply(data.names, | |
FUN=function(x){summary(final.models[[x]])$coefficients}) | |
names(coef.list)<-c(species.names, "richness") | |
coef.dframe<-as.data.frame(do.call("rbind", coef.list)) | |
coef.dframe$variable<-rownames(coef.dframe) | |
name.vector<-unlist(lapply(names(coef.list), FUN=function(x){rep(x, nrow(coef.list[[x]]))})) | |
coef.dframe$response<-name.vector | |
rownames(coef.dframe)<-NULL | |
coef.dframe<-coef.dframe[, c(6, 5, 1:4)] | |
# write.csv(coef.dframe, "coefficients.csv", row.names=FALSE) # export | |
# 3. PLOT MODEL RESULTS | |
# basic specifications | |
# colors | |
light<-rgb(0.85, 0.6, 0.3, 1, maxColorValue=1) | |
dark<-rgb(0.1, 0.05, 0.6, 1, maxColorValue=1) | |
# draw static effects | |
pdf("Fig2.pdf", width=7.5, height=3.5) | |
plot.static.effects(spp=c(species.names, "richness"), model.list=final.models, | |
fx.names=c("(Intercept)", "pc.urban", "log.size", "water.typestill", "veg.local", "pc.canopy", "year"), | |
fx.labels=c("Intercept", "% Urban", "Waterbody size", "Type = still", "Vegetation", "Year", "% Canopy"), | |
label.offsets=c(0.3, 0.3, 0.55, 0.35, 0.32, 0.17, 0.32), #rep(0.2, 7), | |
stagger=c(FALSE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE), | |
cols=c(light, dark), | |
x.limits=data.frame( | |
#min=c(-5, -0.7, -0.6, -2.2, -0.5, -0.3, -0.7), | |
#max=c(1, 0.3, 1.4, 3.5, 1.5, 0.6, 0.6), | |
axis.min=c(-4, -1, -0.5, -2, -0.5, -0.25, -1), | |
axis.max=c(2, 0.25, 1.5, 2, 1.5, 0.5, 0.5), | |
axis.step=c(2, 0.25, 0.5, 2, 0.5, 0.25, 0.25)), | |
spp.labels=c(species.names.pretty, "Species richness")) | |
dev.off() | |
# draw predictons plot | |
# Note: panel labels render ok in PDF on my machine, but don't plot correctly in the window. | |
screen.matrix<-matrix(data=c( | |
0.03, 0.25, 0.5, 1, | |
0.25, 0.47, 0.5, 1, | |
0.47, 0.69, 0.5, 1, | |
0.71, 0.86, 0.5, 1, | |
0.85, 1.0, 0.5, 1, | |
0.03, 0.25, 0.01, 0.51, | |
0.25, 0.47, 0.01, 0.51, | |
0.47, 0.69, 0.01, 0.51, | |
0.71, 0.86, 0.01, 0.51, | |
0.85, 1.0, 0.01, 0.51, | |
0, 1, 0, 1), | |
nrow=11, ncol=4, byrow=TRUE) | |
pdf("Fig3.pdf", width=7.5, height=5.5) | |
par(cex=0.52, mar=c(3, 2, 2, 1)) | |
split.screen(screen.matrix) | |
# panel d | |
plot.full.model(final.models$cpar.model, screens=c(4, 5), | |
cols=c(rep(light, 2), rep(dark, 2)), | |
ltys=c(2, 1, 2, 1), | |
link="logit", | |
low.key=TRUE, | |
labels=list(group1=c("Low", "High"), | |
Vegetation=rep(c("Low", "High"), 2)), | |
main="d)") | |
mtext("C. parinsignifera", side=3, adj=-8.3, font=3, cex=0.55) | |
text(x=2004, y=0.3, labels="% Canopy", pos=2, cex=0.8, srt=90) | |
# panel h | |
plot.full.model(richness.models$richness.model, link="log", ymax=5, screens=c(9, 10), | |
cols=c(rep(light, 2), rep(dark, 2)), | |
ltys=c(2, 1, 2, 1), | |
labels=list(group1=c("Low", "High"), | |
Vegetation=rep(c("Low", "High"), 2)), | |
main="h)") | |
mtext("Species richness", side=3, adj=-11.5, font=1, cex=0.55) | |
text(x=2004, y=4.8, labels="% Canopy", pos=2, cex=0.8, srt=90) | |
# canopy model | |
for(i in 1:4){ | |
screen(c(2, 3, 7, 8)[i]) | |
plot.partial.model( | |
final.models[[c("limtas.model", "csig.model", "lper.model", "ulae.model")[i]]], | |
col.vals=rep(c(light, dark), 2), | |
line.order=c(1, 3, 2, 4), | |
low.key=c(TRUE, TRUE, FALSE, FALSE)[i], | |
ymax=c(1, 1, 1, 1)[i], | |
lty.vals=c(2, 2, 1, 1), | |
link="logit", | |
main=list(bold=c("b)", "c)", "f)", "g)")[i], | |
italic=c("Lim. tasmaniensis", "C. signifera", "L. peronii", "U. laevigata")[i], | |
offset=c(0.35, 0.22, 0.18, 0.23)[i]), | |
labels=list(group1=c("Rural", "Urban"), | |
"% Canopy"=rep(c("Low", "High"), each=2)), | |
factor.vars=c("pc.urban", "pc.canopy")) | |
} | |
for(i in 1:2){ | |
screen(c(1, 6)[i]) | |
plot.partial.model( | |
final.models[[c("limper.model", "lver.model")[i]]], | |
col.vals=rep(c(light, dark), each=2), | |
lty.vals=c(2, 1, 2, 1), | |
link="logit", | |
ymax=0.2, | |
main=list(bold=c("a)", "e)")[i], italic=c("Lim. peronii", "L. verreauxii")[i], | |
offset=c(0.22, 0.22)[i]), | |
labels=list(group1=c("Rural", "Urban"), | |
Vegetation=rep(c("Low", "High"), 2)), | |
factor.vars=c("pc.urban", "veg.local")) | |
} | |
screen(11) | |
plot(1~1, ann=FALSE, axes=FALSE, type="n") | |
mtext("Year", side=1, cex=0.55, font=2, line=2)#4) | |
mtext("Expected Response", side=2, cex=0.55, font=2, line=1)#2.8) | |
close.screen(all=TRUE) | |
# end plot | |
dev.off() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# This script contains functions for analysing the ACT Frogwatch dataset. They may not be useful or sensible in other contexts. | |
## GLMM FUNCTIONS | |
# calculate effect of habitat suitability, urbanisation, and year for each site | |
binomial.model<- function(x, data){ | |
model<-glmer(formula(x), control=glmerControl(optimizer="bobyqa"), | |
data= data, family=binomial(link="logit")) | |
return(model)} | |
poisson.model<- function(x, data){ | |
model<-glmer(formula(x), control=glmerControl(optimizer="bobyqa"), | |
data= data, family=poisson(link="log")) | |
return(model)} | |
# calculate a specified set of models for a given species | |
urban.models<-function(x, source, model.list, richness){ # species name | |
if(missing(richness))richness<-FALSE | |
# put data in correct order | |
data.thisrun<-source[, c(x, "year", "pc.urban", "veg.local", "pc.canopy", "log.size", "water.type", "site.code", | |
paste("ac", x, sep="."))] | |
colnames(data.thisrun)[c(1, ncol(data.thisrun))]<-c("spp", "ac") | |
# run an appropriate set of models | |
if(richness){ | |
model.results<-lapply(model.list, FUN=function(x, data){poisson.model(x, data)}, data=data.thisrun) | |
}else{ | |
model.results<-lapply(model.list, FUN=function(x, data){binomial.model(x, data)}, data=data.thisrun)} | |
# rank and return | |
model.ranks<-aictab(model.results, modnames=names(model.list)) | |
best.model<-model.results[[as.numeric(rownames(model.ranks)[1])]] | |
return(list(aictab=model.ranks, model=best.model)) | |
} | |
# function to extract usable information from a list of models - called by plot.static.effects() | |
get.coefs<-function(x, model.list, fx.names){ | |
entry<-which(names(model.list)==paste(x, "model", sep=".")) | |
model.thisrun<-model.list[[entry]] | |
coef.matrix<-summary(model.thisrun)$coefficients | |
result<-sapply(fx.names, FUN=function(x){ | |
if(any(rownames(coef.matrix)==x)){ | |
as.numeric(coef.matrix[which(rownames(coef.matrix)==x), 1:2])}else{rep(NA, 2)}}) | |
result<-data.frame(t(result)) | |
colnames(result)<-c("coef", "se") | |
result$variable<-rownames(result) | |
rownames(result)<-NULL | |
return(result)} | |
# get predictions from frogwatch models | |
predict.model<-function(model, continuous.var="year", factor.vars=c("veg.local", "pc.canopy", "pc.urban")){ | |
dataset<-model@frame | |
covar.matrix<-attr(terms(model), "factors") | |
# which of the grouping variables are in this model? | |
keep.selector<-sapply(factor.vars, FUN=function(x){ | |
if(any(rownames(covar.matrix)==x)){return(1)}else{return(0)}}) | |
factor.vars<-factor.vars[which(keep.selector==1)] | |
# set up factors | |
factor.list<-lapply(factor.vars, FUN=function(x, dset){ | |
as.numeric(quantile(dset[, x], c(0.1, 0.9)))}, dset=dataset) | |
names(factor.list)<-factor.vars | |
#quantile(dataset$pc.urban, c(0.1, 0.9)) # check = OK | |
# set up continous variables | |
range.var<-range(dataset[, continuous.var]) | |
variable.list<-list(seq(range.var[1], range.var[2], length.out=100)) | |
names(variable.list)<-continuous.var | |
# merge with factors | |
variable.list<-append(variable.list, factor.list) | |
# use this to predict model effects | |
keep.selector<-sapply(rownames(covar.matrix), FUN=function(x){ | |
if(any(names(variable.list)==x)){return(1)}else{return(0)}}) | |
covar.matrix<-covar.matrix[which(keep.selector==1), ] | |
covar.matrix<-covar.matrix[, which(apply(covar.matrix, 2, sum)>0)] | |
# get coefficients | |
coef.matrix<-fixef(model) | |
keep.selector<-sapply(names(coef.matrix), FUN=function(x){ | |
if(any(colnames(covar.matrix)==x)){return(1)}else{return(0)}}) | |
coef.matrix<-coef.matrix[c(1, which(keep.selector==1))] | |
# make data.frame of new data | |
variables.raw<-expand.grid(variable.list) | |
col.order<-as.numeric(sapply(rownames(covar.matrix), FUN=function(x){ | |
which(colnames(variables.raw)==x)})) | |
variables.raw<-variables.raw[, col.order] | |
# multiply relevant columns together, guided by covar.matrix[-1, ] | |
new.data<-cbind(rep(1, nrow(variables.raw)), | |
apply(covar.matrix, 2, FUN=function(x, dset){ | |
x<-as.numeric(x) | |
cols<-which(x>0) | |
if(length(cols)>1){apply(dset[, cols], 1, prod)}else{dset[, cols]} | |
}, dset=variables.raw)) | |
# multiply by coefficients | |
for(i in 1:ncol(new.data)){new.data[, i]<-new.data[, i]*coef.matrix[i]} | |
variables.raw$fit<-apply(new.data, 1, sum) | |
for(i in 1:length(factor.vars)){ | |
col<-which(colnames(variables.raw)==factor.vars[i]) | |
new.var<-factor(variables.raw[, col], | |
levels=sort(unique(variables.raw[, col])), | |
labels=paste(factor.vars[i], c("Low", "High"), sep="")) | |
variables.raw[, col]<-new.var} | |
# rescale continuous variables as needed | |
if(is.null(attr(dataset[, continuous.var], "scaled:center"))==FALSE){ | |
subtractor<-attr(dataset[, continuous.var], "scaled:center") | |
multiplier<-attr(dataset[, continuous.var], "scaled:scale") | |
variables.raw[, continuous.var]<-(variables.raw[, continuous.var]*multiplier)+subtractor} | |
return(variables.raw) | |
} | |
## PLOT FUNCTIONS | |
# wrapper function to plot static effects as horizontal bar charts | |
plot.static.effects<-function(spp, model.list, fx.names, | |
fx.labels, spp.labels, label.offsets, x.limits, stagger, cols){ # species name | |
if(missing(cols))cols<-c("steelblue4", "black") | |
# extract coefficients | |
result<-lapply(spp, FUN=function(x){get.coefs(x, model.list, fx.names)}) | |
result<-do.call("rbind", result) | |
result$spp<-rep(spp, each=length(fx.names)) | |
result$lci<-result$coef-(1.96*result$se) | |
result$uci<-result$coef+(1.96*result$se) | |
# get into a list, reorder as needed | |
variable.plot<-split(result, result$variable) | |
effect.size<-lapply(variable.plot, FUN=function(x){ | |
effect.size<-sqrt((x$coef/x$se)^2) | |
mean(effect.size, na.rm=TRUE)}) | |
variable.plot<-variable.plot[order(as.numeric(effect.size), decreasing=TRUE)] | |
spp.order<-order(variable.plot[[1]]$coef, decreasing=TRUE) | |
variable.plot<-lapply(variable.plot, FUN=function(x, order.var){x[order.var, ]}, | |
order.var= spp.order) | |
if(missing(fx.labels)==FALSE){names(variable.plot)<-fx.labels} | |
# draw | |
par(mfrow=c(1, length(fx.names)), mar=c(2, 0.5, 2, 0.5), oma=c(1, 5, 0, 0)) | |
for(i in 1:length(variable.plot)){ | |
data.thisrun<-variable.plot[[i]] | |
plot(1~1, type="n", ann=FALSE, axes=FALSE, ylim=c(0.7, (length(spp)+0.2)), | |
xlim= c(min(data.thisrun$lci, na.rm=TRUE), max(data.thisrun$uci, na.rm=TRUE))) | |
axis(1, at=c(-10, 10), labels=NULL, tcl=0) | |
#axis(1, at=seq(-10, 10, 0.5), labels=NA) | |
axis(1, at=seq(x.limits$axis.min[i], x.limits$axis.max[i], x.limits$axis.step[i]), labels=NA, | |
lend="square", tcl=-0.4) | |
if(stagger[i]){ | |
initial.seq<-seq(x.limits$axis.min[i], x.limits$axis.max[i], x.limits$axis.step[i]) | |
odd.sel<-seq(1, 11, 2)[which(seq(1, 11, 2)<=length(initial.seq))] | |
even.sel<-seq(2, 12, 2)[which(seq(2, 12, 2)<=length(initial.seq))] | |
initial.seq<-format(initial.seq, digits=2) | |
odd.seq<-initial.seq[odd.sel] | |
even.seq<-initial.seq[even.sel] | |
axis(1, at=odd.seq, labels=format(odd.seq, digits=2), lwd=0, cex.axis=0.6, line=-0.7) | |
axis(1, at= even.seq, labels=format(even.seq, digits=2),, lwd=0, cex.axis=0.6, line=-0.2) | |
}else{ | |
axis(1, seq(x.limits$axis.min[i], x.limits$axis.max[i], x.limits$axis.step[i]), lwd=0, cex.axis=0.6, line=-0.6)} | |
if(i==1){ | |
axis(2, las=1, at=c(2:nrow(data.thisrun)), labels= spp.labels[spp.order][2:9], | |
lwd=0, font=3, hadj=0, line=3, cex.axis=0.8) | |
axis(2, at=1, labels="Species richness", las=1, hadj=0, lwd=0, line=3, cex.axis=0.8) | |
} | |
if(any(c(2, 4, 6)==i)){col<-cols[1]}else{col<-cols[2]} | |
abline(v=0, lty=2, col=col) | |
for(j in 1:nrow(data.thisrun)){ | |
lines(x=c(data.thisrun$lci[j], data.thisrun$uci[j]), y=rep(j, 2), col=col)} | |
points(x=data.thisrun$coef, y=c(1:nrow(data.thisrun)), pch=15, col=col, cex=0.8) | |
mtext(paste(letters[i], ") ", sep=""), side=3, adj=0.05, cex=0.5, font=2, line=1) | |
mtext(names(variable.plot)[i], side=3, adj=label.offsets[i], cex=0.5, font=1, line=1) | |
} # end for-loop | |
par(mfrow=c(1, 1), mar=c(4, 4, 2, 1), oma=rep(0, 4)) | |
mtext("Variable Coefficient", side=1, adj=0.55, cex=0.6, line=3.0) | |
mtext("Response variable", side=2, adj=0.5, cex=0.6, line=3.2) | |
}# end function | |
# function to draw models with two factors | |
plot.partial.model<-function(model, | |
col.vals, lty.vals, ymax=1, main="Plot Title", line.order, # plot attr | |
labels, factor.vars, link="none", low.key=FALSE # predict vars | |
) | |
{ | |
if(missing(col.vals))col.vals<-rep("black", 4) | |
if(missing(lty.vals))lty.vals<-c(1:4) | |
if(missing(line.order))line.order<-c(1:4) | |
# get inputs | |
predict.dframe<-predict.model(model, continuous.var="year", factor.vars=factor.vars) | |
# add effect of type=still | |
add.val<-as.numeric(fixef(model)[which(names(fixef(model))=="water.typestill")]) | |
predict.dframe$fit<-predict.dframe$fit+add.val | |
# note auto-detection of factor vars (Actually interacting continuous variables) should be possible using terms() | |
if(link=="logit"){predict.dframe$fit<-inv.logit(predict.dframe$fit)} | |
if(link=="log"){predict.dframe$fit<-exp(predict.dframe$fit)} | |
cols<-sapply(colnames(predict.dframe), FUN=function(x, comparison){ | |
if(any(x==comparison)){return(1)}else{return(0)}}, comparison=factor.vars) | |
split.list<-as.list(predict.dframe[, which(cols==1)]) | |
predict.list<-split(predict.dframe, split.list) | |
# set plot attributes | |
yvals<-seq(ymax*0.95, ymax*0.8, length.out=length(predict.list))[line.order] | |
if(low.key){key.top<-ymax*0.3; yvals<-yvals-(ymax-key.top)}else{key.top<-ymax} | |
# draw | |
plot(1~1, xlim=c(range(predict.dframe$year)), ylim=c(0, ymax), type="n", ann=FALSE, axes=FALSE) | |
axis(1) | |
axis(2, las=2); box(bty="l") | |
for(h in 1:length(predict.list)){ | |
lines(predict.list[[h]]$year, predict.list[[h]]$fit, col=col.vals[h], lty= lty.vals[h]) | |
lines(x=quantile(predict.list[[h]]$year, c(0.9, 1)), y=rep(yvals[h], 2), | |
col=col.vals[h], lty= lty.vals[h]) | |
if(missing(labels)){ | |
labels<-names(predict.list) | |
text(x=quantile(predict.list[[h]]$year, 0.9), y=yvals[h], labels=labels[h], pos=2, cex=0.8) | |
}else{ | |
text(x=quantile(predict.list[[h]]$year, 0.9), y=yvals[h], labels=labels[[2]][h], pos=2, cex=0.8) | |
} | |
} # end j | |
if(missing(labels)==FALSE){ | |
text(x=quantile(predict.list[[1]]$year, 0.9), y= key.top, labels=names(labels)[2], pos=2, cex=0.8) | |
text(x=rep(quantile(predict.list[[1]]$year, 0.69), 2), y=c(mean(sort(yvals)[3:4]), mean(sort(yvals)[1:2])), | |
labels=labels[[1]], pos=2, cex=0.8) | |
lines(x=rep(quantile(predict.list[[1]]$year, 0.69), 2), y=sort(yvals)[3:4], lend="square", lwd=2, | |
col=unique(col.vals)[1]) | |
lines(x=rep(quantile(predict.list[[1]]$year, 0.69), 2), y=sort(yvals)[1:2], lend="square", lwd=2, | |
col=unique(col.vals)[2]) | |
} | |
mtext(main[[1]], side=3, adj=0.05, font=2, cex=0.55) | |
mtext(main[[2]], side=3, adj=main[[3]], font=3, cex=0.55) | |
} # end function | |
# function to draw model with three factors | |
plot.full.model<-function(model, screens=c(1, 2), | |
cols=rep("black", 4), ltys=c(1:4), ymax=1, labels, link="none", main="Plot title", low.key=FALSE) | |
{ | |
predict.dframe<-predict.model(model, | |
continuous.var="year", factor.vars=c("veg.local", "pc.canopy", "pc.urban")) | |
# add effect of type=still | |
add.val<-as.numeric(fixef(model)[which(names(fixef(model))=="water.typestill")]) | |
predict.dframe$fit<-predict.dframe$fit+add.val | |
if(link=="logit"){predict.dframe$fit<-inv.logit(predict.dframe$fit)} | |
if(link=="log"){predict.dframe$fit<-exp(predict.dframe$fit)} | |
predict.list<-split(predict.dframe, predict.dframe$pc.urban) | |
predict.list<-lapply(predict.list, FUN=function(x){split(x, list(x$veg.local, x$pc.canopy))}) | |
yvals<-seq(ymax*0.95, ymax*0.8, length.out=4) | |
if(low.key){key.top<-ymax*0.35; yvals<-yvals-(ymax-key.top)}else{key.top<-ymax} | |
if(missing(labels)){labels<-names(predict.list[[1]])} | |
# draw | |
for(i in 1:2){ | |
screen(screens[i]) | |
dataset<-predict.list[[i]] | |
plot(1~1, xlim=c(range(predict.dframe$year)), ylim=c(0, ymax), type="n", ann=FALSE, axes=FALSE) | |
axis(1) | |
if(i==1){ | |
axis(2, las=2); box(bty="l") | |
for(j in 1:4){lines(dataset[[j]]$year, dataset[[j]]$fit, col=cols[j], lty=ltys[j])} | |
mtext(main, side=3, adj=0.05, font=2, cex=0.55) | |
}else{for(j in 1:4){ | |
lines(dataset[[j]]$year, dataset[[j]]$fit, col=cols[j], lty=ltys[j]) | |
lines(x=quantile(dataset[[j]]$year, c(0.8, 1)), y=rep(yvals[j], 2), col=cols[j], lty=ltys[j]) | |
text(x=quantile(dataset[[j]]$year, 0.8), y=yvals[j], labels=labels[[2]][j], pos=2, cex=0.8) | |
} # end j | |
text(x=quantile(dataset[[2]]$year, 0.8), y= key.top, labels=names(labels)[2], pos=2, cex=0.8) | |
text(x=rep(quantile(dataset[[2]]$year, 0.43), 2), y=c(mean(sort(yvals)[3:4]), mean(sort(yvals)[1:2])), | |
labels=labels[[1]], pos=2, cex=0.8) | |
lines(x=rep(quantile(dataset[[2]]$year, 0.45), 2), y=sort(yvals)[3:4], | |
lend="square", lwd=2, col=unique(cols)[1]) | |
lines(x=rep(quantile(dataset[[2]]$year, 0.45), 2), y=sort(yvals)[1:2], | |
lend="square", lwd=2, col=unique(cols)[2]) | |
} # end if i==2 | |
mtext(c("Rural", "Urban")[i], side=1, adj=0.5, font=2, cex=0.5, line=-2) | |
} # end i | |
} # end function |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment