#Analysis section of USGS report library(arm) dataDIR <- "z:/EUSE/Data" ## the data EUSE <- read.csv(paste(dataDIR,"EUSE_USGSReportData.csv", sep="/"), header=T) attach(EUSE) n<-9 REGION.name <- as.vector(REGION) #length(REGION.name)=261 (datapoints)-no Denver outlier uniqREGION <- unique(REGION.name) REGIONindex <- as.numeric(ordered(REGION)) REGIONprecip<-tapply(AnnMeanP,REGION,mean) REGIONtemp<-tapply(AnnMeanT,REGION,mean) REGIONbackag<-rep(NA,n) P.NLCD78<-P.NLCD7+P.NLCD8 for (i in 1:9){ REGIONbackag[i]<-mean(P.NLCD78[REGION.UII<=10®IONindex==i]) } REGIONbackag.cat<-rep(NA,n) for (i in 1:9){ if (REGIONbackag[i]<=50){ REGIONbackag.cat[i]<-0 } else { REGIONbackag.cat[i]<-1 } } ## the models #set urban predictor URB<-P.NLCD2/100 #set eco response variable-- pick one and uncomment it #REGULAR LINEAR REGRESSION #ECO<-NMDS1; ECO.name<-"NMDS1" #ECO<-RICHTOL; ECO.name<-"RICHTOL" #POISSON GENERALIZED LINEAR REGRESSION ECO<-RICH; ECO.name<-"RICH" #ECO<-EPTRICH; ECO.name<-"EPTRICH" ########################################################## ## REGULAR LINEAR REGRESSION-- NMDS1 and RICHTOL #complete pooling lm.pooled <- lm(ECO ~ URB) #one slope and one intercept #no pooling (separate slopes and intercepts) ab.hat.unpooled <- array (NA, c(n,2)) for (j in 1:n){ lm.unpooled <- lm (ECO ~ URB, subset=(REGIONindex==j)) #9 (separate) slopes and 9 (separate) intercepts ab.hat.unpooled[j,] <- summary(lm.unpooled)$coef[,1] } #partial pooling- varying intercept and slope by group #MODEL 1: no group-level predictor M1 <- lmer(ECO ~ URB + (1+URB|REGION)) #MODEL 2: precip group-level predictor M2 <- lmer(ECO ~ URB + REGIONprecip[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION)) #MODEL 3: temp group-level predictor M3 <- lmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONtemp[REGION]+(1+URB|REGION)) #MODEL 4: temp group-level predictor for intercept; precip predictor for slope M4 <- lmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION)) #MODEL 5: antecedent ag+past group-level predictor M5 <- lmer(ECO ~ URB + REGIONbackag[REGION] + URB:REGIONbackag[REGION]+(1+URB|REGION)) #MODEL 6: categorical ag, precip group-level predictor M6 <- lmer(ECO ~ URB + REGIONprecip[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION)+(1+URB|REGIONbackag.cat[REGION])) #MODEL 7: categorical ag, temp group-level predictor M7 <- lmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONtemp[REGION]+(1+URB|REGION)+(1+URB|REGIONbackag.cat[REGION])) #MODEL 8: categorical ag, temp group-level predictor for intercept; precip predictor for slope M8 <- lmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION)+(1+URB|REGIONbackag.cat[REGION])) ########################################################## ## GENERALIZED (POISSON) LINEAR REGRESSION-- RICH and EPTRICH #complete pooling lm.pooled <- glm (ECO ~ URB,family=poisson) #one slope and one intercept #no pooling (separate slopes and intercepts) ab.hat.unpooled <- array (NA, c(n,2)) for (j in 1:n){ lm.unpooled <- glm (ECO ~ URB, subset=(REGIONindex==j),family=poisson) #9 (separate) slopes and 9 (separate) intercepts ab.hat.unpooled[j,] <- summary(lm.unpooled)$coef[,1] } #partial pooling- varying intercept and slope by group #MODEL 1: no group-level predictor M1 <- glmer(ECO ~ URB + (1+URB|REGION),family=poisson) #MODEL 2: precip group-level predictor M2 <- glmer(ECO ~ URB + REGIONprecip[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION),family=poisson) #MODEL 3: temp group-level predictor M3 <- glmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONtemp[REGION]+(1+URB|REGION),family=poisson) #MODEL 4: temp group-level predictor for intercept; precip predictor for slope M4 <- glmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION),family=poisson) #MODEL 5: antecedent ag+past group-level predictor M5 <- glmer(ECO ~ URB + REGIONbackag[REGION] + URB:REGIONbackag[REGION]+(1+URB|REGION),family=poisson) #MODEL 6: categorical ag, precip group-level predictor M6 <- glmer(ECO ~ URB + REGIONprecip[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION)+(1+URB|REGIONbackag.cat[REGION]),family=poisson) #MODEL 7: categorical ag, temp group-level predictor M7 <- glmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONtemp[REGION]+(1+URB|REGION)+(1+URB|REGIONbackag.cat[REGION]),family=poisson) #MODEL 8: categorical ag, temp group-level predictor for intercept; precip predictor for slope M8 <- glmer(ECO ~ URB + REGIONtemp[REGION] + URB:REGIONprecip[REGION]+(1+URB|REGION)+(1+URB|REGIONbackag.cat[REGION]),family=poisson) ########################################################## #MODEL COEFFICIENTS--can use for all ECO response variables, one at a time a.hat.M1 <- coef(M1)$REGION[,1] b.hat.M1 <- coef(M1)$REGION[,2] a.hat.M2 <- coef(M2)$REGION[,1] + coef(M2)$REGION[,3]*REGIONprecip b.hat.M2 <- coef(M2)$REGION[,2] + coef(M2)$REGION[,4]*REGIONprecip a.se.M2 <- se.ranef(M2)$REGION[,1] b.se.M2 <- se.ranef(M2)$REGION[,2] a.hat.M3 <- coef(M3)$REGION[,1] + coef(M3)$REGION[,3]*REGIONtemp b.hat.M3 <- coef(M3)$REGION[,2] + coef(M3)$REGION[,4]*REGIONtemp a.se.M3 <- se.ranef(M3)$REGION[,1] b.se.M3 <- se.ranef(M3)$REGION[,2] a.hat.M4 <- coef(M4)$REGION[,1] + coef(M4)$REGION[,3]*REGIONtemp b.hat.M4 <- coef(M4)$REGION[,2] + coef(M4)$REGION[,4]*REGIONprecip a.se.M4 <- se.ranef(M4)$REGION[,1] b.se.M4 <- se.ranef(M4)$REGION[,2] a.hat.M5 <- coef(M5)$REGION[,1] + coef(M5)$REGION[,3]*REGIONbackag b.hat.M5 <- coef(M5)$REGION[,2] + coef(M5)$REGION[,4]*REGIONbackag a.se.M5 <- se.ranef(M5)$REGION[,1] b.se.M5 <- se.ranef(M5)$REGION[,2] M6.fixef <- fixef(M6) M6.ranef <- ranef(M6) M6.seranef <- se.ranef(M6) a.hat.M6 <- M6.fixef[1] + M6.ranef[[1]][,1] + M6.ranef[[2]][c(1,1,1,2,2,2,1,1,1),1] + M6.fixef[3]*REGIONprecip b.hat.M6 <- M6.fixef[2] + M6.ranef[[1]][,2] + M6.ranef[[2]][c(1,1,1,2,2,2,1,1,1),2] + M6.fixef[4]*REGIONprecip a.se.M6 <- M6.seranef[[1]][,1] b.se.M6 <- M6.seranef[[1]][,2] M7.fixef <- fixef(M7) M7.ranef <- ranef(M7) M7.seranef <- se.ranef(M7) a.hat.M7 <- M7.fixef[1] + M7.ranef[[1]][,1] + M7.ranef[[2]][c(1,1,1,2,2,2,1,1,1),1] + M7.fixef[3]*REGIONtemp b.hat.M7 <- M7.fixef[2] + M7.ranef[[1]][,2] + M7.ranef[[2]][c(1,1,1,2,2,2,1,1,1),2] + M7.fixef[4]*REGIONtemp a.se.M7 <- M7.seranef[[1]][,1] b.se.M7 <- M7.seranef[[1]][,2] M8.fixef <- fixef(M8) M8.ranef <- ranef(M8) M8.seranef <- se.ranef(M8) a.hat.M8 <- M8.fixef[1] + M8.ranef[[1]][,1] + M8.ranef[[2]][c(1,1,1,2,2,2,1,1,1),1] + M8.fixef[3]*REGIONtemp b.hat.M8 <- M8.fixef[2] + M8.ranef[[1]][,2] + M8.ranef[[2]][c(1,1,1,2,2,2,1,1,1),2] + M8.fixef[4]*REGIONprecip a.se.M8 <- M8.seranef[[1]][,1] b.se.M8 <- M8.seranef[[1]][,2] #GRAPHS################################################### #ONE 9-REGION x vs. y--- for regular linear regressions (AV1 or RICHTOL) #M1: no group-level predictors OUT<-"z:/EUSE/USGS Report/Inverts/USGS Final report drafts and comments/Final NMDS1 Figs" postscript(file=paste(OUT, "USGS.M1.eps",sep="/"), width=8, height=8.5, horizontal=F, paper="special") par (mfrow=c(3,3), mar=c(4,4,3,1), oma=c(1,1,2,1)) for (j in 1:n){ plot (URB[REGIONindex==j]*100, ECO[REGIONindex==j], xlim=c(min(URB*100),max(URB*100)), ylim=c(min(ECO[!is.na(ECO)==TRUE]),max(ECO[!is.na(ECO)==TRUE])),xlab="URB", ylab=ECO.name, cex=1.5, cex.lab=1.2, cex.axis=1.2,main=uniqREGION[j],cex.main=1.5) curve (coef(lm.pooled)[1] + coef(lm.pooled)[2]/100*x, lwd=.5, lty=2, add=T) #completely pooled=dashed line curve (ab.hat.unpooled[j,1] + ab.hat.unpooled[j,2]/100*x, lwd=.5, col="gray10", add=T) #unpooled=black line curve (a.hat.M1[j] + b.hat.M1[j]/100*x, lwd=1, col=3, add=T) #partially pooled- green line } dev.off() #ONE 9-REGION x vs. y--- for Poission generalized linear regressions (TR or EPT) #M1: no group-level predictors OUT<-"z:/EUSE/USGS Report/Inverts/USGS Final report drafts and comments/Final EPTRICH Figs" postscript(file=paste(OUT, "USGS.M1.eps",sep="/"), width=8, height=8.5, horizontal=F, paper="special") par (mfrow=c(3,3), mar=c(4,4,3,1), oma=c(1,1,2,1)) for (j in 1:n){ plot (URB[REGIONindex==j]*100, ECO[REGIONindex==j], xlim=c(0,100), ylim=c(min(ECO),max(ECO)),xlab="URB", ylab=ECO.name, cex=1.5, cex.lab=1.2, cex.axis=1.2, main=uniqREGION[j], cex.main=1.5) curve (exp(coef(lm.pooled)[1] + coef(lm.pooled)[2]/100*x), lwd=.5, lty=2, add=T) #completely pooled=dashed line curve (exp(ab.hat.unpooled[j,1] + ab.hat.unpooled[j,2]/100*x), lwd=.5, col="gray10", add=T) #unpooled=black line curve (exp(a.hat.M1[j] + b.hat.M1[j]/100*x), lwd=1, col=3, add=T) #partially pooled- green line } dev.off() #NINE estimated intercept and slopes across REGION-level predictor--for all ECO variables #M2: precip REGION-level predictor lower.aM2 <- a.hat.M2 - a.se.M2 upper.aM2 <- a.hat.M2 + a.se.M2 postscript(file=paste(OUT, "USGS.M2.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONprecip, a.hat.M2, ylim=range(lower.aM2,upper.aM2+0.1), xlab="REGION-level annual precipitation", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19,xlim=c(35,170), main=paste("Multilevel Model 2 (",ECO.name,"), intercept with precipitation")) curve (fixef(M2)["(Intercept)"] + fixef(M2)["REGIONprecip[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONprecip, lower.aM2, REGIONprecip, upper.aM2, lwd=.5, col="gray10") text(x=REGIONprecip,y=a.hat.M2,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4)) lower.bM2 <- b.hat.M2 - b.se.M2 upper.bM2 <- b.hat.M2 + b.se.M2 plot (REGIONprecip, b.hat.M2, ylim=range(lower.bM2,upper.bM2), xlab="REGION-level annual precipitation", ylab=expression(paste("Estimated slope, ", beta[j])), pch=19, xlim=c(35,170),main=paste("Multilevel Model 2(",ECO.name,"), slope with precipitation")) curve (fixef(M2)["URB"] + fixef(M2)["URB:REGIONprecip[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONprecip, lower.bM2, REGIONprecip, upper.bM2, lwd=.5, col="gray10") text(x=REGIONprecip,y=b.hat.M2,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4)) dev.off() #M3: temp REGION-level predictor lower.aM3 <- a.hat.M3 - a.se.M3 upper.aM3 <- a.hat.M3 + a.se.M3 postscript(file=paste(OUT, "USGS.M3.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONtemp, a.hat.M3, ylim=range(lower.aM3,upper.aM3), xlab="REGION-level annual temperature", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19, xlim=c(7,20),main=paste("Multilevel Model 3 (",ECO.name,"), intercept with temperature")) curve (fixef(M3)["(Intercept)"] + fixef(M3)["REGIONtemp[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONtemp, lower.aM3, REGIONtemp, upper.aM3, lwd=.5, col="gray10") text(x=REGIONtemp,y=a.hat.M3,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,2,4)) lower.bM3 <- b.hat.M3 - b.se.M3 upper.bM3 <- b.hat.M3 + b.se.M3 plot (REGIONtemp, b.hat.M3, ylim=range(lower.bM3,upper.bM3), xlab="REGION-level annual temperature", ylab=expression(paste("Estimated slope, ", beta[j])), pch=19, xlim=c(7,20),main=paste("Multilevel Model 3 (",ECO.name,"), slope with temperature")) curve (fixef(M3)["URB"] + fixef(M3)["URB:REGIONtemp[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONtemp, lower.bM3, REGIONtemp, upper.bM3, lwd=.5, col="gray10") text(x=REGIONtemp,y=b.hat.M3,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4)) dev.off() #M4: temp group-level predictor for intercept; precip predictor for slope lower.aM4 <- a.hat.M4 - a.se.M4 upper.aM4 <- a.hat.M4 + a.se.M4 postscript(file=paste(OUT, "USGS.M4.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONtemp, a.hat.M4, ylim=range(lower.aM4,upper.aM4), xlab="REGION-level annual temperature", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19, xlim=c(7,20),main=paste("Multilevel Model 4 (",ECO.name,"), intercept with temperature")) curve (fixef(M4)["(Intercept)"] + fixef(M4)["REGIONtemp[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONtemp, lower.aM4, REGIONtemp, upper.aM4, lwd=.5, col="gray10") text(x=REGIONtemp,y=a.hat.M4,labels=uniqREGION,pos=c(4,4,2,4,4,4,4,4,4)) lower.bM4 <- b.hat.M4 - b.se.M4 upper.bM4 <- b.hat.M4 + b.se.M4 plot (REGIONprecip, b.hat.M4, ylim=range(lower.bM4,upper.bM4), xlab="REGION-level annual precipitation", ylab=expression(paste("Estimated slope, ", beta[j])), pch=19, xlim=c(35,170),main=paste("Multilevel Model 4(",ECO.name,"), slope with precipitation")) curve (fixef(M4)["URB"] + fixef(M4)["URB:REGIONprecip[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONprecip, lower.bM4, REGIONprecip, upper.bM4, lwd=.5, col="gray10") text(x=REGIONprecip,y=b.hat.M4,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4)) dev.off() #M5: antecedent ag+past group-level predictor lower.aM5 <- a.hat.M5 - a.se.M5 upper.aM5 <- a.hat.M5 + a.se.M5 postscript(file=paste(OUT, "USGS.M5.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONbackag, a.hat.M5, ylim=range(lower.aM5,upper.aM5), xlab="REGION-level antecedent ag+pasture", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19,xlim=c(-2,100), main=paste("Multilevel Model 5 (",ECO.name,"), intercept with antecedent ag+pasture")) curve (fixef(M5)["(Intercept)"] + fixef(M5)["REGIONbackag[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONbackag, lower.aM5, REGIONbackag, upper.aM5, lwd=.5, col="gray10") text(x=REGIONbackag,y=a.hat.M5,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,2)) lower.bM5 <- b.hat.M5 - b.se.M5 upper.bM5 <- b.hat.M5 + b.se.M5 plot (REGIONbackag, b.hat.M5, ylim=range(lower.bM5,upper.bM5), xlab="REGION-level antecedent ag+pasture", ylab=expression(paste("Estimated slope, ", beta[j])), pch=19, xlim=c(-2,100),main=paste("Multilevel Model 5(",ECO.name,"), slope with antecedent ag+pasture")) curve (fixef(M5)["URB"] + fixef(M5)["URB:REGIONbackag[REGION]"]*x, lwd=1, col="black", add=TRUE) segments (REGIONbackag, lower.bM5, REGIONbackag, upper.bM5, lwd=.5, col="gray10") text(x=REGIONbackag,y=b.hat.M5,labels=uniqREGION,pos=c(4,4,4,4,4,4,2,4,2)) dev.off() #M6: categorical ag, precip REGION-level predictor lower.aM6 <- a.hat.M6 - a.se.M6 upper.aM6 <- a.hat.M6 + a.se.M6 postscript(file=paste(OUT, "USGS.M6.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONprecip, a.hat.M6, ylim=range(lower.aM6,upper.aM6), xlab="REGION-level annual precipitation", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19, xlim=c(35,170), col=c(4,4,4,2,2,2,4,4,4), main=paste("Multilevel Model 6 (",ECO.name,"), intercept with ag+grassland and precipitation ")) curve (fixef(M6)["(Intercept)"] + M6.ranef[[2]][1,1] + fixef(M6)["REGIONprecip[REGION]"]*x, lwd=1,col="blue", add=TRUE) curve (fixef(M6)["(Intercept)"] + M6.ranef[[2]][2,1] + fixef(M6)["REGIONprecip[REGION]"]*x, lwd=1,col="red", add=TRUE) segments (REGIONprecip, lower.aM6, REGIONprecip, upper.aM6, lwd=.5, col=c(4,4,4,2,2,2,4,4,4)) text(x=REGIONprecip,y=a.hat.M6,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4),col=c(4,4,4,2,2,2,4,4,4)) lower.bM6 <- b.hat.M6 - b.se.M6 upper.bM6 <- b.hat.M6 + b.se.M6 plot (REGIONprecip, b.hat.M6, ylim=range(lower.bM6,upper.bM6), xlab="REGION-level annual precipitation", ylab=expression(paste("Estimated slope, ", beta[j])),pch=19, xlim=c(35,170),col=c(4,4,4,2,2,2,4,4,4),main=paste("Multilevel Model 6 (",ECO.name,"), slope with ag+grassland and precipitation")) curve (fixef(M6)["URB"] + M6.ranef[[2]][1,2] + fixef(M6)["URB:REGIONprecip[REGION]"]*x,lwd=1,col="blue", add=TRUE) curve (fixef(M6)["URB"] + M6.ranef[[2]][2,2] + fixef(M6)["URB:REGIONprecip[REGION]"]*x, lwd=1, col="red", add=TRUE) segments (REGIONprecip, lower.bM6, REGIONprecip, upper.bM6, lwd=.5, col=c(4,4,4,2,2,2,4,4,4)) text(x=REGIONprecip,y=b.hat.M6,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4),col=c(4,4,4,2,2,2,4,4,4)) dev.off() #M7: categorical ag, temp REGION-level predictor lower.aM7 <- a.hat.M7 - a.se.M7 upper.aM7 <- a.hat.M7 + a.se.M7 postscript(file=paste(OUT, "USGS.M7.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONtemp, a.hat.M7, ylim=range(lower.aM7,upper.aM7), xlab="REGION-level annual temperature", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19,xlim=c(7,20), col=c(4,4,4,2,2,2,4,4,4), main=paste("Multilevel Model 7 (",ECO.name,"), intercept with ag+grassland and temperature")) curve (fixef(M7)["(Intercept)"] + M7.ranef[[2]][1,1] + fixef(M7)["REGIONtemp[REGION]"]*x, lwd=1, col="blue", add=TRUE) curve (fixef(M7)["(Intercept)"] + M7.ranef[[2]][2,1] + fixef(M7)["REGIONtemp[REGION]"]*x, lwd=1, col="red", add=TRUE) segments (REGIONtemp, lower.aM7, REGIONtemp, upper.aM7, lwd=.5, col=c(4,4,4,2,2,2,4,4,4)) text(x=REGIONtemp,y=a.hat.M7,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4),col=c(4,4,4,2,2,2,4,4,4)) lower.bM7 <- b.hat.M7 - b.se.M7 upper.bM7 <- b.hat.M7 + b.se.M7 plot (REGIONtemp, b.hat.M7, ylim=range(lower.bM7,upper.bM7), xlab="REGION-level annual temperature", ylab=expression(paste("Estimated slope, ", beta[j])), pch=19, xlim=c(7,20),col=c(4,4,4,2,2,2,4,4,4),main=paste("Multilevel Model 7 (",ECO.name,"), slope with ag+grassland and temperature")) curve (fixef(M7)["URB"] + M7.ranef[[2]][1,2] + fixef(M7)["URB:REGIONtemp[REGION]"]*x, lwd=1, col="blue", add=TRUE) curve (fixef(M7)["URB"] + M7.ranef[[2]][2,2] + fixef(M7)["URB:REGIONtemp[REGION]"]*x, lwd=1, col="red", add=TRUE) segments (REGIONtemp, lower.bM7, REGIONtemp, upper.bM7, lwd=.5, col=c(4,4,4,2,2,2,4,4,4)) text(x=REGIONtemp,y=b.hat.M7,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4),col=c(4,4,4,2,2,2,4,4,4)) dev.off() #M8: categorical ag, temp REGION-level predictor lower.aM8 <- a.hat.M8 - a.se.M8 upper.aM8 <- a.hat.M8 + a.se.M8 postscript(file=paste(OUT, "USGS.M8.eps",sep="/"), width=11.5, height=6, horizontal=F,paper="special") par (mfrow=c(1,2),mgp=c(1.5,0.5,0),tck=-0.02,mar=c(3,3,3,1)) plot (REGIONtemp, a.hat.M8, ylim=range(lower.aM8,upper.aM8), xlab="REGION-level annual temperature", ylab=expression(paste("Estimated intercept, ", alpha[j])), pch=19,xlim=c(7,20), col=c(4,4,4,2,2,2,4,4,4), main=paste("Multilevel Model 8 (",ECO.name,"), intercept with ag+grassland and temperature")) curve (fixef(M8)["(Intercept)"] + M8.ranef[[2]][1,1] + fixef(M8)["REGIONtemp[REGION]"]*x, lwd=1, col="blue", add=TRUE) curve (fixef(M8)["(Intercept)"] + M8.ranef[[2]][2,1] + fixef(M8)["REGIONtemp[REGION]"]*x, lwd=1, col="red", add=TRUE) segments (REGIONtemp, lower.aM8, REGIONtemp, upper.aM8, lwd=.5, col=c(4,4,4,2,2,2,4,4,4)) text(x=REGIONtemp,y=a.hat.M8,labels=uniqREGION,pos=c(4,2,4,4,4,4,4,4,4),col=c(4,4,4,2,2,2,4,4,4)) lower.bM8 <- b.hat.M8 - b.se.M8 upper.bM8 <- b.hat.M8 + b.se.M8 plot (REGIONprecip, b.hat.M8, ylim=range(lower.bM8,upper.bM8), xlab="REGION-level annual precipitation", ylab=expression(paste("Estimated slope, ", beta[j])), pch=19, xlim=c(35,170),col=c(4,4,4,2,2,2,4,4,4),main=paste("Multilevel Model 8 (",ECO.name,"), slope with ag+grassland and precipitation")) curve (fixef(M8)["URB"] + M8.ranef[[2]][1,2] + fixef(M8)["URB:REGIONprecip[REGION]"]*x, lwd=1, col="blue", add=TRUE) curve (fixef(M8)["URB"] + M8.ranef[[2]][2,2] + fixef(M8)["URB:REGIONprecip[REGION]"]*x, lwd=1, col="red", add=TRUE) segments (REGIONprecip, lower.bM8, REGIONprecip, upper.bM8, lwd=.5, col=c(4,4,4,2,2,2,4,4,4)) text(x=REGIONprecip,y=b.hat.M8,labels=uniqREGION,pos=c(4,4,4,4,4,4,4,4,4),col=c(4,4,4,2,2,2,4,4,4)) dev.off()