mjwestgate/Frogwatch_2015_analysis.R

## Frogwatch_2015_analysis.R
# 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()

## Frogwatch_2015_functions.R
# 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
	# 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 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(ymax0.95, ymax0.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(ymax0.95, ymax0.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