update the inst/doc folder with 3.11.4 files

Author: eeholmes

inst/doc/Chapter_AnimalTracking.R

@@ -229,20 +229,7 @@ levels(loggerheadNoisy$turtle)
 
 
 ###################################################
-### code chunk number 19: Cs17_kftrack (eval = FALSE)
-###################################################
-## library(kftrack) # must be installed from github
-## loggerhead <- loggerhead
-## turtlename <- "BigMama"
-## dat <- loggerhead[which(loggerhead$turtle == turtlename), 2:6]
-## model <- kftrack(dat,
-##   fix.first = FALSE, fix.last = FALSE,
-##   var.struct = "uniform"
-## )
-
-
-###################################################
-### code chunk number 21: Cs18_code
+### code chunk number 20: Cs18_code
 ###################################################
 ###############################################################
 #    GCDF FUNCTION

inst/doc/Chapter_CombiningTrendData.R

@@ -8,83 +8,87 @@ library(MARSS)
 ### code chunk number 3: Cs001_readinredddata
 ###################################################
 head(okanaganRedds)
-logRedds = log(t(okanaganRedds)[2:3,])
+logRedds <- log(t(okanaganRedds)[c("aerial", "ground"), ])
 
 
 ###################################################
 ### code chunk number 4: Cs002_fig1
 ###################################################
 # Code for plotting raw Okanagan redd counts
-plot(okanaganRedds[,1], okanaganRedds[,2],
-    xlab = "", ylab="Redd counts",main="", col="red", pch=1)
-points(okanaganRedds[,1], okanaganRedds[,3], col="blue", pch=2)
-legend('topleft', inset=0.1, legend=c("Aerial survey","Ground survey"),
- col=c("red","blue"), pch=c(1,2))
+plot(okanaganRedds[, 1], okanaganRedds[, 2],
+  xlab = "", ylab = "Redd counts", main = "", col = "red", pch = 1
+)
+points(okanaganRedds[, 1], okanaganRedds[, 3], col = "blue", pch = 2)
+legend("topleft",
+  inset = 0.1, legend = c("Aerial survey", "Ground survey"),
+  col = c("red", "blue"), pch = c(1, 2)
+)
 
 
 ###################################################
 ### code chunk number 5: Cs003_reddmodel1
 ###################################################
-model1=list()
-model1$R="diagonal and equal"
-model1$Z=matrix(1,2,1)
-model1$A="scaling" 
-kem1 = MARSS(logRedds, model=model1)
+model1 <- list()
+model1$R <- "diagonal and equal"
+model1$Z <- matrix(1, 2, 1)
+model1$A <- "scaling"
+kem1 <- MARSS(logRedds, model = model1)
 
 
 ###################################################
 ### code chunk number 6: Cs004_reddmodel2
 ###################################################
-model2=model1 #model2 is based on model1
-model2$R="diagonal and unequal"
-kem2 = MARSS(logRedds, model=model2)
+model2 <- model1 # model2 is based on model1
+model2$R <- "diagonal and unequal"
+kem2 <- MARSS(logRedds, model = model2)
 
 
 ###################################################
 ### code chunk number 7: Cs005_reddmodel3
 ###################################################
-model3=list()
-model3$Q="diagonal and equal"
-model3$R="diagonal and equal"
-model3$U="equal"
-model3$Z="identity"
-model3$A="zero"
-kem3 = MARSS(logRedds, model=model3)
+model3 <- list()
+model3$Q <- "diagonal and equal"
+model3$R <- "diagonal and equal"
+model3$U <- "equal"
+model3$Z <- "identity"
+model3$A <- "zero"
+kem3 <- MARSS(logRedds, model = model3)
 
 
 ###################################################
 ### code chunk number 8: Cs005b_aic
 ###################################################
-c(mod1=kem1$AICc, mod2=kem2$AICc, mod3=kem3$AICc)
+c(mod1 = kem1$AICc, mod2 = kem2$AICc, mod3 = kem3$AICc)
 
 
 ###################################################
 ### code chunk number 9: Cs006_fig2
 ###################################################
 # Code for plotting the fit from the best model
