Load libraries
library(nlme)
library(ggplot2)
library(MASS)
library(stargazer)
library(ezknitr)
Read in data
mwdat.t<-read.csv("data/MooseWolfDataUpdated.csv")
n<-nrow(mwdat.t)
Correlation between calf:pop and calf:cow 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=-122.95
## ca.total ~ pvwolf + survey + survey:pvwolf - 1 - pvwolf
##
## Df AIC
## <none> -122.95
## - pvwolf:survey 3 -120.47
## Generalized least squares fit by maximum likelihood
## Model: ca.total ~ pvwolf + survey + survey:pvwolf - 1 - pvwolf
## Data: mwdat
## Log-likelihood: 70.47639
##
## Coefficients:
## surveyfixed.early surveyfixed.late surveyheli
## 0.191653986 0.235734587 0.154813701
## pvwolf:surveyfixed.early pvwolf:surveyfixed.late pvwolf:surveyheli
## -0.028783109 0.002603613 -0.008770150
##
## Correlation Structure: AR(1)
## Formula: ~yr
## Parameter estimate(s):
## Phi
## -0.004219534
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 2.141391
## Degrees of freedom: 32 total; 26 residual
## Residual standard error: 1.018344
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
## -81.85688 -70.53401 49.92844
##
## Correlation Structure: AR(1)
## Formula: ~yr
## Parameter estimate(s):
## Phi
## 0.1840362
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 2.153253
##
## Coefficients:
## Value Std.Error t-value p-value
## surveyfixed.early 0.19078171 0.01134441 16.817237 0.0000
## surveyfixed.late 0.24493681 0.05013372 4.885670 0.0000
## surveyheli 0.15649189 0.00831678 18.816399 0.0000
## pvwolf:surveyfixed.early -0.02650804 0.01228918 -2.157024 0.0404
## pvwolf:surveyfixed.late 0.01845220 0.07809078 0.236292 0.8151
## pvwolf:surveyheli -0.01039358 0.00546959 -1.900247 0.0685
##
## Correlation:
## srvyfxd.r srvyfxd.l srvyhl pvwlf:srvyfxd.r
## surveyfixed.late 0.044
## surveyheli 0.001 0.027
## pvwolf:surveyfixed.early -0.085 0.027 0.001
## pvwolf:surveyfixed.late 0.031 0.870 0.011 0.019
## pvwolf:surveyheli -0.001 -0.019 -0.645 -0.001
## pvwlf:srvyfxd.l
## surveyfixed.late
## surveyheli
## pvwolf:surveyfixed.early
## pvwolf:surveyfixed.late
## pvwolf:surveyheli -0.008
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.05214472 -0.53933733 -0.05859649 0.41346540 2.16441870
##
## Residual standard error: 1.187416
## Degrees of freedom: 32 total; 26 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)
In response to reviewer no. 1, fit and summarize the ANCOVA model
ca.ancova<-gls(ca.total~pvwolf+survey, correlation = corAR1(value=0.2, form=~yr),
weights = varPower(form=~fitted(.)), data=mwdat, method="ML", na.action=na.omit)
summary(ca.ancova)
## Generalized least squares fit by maximum likelihood
## Model: ca.total ~ pvwolf + survey
## Data: mwdat
## AIC BIC logLik
## -124.2795 -114.0193 69.13975
##
## Correlation Structure: AR(1)
## Formula: ~yr
## Parameter estimate(s):
## Phi
## 0.0915449
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 2.095532
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.19429995 0.010641331 18.258989 0.0000
## pvwolf -0.01133751 0.004846058 -2.339532 0.0267
## surveyfixed.late 0.03481484 0.024663577 1.411589 0.1691
## surveyheli -0.03687112 0.013417126 -2.748064 0.0104
##
## Correlation:
## (Intr) pvwolf srvyf.
## pvwolf 0.088
## surveyfixed.late -0.410 0.061
## surveyheli -0.825 -0.435 0.307
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.2205499 -0.6152457 -0.1263926 0.5695437 2.2689521
##
## Residual standard error: 0.9888305
## Degrees of freedom: 32 total; 28 residual
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")
Use backwards AIC to choose the model.
lambda.mod.r<-stepAIC(lambda.mod)
## Start: AIC=19.32
## log.lam ~ pvwolf + survey + survey:pvwolf
##
## Df AIC
## - pvwolf:survey 2 16.738
## <none> 19.320
##
## Step: AIC=16.74
## log.lam ~ pvwolf + survey
##
## Df AIC
## - survey 2 13.364
## <none> 16.738
## - pvwolf 1 19.789
##
## Step: AIC=13.36
## log.lam ~ pvwolf
##
## Df AIC
## <none> 13.364
## - pvwolf 1 16.080
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
## 21.85048 27.31967 -6.925241
##
## Correlation Structure: Continuous AR(1)
## Formula: ~yr
## Parameter estimate(s):
## Phi
## 1.404965e-10
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.02974038 0.04899329 -0.6070296 0.5486
## pvwolf -0.10532322 0.04824866 -2.1829255 0.0373
##
## Correlation:
## (Intr)
## pvwolf 0.001
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.6051042 -0.7054271 -0.1015936 0.7244349 1.9402045
##
## Residual standard error: 0.272783
## Degrees of freedom: 31 total; 29 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, ca.ancova, lambda.mod, type="text")
##
## ======================================================
## Dependent variable:
## -----------------------------
## ca.total log.lam
## (1) (2) (3)
## ------------------------------------------------------
## surveyfixed.early 0.191***
## (0.011)
##
## pvwolf -0.011** -0.105**
## (0.005) (0.048)
##
## surveyfixed.late 0.245*** 0.035
## (0.050) (0.025)
##
## surveyheli 0.156*** -0.037***
## (0.008) (0.013)
##
## pvwolf:surveyfixed.early -0.027**
## (0.012)
##
## pvwolf:surveyfixed.late 0.018
## (0.078)
##
## pvwolf:surveyheli -0.010*
## (0.005)
##
## Constant 0.194*** -0.030
## (0.011) (0.049)
##
## ------------------------------------------------------
## Observations 32 32 31
## Log Likelihood 49.928 69.140 -6.925
## Akaike Inf. Crit. -81.857 -124.280 21.850
## Bayesian Inf. Crit. -70.534 -114.019 27.320
## ======================================================
## 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/RegressionRevisited.R”, out_dir = “output”, fig_dir = “figures”, keep_md=FALSE)