Load libraries
library(nlme)
library(ggplot2)
library(MASS)
library(stargazer)
library(ezknitr)
Read in data
mwdat.t<-read.csv("data/MooseWolfDataUpdated.csv")
mwdat.t<-subset(mwdat.t, Year!=1996)
n<-nrow(mwdat.t)
Create a survey variable to indicate the type of survey
mwdat.t$survey=rep("heli", n)
mwdat.t$survey[mwdat.t$Year <= 1997] <- "fixed.early"
mwdat.t$survey[mwdat.t$Year > 1997 & mwdat.t$Year <2004 ] <-"fixed.late"
Calculate lambdas
mwdat.t$lambda<-rep(NA, n)
mwdat.t$lambda[2:n]<-mwdat.t$moose[2:n]/mwdat.t$moose[1:(n-1)]
mwdat.t$log.lam<-log(mwdat.t$lambda)
Drop first and last observations since they won't be used in the regression due to missing values
mwdat<-mwdat.t[-c(1,n),]
mwdat$yr<-mwdat$Year-1983
Center wolf variables
mwdat$pvwolf<-scale(mwdat$prvwfden)
Fit full model with AR1 correlation structure.
ca.totmod<-gls(ca.total~pvwolf+survey+survey:pvwolf-1-pvwolf, correlation = corAR1(value=0.2, form=~yr),
weights = varPower(form=~fitted(.)), data=mwdat, method="ML", na.action=na.omit)
Use backwards AIC to choose the model. This suggests the full model is “best”.
stepAIC(ca.totmod)
## Start: AIC=-118.34
## ca.total ~ pvwolf + survey + survey:pvwolf - 1 - pvwolf
## Error in dropterm.default(fit, scope$drop, scale = scale, trace = max(0, : number of rows in use has changed: remove missing values?
Refit the full model using gls with REML and means coding
ca.totmod<-gls(ca.total~pvwolf+survey+survey:pvwolf-1-pvwolf, correlation = corAR1(value=0.2, form=~yr),
weights = varPower(form=~fitted(.)), data=mwdat, method="REML", na.action=na.omit)
summary(ca.totmod)
## Generalized least squares fit by REML
## Model: ca.total ~ pvwolf + survey + survey:pvwolf - 1 - pvwolf
## Data: mwdat
## AIC BIC logLik
## -77.68388 -66.714 47.84194
##
## Correlation Structure: ARMA(1,0)
## Formula: ~yr
## Parameter estimate(s):
## Phi1
## 0.2062438
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 2.210641
##
## Coefficients:
## Value Std.Error t-value p-value
## surveyfixed.early 0.18879144 0.01189336 15.873681 0.0000
## surveyfixed.late 0.24619673 0.05295091 4.649528 0.0001
## surveyheli 0.15663509 0.00861658 18.178334 0.0000
## pvwolf:surveyfixed.early -0.02497013 0.01275624 -1.957483 0.0615
## pvwolf:surveyfixed.late 0.02051953 0.08176722 0.250951 0.8039
## pvwolf:surveyheli -0.01078373 0.00570820 -1.889163 0.0705
##
## Correlation:
## srvyfxd.r srvyfxd.l srvyhl pvwlf:srvyfxd.r
## surveyfixed.late 0.050
## surveyheli 0.002 0.031
## pvwolf:surveyfixed.early -0.096 0.029 0.001
## pvwolf:surveyfixed.late 0.036 0.869 0.013 0.021
## pvwolf:surveyheli -0.001 -0.022 -0.645 -0.001
## pvwlf:srvyfxd.l
## surveyfixed.late
## surveyheli
## pvwolf:surveyfixed.early
## pvwolf:surveyfixed.late
## pvwolf:surveyheli -0.009
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.97289590 -0.51910891 -0.05751342 0.38978056 2.10920424
##
## Residual standard error: 1.345614
## Degrees of freedom: 31 total; 25 residual
Calculate fitted values and CI for the regression line
V<-vcov(ca.totmod) # variance covariance matrix of beta^
X<-model.matrix(ca.totmod, data=mwdat[is.na(mwdat$ca.total)!=TRUE,])
mwdat$ca.total.hat<-mwdat$ca.total.se<-rep(NA, nrow(mwdat))
ind<-is.na(mwdat$ca.total)!=TRUE
mwdat$ca.total.hat[ind]<-X%*%coef(ca.totmod)
mwdat$ca.total.se[ind]<-sqrt(diag(X%*%V%*%t(X)))
mwdat$lci<-mwdat$ca.total.hat-1.96*mwdat$ca.total.se
mwdat$uci<-mwdat$ca.total.hat+1.96*mwdat$ca.total.se
Now, plot the results
cols<-c("Red", "Blue", "Black")
#windows(width=8, height=6)
par(bty="L", cex.main=1.8, cex.lab=1.8, cex.axis=1.8, mar=c(6.1, 8.7, 4.1, 1.5), lwd=2)
xo<-order(mwdat$ca.total.hat)
with(mwdat, plot(prvwfden, ca.total,col=cols[as.numeric(as.factor(mwdat$survey))], pch=16, xlab="",
ylab="",cex=2))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden, ca.total.hat, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden, lci, lty=2, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden, uci, lty=2, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, ca.total.hat, col=cols[2]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, lci, lty=2, col=cols[2]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, uci, lty=2,, col=cols[2]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, ca.total.hat, col=cols[3]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, lci, lty=2, col=cols[3]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, uci, lty=2, col=cols[3]))
mtext(side=2, "Calf:Population (t+1)", cex=1.8, line=3)
mtext(side=1, "Wolves (t)", cex=1.8, line=3.2)
Residuals depend on survey method?
boxplot(resid(ca.totmod, type="normalized")~ mwdat$survey[ind])
Residual versus fitted value plot
plot(ca.totmod, pch=16)
Look at residuals over time
plot(mwdat$Year[ind],resid(ca.totmod, type="normalized"), pch=16, xlab="Year", ylab="Residual")
abline(h=0)
Fit full model. In this case, note that there is a missing value for one of the years, so it is better to use the continuous time AR1 model.
lambda.mod<-gls(log.lam~pvwolf+survey+survey:pvwolf, correlation = corCAR1(value=0.2, form=~yr),
data=mwdat[is.na(mwdat$log.lam)!=T,], method="ML", na.action=na.omit)
Use backwards AIC to choose the model.
lambda.mod.r<-stepAIC(lambda.mod)
## Start: AIC=8.22
## log.lam ~ pvwolf + survey + survey:pvwolf
## Error in dropterm.default(fit, scope$drop, scale = scale, trace = max(0, : number of rows in use has changed: remove missing values?
Refit using REML and means parameterization
lambda.mod<-gls(log.lam~pvwolf, correlation = corCAR1(value=0.2, form=~yr),
data=mwdat, na.action=na.omit, method="REML")
summary(lambda.mod)
## Generalized least squares fit by REML
## Model: log.lam ~ pvwolf
## Data: mwdat
## AIC BIC logLik
## 10.34238 15.6712 -1.17119
##
## Correlation Structure: Continuous AR(1)
## Formula: ~yr
## Parameter estimate(s):
## Phi
## 7.379663e-11
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.03066080 0.04077256 -0.7519959 0.4583
## pvwolf -0.09101867 0.04013192 -2.2679866 0.0312
##
## Correlation:
## (Intr)
## pvwolf 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.7216839 -0.7705192 -0.1346553 0.7316979 2.0428126
##
## Residual standard error: 0.2233205
## Degrees of freedom: 30 total; 28 residual
Residual plots
plot(lambda.mod, pch=16)
Residuals depend on survey method?
ind2<-is.na(mwdat$log.lam)!=T
boxplot(resid(lambda.mod, type="normalized")~ mwdat$survey[ind2])
Residual versus fitted value plot
plot(lambda.mod, pch=16)
Look at residuals over time
plot(mwdat$Year[ind2],resid(lambda.mod, type="normalized"), pch=16, xlab="Year", ylab="Residual")
abline(h=0)
Calculate fitted values and CI for the regression line
V<-vcov(lambda.mod) # variance covariance matrix of beta^
X<-model.matrix(lambda.mod, data=mwdat[ind2,])
mwdat$lambda.hat<-mwdat$ca.total.se<-rep(NA, nrow(mwdat))
mwdat$lambda.hat[ind2]<-X%*%coef(lambda.mod)
mwdat$lambda.se[ind2]<-sqrt(diag(X%*%V%*%t(X)))
mwdat$lambda.lci<-mwdat$lambda.hat-1.96*mwdat$lambda.se
mwdat$lambda.uci<-mwdat$lambda.hat+1.96*mwdat$lambda.se
Plot the results
cols<-c("Red", "Blue", "Black")
#windows(width=8, height=6)
par(bty="L", cex.main=1.8, cex.lab=1.8, cex.axis=1.8, mar=c(6.1, 8.7, 4.1, 1.5), lwd=2)
xo<-order(mwdat$lambda.hat)
with(mwdat, plot(prvwfden,log.lam,col=cols[as.numeric(as.factor(mwdat$survey))], pch=16, xlab="",
ylab="",cex=2))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden,lambda.hat, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden, lambda.lci, lty=2, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden, lambda.uci, lty=2, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, lambda.hat, col=cols[2]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, lambda.lci, lty=2, col=cols[2]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, lambda.uci, lty=2,, col=cols[2]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, lambda.hat, col=cols[3]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, lambda.lci, lty=2, col=cols[3]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, lambda.uci, lty=2, col=cols[3]))
mtext(side=2, expression(paste("Moose log(", lambda[t], ")")), cex=1.8, line=3)
mtext(side=1, "Wolves (t)", cex=1.8, line=3.2)
Now, both, but without the CIs
Without the CIs
#windows(width=16, height=8)
par(mfrow=c(1,2), bty="L", cex.main=1.8, cex.lab=1.8, cex.axis=1.8, mar=c(6.1, 8.7, 4.1, 1.5), lwd=2)
xo<-order(mwdat$ca.total.hat)
with(mwdat, plot(prvwfden, ca.total,col=cols[as.numeric(as.factor(mwdat$survey))], pch=16, xlab="",
ylab="",cex=2))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden, ca.total.hat, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, ca.total.hat, col=cols[2]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, ca.total.hat, col=cols[3]))
mtext(side=2, "Calf:Population (t+1)", cex=1.8, line=3)
mtext(side=1, "Wolves (t)", cex=1.8, line=3.2)
xo<-order(mwdat$lambda.hat)
with(mwdat, plot(prvwfden,log.lam,col=cols[as.numeric(as.factor(mwdat$survey))], pch=16, xlab="",
ylab="",cex=2))
with(subset(mwdat[xo,], survey=="fixed.early"), lines(prvwfden,lambda.hat, col=cols[1]))
with(subset(mwdat[xo,], survey=="fixed.late"), lines(prvwfden, lambda.hat, col=cols[2]))
with(subset(mwdat[xo,], survey=="heli"), lines(prvwfden, lambda.hat, col=cols[3]))
mtext(side=2, expression(paste("Moose log(", lambda[t], ")")), cex=1.8, line=3)
mtext(side=1, "Wolves (t)", cex=1.8, line=3.2)
Table of regression results
stargazer(ca.totmod, lambda.mod, type="text")
##
## =====================================================
## Dependent variable:
## ----------------------------
## ca.total log.lam
## (1) (2)
## -----------------------------------------------------
## surveyfixed.early 0.189***
## (0.012)
##
## surveyfixed.late 0.246***
## (0.053)
##
## surveyheli 0.157***
## (0.009)
##
## pvwolf:surveyfixed.early -0.025*
## (0.013)
##
## pvwolf:surveyfixed.late 0.021
## (0.082)
##
## pvwolf:surveyheli -0.011*
## (0.006)
##
## pvwolf -0.091**
## (0.040)
##
## Constant -0.031
## (0.041)
##
## -----------------------------------------------------
## Observations 31 30
## Log Likelihood 47.842 -1.171
## Akaike Inf. Crit. -77.684 10.342
## Bayesian Inf. Crit. -66.714 15.671
## =====================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
Add back data from 1983 for plot & write out data
mwdat.1983<-c(as.data.frame(mwdat.t[1,]), rep(NA, 10))
names(mwdat.1983)<-names(mwdat)
mwdat2<-rbind(mwdat.1983, mwdat)
Write out data with regression results for plotting
write.csv(mwdat2, file="data/moosewolf.csv")
sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows >= 8 x64 (build 9200)
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] stargazer_5.2 MASS_7.3-45 ggplot2_2.2.1 nlme_3.1-131 knitr_1.16
## [6] ezknitr_0.6
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.12 lattice_0.20-34 mime_0.5
## [4] R.methodsS3_1.7.1 plyr_1.8.4 grid_3.3.3
## [7] gtable_0.2.0 magrittr_1.5 evaluate_0.10
## [10] scales_0.4.1 highr_0.6 rlang_0.1.2
## [13] stringi_1.1.5 lazyeval_0.2.0 R.oo_1.21.0
## [16] R.utils_2.5.0 tools_3.3.3 stringr_1.2.0
## [19] markdown_0.7.7 munsell_0.4.3 colorspace_1.3-2
## [22] tibble_1.3.3
Spun using: ezspin(“Scripts/RegressionRevisitedMinus1996.R”, out_dir = “output”, fig_dir = “figures”, keep_md=FALSE)