-plot(okanaganRedds[,1], logRedds[1,],
-xlab = "", ylab="Redd counts",main="", col="red", ylim=c(0,8))
-points(okanaganRedds[,1], logRedds[2,], col="blue", pch=2)
-lines(okanaganRedds[,1], c(kem1$states), lty=1, lwd=2)
-lines(okanaganRedds[,1], c(kem1$states + 2*kem1$states.se), lty=1, lwd=1, col="grey40")
-lines(okanaganRedds[,1], c(kem1$states - 2*kem1$states.se), lty=1, lwd=1, col="grey40")
+plot(okanaganRedds[, 1], logRedds[1, ], xlab = "", ylab = "Redd counts", 
+     main = "", col = "red", ylim = c(0, 8)
+)
+points(okanaganRedds[, 1], logRedds[2, ], col = "blue", pch = 2)
+lines(okanaganRedds[, 1], c(kem1$states), lty = 1, lwd = 2)
+lines(okanaganRedds[, 1], c(kem1$states + 2 * kem1$states.se), lty = 1, lwd = 1, col = "grey40")
+lines(okanaganRedds[, 1], c(kem1$states - 2 * kem1$states.se), lty = 1, lwd = 1, col = "grey40")
 
 
 ###################################################
 ### code chunk number 26: Cs007_birddata
 ###################################################
-birddat = t(kestrel[,2:4]) 
+birddat <- t(kestrel[, c("British.Columbia", "Alberta", "Saskatchewan")])
 head(kestrel)
 
 
 ###################################################
 ### code chunk number 27: Cs008_plot-bird-data
 ###################################################
 # Make a plot of the three time series
-plot(kestrel[,1], kestrel[,2], xlab = "", ylab="Index of kestrel abundance",main="", col="red",ylim=c(0,2), pch=21)
-points(kestrel[,1], kestrel[,3], col="blue", pch=22)
-points(kestrel[,1], kestrel[,4], col="purple", pch=25)
-legend('topright',inset=0.1, legend=c("British Columbia","Alberta","Saskatchewan"), col=c("red","blue","purple"), pch=c(21,22,25))
+plot(kestrel[, 1], kestrel[, 2], xlab = "", ylab = "Index of kestrel abundance", main = "", col = "red", ylim = c(0, 2), pch = 21)
+points(kestrel[, 1], kestrel[, 3], col = "blue", pch = 22)
+points(kestrel[, 1], kestrel[, 4], col = "purple", pch = 25)
+legend("topright", inset = 0.1, legend = c("British Columbia", "Alberta", "Saskatchewan"), col = c("red", "blue", "purple"), pch = c(21, 22, 25))
 
 
 ###################################################
@@ -99,53 +103,56 @@ kem.b1 = MARSS(birddat, model=model.b1, control=list(minit=100) )
 ###################################################
 ### code chunk number 29: Cs011_fit-bird-model-2
 ###################################################
-model.b2=list()
-model.b2$Q="diagonal and equal"
-model.b2$R="diagonal and equal"
-model.b2$Z="identity"
-model.b2$A="zero"
-model.b2$U="equal"
-kem.b2 = MARSS(birddat, model=model.b2)
+model.b2 <- list()
+model.b2$Q <- "diagonal and equal"
+model.b2$R <- "diagonal and equal"
+model.b2$Z <- "identity"
+model.b2$A <- "zero"
+model.b2$U <- "equal"
+kem.b2 <- MARSS(birddat, model = model.b2)
 
 
 ###################################################
 ### code chunk number 30: Cs013_fit-bird-model-3
 ###################################################
-model.b3=model.b2 #is is based on model.b2
-#all we change is the structure of Q
-model.b3$Q="diagonal and unequal"
-model.b3$U="unequal"
-kem.b3 = MARSS(birddat, model=model.b3)
+model.b3 <- model.b2 # is is based on model.b2
+# all we change is the structure of Q
+model.b3$Q <- "diagonal and unequal"
+model.b3$U <- "unequal"
+kem.b3 <- MARSS(birddat, model = model.b3)
 
 
 ###################################################
 ### code chunk number 31: Cs015_fit-bird-model-4
 ###################################################
-model.b4=list()
-model.b4$Q="diagonal and unequal"
-model.b4$R="diagonal and equal"
-model.b4$Z=factor(c("BC","AB-SK","AB-SK"))
-model.b4$A="scaling"
-model.b4$U="unequal"
-kem.b4 = MARSS(birddat, model=model.b4)
+model.b4 <- list()
+model.b4$Q <- "diagonal and unequal"
+model.b4$R <- "diagonal and equal"
+model.b4$Z <- factor(c("BC", "AB-SK", "AB-SK"))
+model.b4$A <- "scaling"
+model.b4$U <- "unequal"
+kem.b4 <- MARSS(birddat, model = model.b4)
 
 
 ###################################################
 ### code chunk number 32: Cs016_aics
 ###################################################
-c(mod1=kem.b1$AICc, mod2=kem.b2$AICc, mod3=kem.b3$AICc, mod4=kem.b4$AICc)
+c(mod1 = kem.b1$AICc, mod2 = kem.b2$AICc, mod3 = kem.b3$AICc, mod4 = kem.b4$AICc)
 
 
 ###################################################
 ### code chunk number 33: Cs017_plot-bird-model-4-fits
 ###################################################
 # Make a plot of the predicted trajectory, confidence intervals, and the raw data in log-space
-plot(kestrel[,1], kestrel[,2], xlab = "", ylab="Index of kestrel abundance",main="", col="red", ylim=c(0,2), pch=21)
-points(kestrel[,1], kestrel[,3], col="blue", pch=22)
-points(kestrel[,1], kestrel[,4], col="purple", pch=25)
-lines(kestrel[,1], c(kem.b4$states[1,]), lty=3, lwd=2, col="red")
-lines(kestrel[,1], c(kem.b4$states[2,]), lty=3, lwd=2, col="blue")
-lines(kestrel[,1], c(kem.b4$states[2,]+coef(kem.b4,type="matrix")$A[3,1]), lty=3, lwd=2, col="purple")
-legend('topright',inset=0.1, legend=c("British Columbia","Alberta","Saskatchewan"), col=c("red","blue","purple"), pch=c(21,22,25))
+plot(kestrel[, 1], kestrel[, 2], xlab = "", ylab = "Index of kestrel abundance", 
+     main = "", col = "red", ylim = c(0, 2), pch = 21)
+points(kestrel[, 1], kestrel[, 3], col = "blue", pch = 22)
+points(kestrel[, 1], kestrel[, 4], col = "purple", pch = 25)
+lines(kestrel[, 1], c(kem.b4$states[1, ]), lty = 3, lwd = 2, col = "red")
+lines(kestrel[, 1], c(kem.b4$states[2, ]), lty = 3, lwd = 2, col = "blue")
+lines(kestrel[, 1], c(kem.b4$states[2, ] + coef(kem.b4, type = "matrix")$A[3, 1]), 
+      lty = 3, lwd = 2, col = "purple")
+legend("topright", inset = 0.1, legend = c("British Columbia", "Alberta", "Saskatchewan"), 
+       col = c("red", "blue", "purple"), pch = c(21, 22, 25))
 
 

inst/doc/Chapter_Covariates.R

@@ -22,7 +22,6 @@ covariates <- rbind(
   Temp = fulldat[years, "Temp"],
   TP = fulldat[years, "TP"]
 )
-
 # z.score the covariates
 covariates <- zscore(covariates)
 

inst/doc/Chapter_KFAS.R

@@ -174,7 +174,7 @@ cbind(
 ###################################################
 ### code chunk number 19: Cs301_obs-filtering
 ###################################################
-f_kfas <- KFS(fit_kfas$model,
+kf_kfas <- KFS(fit_kfas$model,
   filtering = "signal",
   smoothing = "signal", simplify = FALSE
 )
@@ -183,7 +183,7 @@ f_kfas <- KFS(fit_kfas$model,
 ###################################################
 ### code chunk number 20: Cs302_obs-filtering
 ###################################################
-f_marss <- MARSSkf(fit_marss)
+kf_marss <- MARSSkf(fit_marss)
 
 
 ###################################################
@@ -197,7 +197,7 @@ n <- 10
 ###################################################
 ytt1_fit <- fitted(fit_marss, type = "ytt1")$.fitted
 ytt1_hatyt <- MARSShatyt(fit_marss, only.kem = FALSE)$ytt1
-cbind(m = f_kfas$m[1:n], fitted = ytt1_fit[1:n], MARSShatyt = ytt1_hatyt[1:n])
+cbind(m = kf_kfas$m[1:n], fitted = ytt1_fit[1:n], MARSShatyt = ytt1_hatyt[1:n])
 
 
 ###################################################
@@ -208,7 +208,7 @@ var.Eytt1_fit <-
 var.Eytt1_hatyt <-
   MARSShatyt(fit_marss, only.kem = FALSE)$var.Eytt1
 cbind(
-  P_mu = f_kfas$P_mu[1:n], fitted = var.Eytt1_fit[1:n],
+  P_mu = kf_kfas$P_mu[1:n], fitted = var.Eytt1_fit[1:n],
   MARSShatyt = var.Eytt1_hatyt[1:n]
 )
 
@@ -219,7 +219,7 @@ cbind(
 ytT_fit <- fitted(fit_marss, type = "ytT")$.fitted
 ytT_hatyt <- MARSShatyt(fit_marss)$ytT
 cbind(
-  a = f_kfas$muhat[1:n], fitted = ytT_fit[1:n],
+  a = kf_kfas$muhat[1:n], fitted = ytT_fit[1:n],
   MARSShatyt = ytT_hatyt[1:n], Nile = Nile[1:n]
 )
 
@@ -232,7 +232,7 @@ var.EytT_fit <-
 var.EytT_hatyt <-
   MARSShatyt(fit_marss, only.kem = FALSE)$var.EytT
 cbind(
-  V_mu = f_kfas$V_mu[1:n], fitted = var.EytT_fit[1:n],
+  V_mu = kf_kfas$V_mu[1:n], fitted = var.EytT_fit[1:n],
   MARSShatyt = var.EytT_hatyt[1:n]
 )
 
@@ -371,7 +371,8 @@ head(df)
 ###################################################
 ### code chunk number 50: Cs502_plotting
 ###################################################
-plot(fit_marss, plot.type = "fitted.ytT", pi.int = TRUE)
+plot.type <- ifelse(packageVersion("MARSS") < '3.11.4', "model.ytT", "fitted.ytT")
+plot(fit_marss, plot.type = plot.type, pi.int = TRUE)
 
 
 ###################################################

inst/doc/Chapter_MLR.R

@@ -244,10 +244,8 @@ xyplot(Reaction ~ Days | Subject, sleepstudy,
 # number of subjects
 nsub <- length(unique(sleepstudy$Subject))
 ndays <- length(sleepstudy$Days) / nsub
-# each subject is a row with day across the columns
 dat <- matrix(sleepstudy$Reaction, nsub, ndays, byrow = TRUE)
 rownames(dat) <- paste("sub", unique(sleepstudy$Subject), sep = ".")
-# the day number 0 to 9 is the explanatory variable
 exp.var <- matrix(sleepstudy$Days, 1, ndays, byrow = TRUE)
 
 

inst/doc/Chapter_PVA.R

@@ -93,59 +93,50 @@ params[nsim + 2, ] <- c(sim.u, sim.u, sim.Q, sim.R, sim.Q)
 ###################################################
 # Needs Example 2 to be run first
 par(mfrow = c(3, 3))
-pd <- 0.1
-xd <- -log(pd) # decline threshold
-te <- 100
-tyrs <- 1:te # extinction time horizon
+pd <- 0.1; xd <- -log(pd) # decline threshold
+te <- 100; tyrs <- 1:te # extinction time horizon
 for (j in c(10, 1:8)) {
   real.ex <- denn.ex <- kal.ex <- matrix(nrow = te)
 
   # MARSS parameter estimates
-  u <- params[j, 1]
-  Q <- params[j, 3]
+  u <- params[j, 1]; Q <- params[j, 3]
   if (Q == 0) Q <- 1e-4 # just so the extinction calc doesn't choke
   p.ever <- ifelse(u <= 0, 1, exp(-2 * u * xd / Q))
   for (i in 1:100) {
     if (is.finite(exp(2 * xd * abs(u) / Q))) {
-      sec.part <- exp(2 * xd * abs(u) / Q) * pnorm((-xd - abs(u) * tyrs[i]) / sqrt(Q * tyrs[i]))
-    } else {
-      sec.part <- 0
-    }
+      sec.part <- exp(2 * xd * abs(u) / Q) * 
+        pnorm((-xd - abs(u) * tyrs[i]) / sqrt(Q * tyrs[i]))
+    } else { sec.part <- 0 }
     kal.ex[i] <- p.ever * pnorm((-xd + abs(u) * tyrs[i]) / sqrt(Q * tyrs[i])) + sec.part
   } # end i loop
 
   # Dennis et al 1991 parameter estimates
-  u <- params[j, 2]
-  Q <- params[j, 5]
+  u <- params[j, 2]; Q <- params[j, 5]
   p.ever <- ifelse(u <= 0, 1, exp(-2 * u * xd / Q))
   for (i in 1:100) {
     denn.ex[i] <- p.ever * pnorm((-xd + abs(u) * tyrs[i]) / sqrt(Q * tyrs[i])) +
-      exp(2 * xd * abs(u) / Q) * pnorm((-xd - abs(u) * tyrs[i]) / sqrt(Q * tyrs[i]))
+      exp(2 * xd * abs(u) / Q) * 
+      pnorm((-xd - abs(u) * tyrs[i]) / sqrt(Q * tyrs[i]))
   } # end i loop
 
   # True parameter values
-  u <- sim.u
-  Q <- sim.Q
+  u <- sim.u; Q <- sim.Q
   p.ever <- ifelse(u <= 0, 1, exp(-2 * u * xd / Q))
   for (i in 1:100) {
     real.ex[i] <- p.ever * pnorm((-xd + abs(u) * tyrs[i]) / sqrt(Q * tyrs[i])) +
-      exp(2 * xd * abs(u) / Q) * pnorm((-xd - abs(u) * tyrs[i]) / sqrt(Q * tyrs[i]))
+      exp(2 * xd * abs(u) / Q) * 
+      pnorm((-xd - abs(u) * tyrs[i]) / sqrt(Q * tyrs[i]))
   } # end i loop
 
-  # plot it
-  plot(tyrs, real.ex,
-    xlab = "Time steps into future",
-    ylab = "Probability of extinction", ylim = c(0, 1), bty = "l"
-  )
+  plot(tyrs, real.ex, xlab = "Time steps into future",
+    ylab = "Probability of extinction", ylim = c(0, 1), bty = "l")
   if (j <= 8) title(paste("simulation ", j))
   if (j == 10) title("average over sims")
   lines(tyrs, denn.ex, type = "l", col = "red", lwd = 2, lty = 1)
   lines(tyrs, kal.ex, type = "l", col = "green", lwd = 2, lty = 2)
 }
-legend("bottomright", c("True", "Dennis", "KalmanEM"),
-  pch = c(1, -1, -1),
-  col = c(1, 2, 3), lty = c(-1, 1, 2), lwd = c(-1, 2, 2), bty = "n"
-)
+legend("bottomright", c("True", "Dennis", "KalmanEM"), pch = c(1, -1, -1),
+  col = c(1, 2, 3), lty = c(-1, 1, 2), lwd = c(-1, 2, 2), bty = "n")
 
 
 ###################################################