Nothing
params <-
list(EVAL = FALSE)
## ----setup, include=FALSE, message=FALSE--------------------------------------
library(knitr)
library(glmmTMB)
library(MASS) ## for mvrnorm()
library(TMB) ## for tmbprofile()
library(mvabund) ## for spider data
## devtools::install_github("kaskr/adcomp/TMB") ## get development version
knitr::opts_chunk$set(echo = TRUE, eval=if (exists("params")) params$EVAL else FALSE)
do_image <- exists("params") && params$EVAL
## want to *store* images within package
save_vig_dir <- file.path("inst","vignette_data")
pkg_dir <- "glmmTMB"
## guess where we are ...
if (grepl("/vignettes$",getwd())) { ## in vignettes dir
save_vig_dir <- file.path("../",save_vig_dir)
} else if (grepl(paste0("/",pkg_dir,"$"),getwd())) { ## in repo head
save_vig_dir <- file.path(pkg_dir,save_vig_dir)
}
## want to *retrieve* images from system files
use_vig_dir <- system.file("vignette_data",package="glmmTMB")
mkfig <- function(expr,fn) {
png(normalizePath(file.path(save_vig_dir,fn)))
eval(substitute(expr))
invisible(dev.off())
}
usefig <- function(fn) {
knitr::include_graphics(file.path(use_vig_dir,fn))
}
## turned off caching for now: got error in chunk 'fit.us.2'
## Error in retape() :
## Error when reading the variable: 'thetaf'. Please check data and parameters.
## In addition: Warning message:
## In retape() : Expected object. Got NULL.
set.seed(1)
## run this in interactive session if you actually want to evaluate chunks ...
## Sys.setenv(NOT_CRAN="true")
## ----covstruct-table, echo = FALSE, eval = TRUE------------------------------
ctab <- read.delim(sep = "#", comment = "",
header = TRUE,
check.names = FALSE,
text = "
Covariance # Notation # no. parameters # Requirement # Parameters
Unstructured (general positive definite) # `us` # $n(n+1)/2$ # # See [Mappings]
Heterogeneous Toeplitz # `toep` # $2n-1$ # # log-SDs ($\\theta_1-\\theta_n$); correlations $\\rho_k = \\theta_{n+k}/\\sqrt{1+\\theta_{n+k}^2}$, $k = \\textrm{abs}(i-j+1)$
Het. compound symmetry # `cs` # $n+1$ # # log-SDs ($\\theta_1-\\theta_n$); correlation $\\rho = \\theta_{n+1}/\\sqrt{1+\\theta_{n+1}^2}$
Het. diagonal # `diag` # $n$ # # log-SDs
AR(1) # `ar1` # $2$ # Unit spaced levels # log-SD; $\\rho = \\left(\\theta_2/\\sqrt{1+\\theta_2^2}\\right)^{d_{ij}}$
Ornstein-Uhlenbeck # `ou` # $2$ # Coordinates # log-SD; log-OU rate ($\\rho = \\exp(-\\exp(\\theta_2) d_{ij})$)
Spatial exponential # `exp` # $2$ # Coordinates # log-SD; log-scale ($\\rho = \\exp(-\\exp(-\\theta_2) d_{ij})$)
Spatial Gaussian # `gau` # $2$ # Coordinates # log-SD; log-scale ($\\rho = \\exp(-\\exp(-2\\theta_2) d_{ij}^2$)
Spatial Matèrn # `mat` # $3$ # Coordinates # log-SD, log-range, log-shape (power)
Reduced rank # `rr` # $nd-d(d-1)/2$ # rank (d)
"
)
knitr::kable(ctab)
## ----sim1, eval=TRUE----------------------------------------------------------
n <- 25 ## Number of time points
x <- MASS::mvrnorm(mu = rep(0,n),
Sigma = .7 ^ as.matrix(dist(1:n)) ) ## Simulate the process using the MASS package
y <- x + rnorm(n) ## Add measurement noise
## ----simtimes-----------------------------------------------------------------
# times <- factor(1:n, levels=1:n)
# head(levels(times))
## ----simgroup-----------------------------------------------------------------
# group <- factor(rep(1,n))
## ----simcomb------------------------------------------------------------------
# dat0 <- data.frame(y, times, group)
## ----fitar1, eval=FALSE-------------------------------------------------------
# glmmTMB(y ~ ar1(times + 0 | group), data=dat0)
## ----ar0fit,echo=FALSE--------------------------------------------------------
# glmmTMB(y ~ ar1(times + 0 | group), data=dat0)
## ----simGroup-----------------------------------------------------------------
# simGroup <- function(g, n=6, phi=0.7) {
# x <- MASS::mvrnorm(mu = rep(0,n),
# Sigma = phi ^ as.matrix(dist(1:n)) ) ## Simulate the process
# y <- x + rnorm(n) ## Add measurement noise
# times <- factor(1:n)
# group <- factor(rep(g,n))
# data.frame(y, times, group)
# }
# simGroup(1)
## ----simGroup2----------------------------------------------------------------
# dat1 <- do.call("rbind", lapply(1:1000, simGroup) )
## ----fit.ar1------------------------------------------------------------------
# (fit.ar1 <- glmmTMB(y ~ ar1(times + 0 | group), data=dat1))
## ----fit.us-------------------------------------------------------------------
# fit.us <- glmmTMB(y ~ us(times + 0 | group), data=dat1, dispformula=~0)
# fit.us$sdr$pdHess ## Converged ?
## ----fit.us.vc----------------------------------------------------------------
# VarCorr(fit.us)
## ----fit.toep-----------------------------------------------------------------
# fit.toep <- glmmTMB(y ~ toep(times + 0 | group), data=dat1,
# dispformula=~0)
# fit.toep$sdr$pdHess ## Converged ?
## ----fit.toep.vc--------------------------------------------------------------
# (vc.toep <- VarCorr(fit.toep))
## ----fit.toep.vc.diag---------------------------------------------------------
# vc1 <- vc.toep$cond[[1]] ## first term of var-cov for RE of conditional model
# summary(diag(vc1))
# summary(vc1[row(vc1)!=col(vc1)])
## ----fit.toep.reml------------------------------------------------------------
# fit.toep.reml <- update(fit.toep, REML=TRUE)
# vc1R <- VarCorr(fit.toep.reml)$cond[[1]]
# summary(diag(vc1R))
# summary(vc1R[row(vc1R)!=col(vc1R)])
## ----fit.cs-------------------------------------------------------------------
# fit.cs <- glmmTMB(y ~ cs(times + 0 | group), data=dat1, dispformula=~0)
# fit.cs$sdr$pdHess ## Converged ?
## ----fit.cs.vc----------------------------------------------------------------
# VarCorr(fit.cs)
## ----anova1-------------------------------------------------------------------
# anova(fit.ar1, fit.toep, fit.us)
## ----anova2-------------------------------------------------------------------
# anova(fit.cs, fit.toep)
## ----sample2------------------------------------------------------------------
# x <- sample(1:2, 10, replace=TRUE)
# y <- sample(1:2, 10, replace=TRUE)
## ----numFactor----------------------------------------------------------------
# (pos <- numFactor(x,y))
## ----parseNumLevels-----------------------------------------------------------
# parseNumLevels(levels(pos))
## ----numFactor2---------------------------------------------------------------
# dat1$times <- numFactor(dat1$times)
# levels(dat1$times)
## ----fit.ou-------------------------------------------------------------------
# fit.ou <- glmmTMB(y ~ ou(times + 0 | group), data=dat1)
# fit.ou$sdr$pdHess ## Converged ?
## ----fit.ou.vc----------------------------------------------------------------
# VarCorr(fit.ou)
## ----fit.mat------------------------------------------------------------------
# fit.mat <- glmmTMB(y ~ mat(times + 0 | group), data=dat1, dispformula=~0)
# fit.mat$sdr$pdHess ## Converged ?
## ----fit.mat.vc---------------------------------------------------------------
# VarCorr(fit.mat)
## ----fit.gau------------------------------------------------------------------
# fit.gau <- glmmTMB(y ~ gau(times + 0 | group), data=dat1, dispformula=~0)
# fit.gau$sdr$pdHess ## Converged ?
## ----fit.gau.vc---------------------------------------------------------------
# VarCorr(fit.gau)
## ----fit.exp------------------------------------------------------------------
# fit.exp <- glmmTMB(y ~ exp(times + 0 | group), data=dat1)
# fit.exp$sdr$pdHess ## Converged ?
## ----fit.exp.vc---------------------------------------------------------------
# VarCorr(fit.exp)
## ----spatial_data-------------------------------------------------------------
# d <- data.frame(z = as.vector(volcano),
# x = as.vector(row(volcano)),
# y = as.vector(col(volcano)))
## ----spatial_sub_sample-------------------------------------------------------
# set.seed(1)
# d$z <- d$z + rnorm(length(volcano), sd=15)
# d <- d[sample(nrow(d), 100), ]
## ----volcano_data_image_fake,eval=FALSE---------------------------------------
# volcano.data <- array(NA, dim(volcano))
# volcano.data[cbind(d$x, d$y)] <- d$z
# image(volcano.data, main="Spatial data", useRaster=TRUE)
## ----volcano_data_image_real,echo=FALSE---------------------------------------
# if (do_image) {
# volcano.data <- array(NA, dim(volcano))
# volcano.data[cbind(d$x, d$y)] <- d$z
# mkfig(image(volcano.data, main="Spatial data"),"volcano_data.png")
# }
## ----volcano_image,eval=TRUE,echo=FALSE---------------------------------------
usefig("volcano_data.png")
## ----spatial_add_pos_and_group------------------------------------------------
# d$pos <- numFactor(d$x, d$y)
# d$group <- factor(rep(1, nrow(d)))
## ----fit_spatial_model, cache=TRUE--------------------------------------------
# f <- glmmTMB(z ~ 1 + exp(pos + 0 | group), data=d)
## ----confint_sigma------------------------------------------------------------
# confint(f, "sigma")
## ----newdata_corner-----------------------------------------------------------
# newdata <- data.frame( pos=numFactor(expand.grid(x=1:3,y=1:3)) )
# newdata$group <- factor(rep(1, nrow(newdata)))
# newdata
## ----predict_corner-----------------------------------------------------------
# predict(f, newdata, type="response", allow.new.levels=TRUE)
## ----predict_column-----------------------------------------------------------
# predict_col <- function(i) {
# newdata <- data.frame( pos = numFactor(expand.grid(1:87,i)))
# newdata$group <- factor(rep(1,nrow(newdata)))
# predict(f, newdata=newdata, type="response", allow.new.levels=TRUE)
# }
## ----predict_all--------------------------------------------------------------
# pred <- sapply(1:61, predict_col)
## ----image_results_fake,eval=FALSE--------------------------------------------
# image(pred, main="Reconstruction")
## ----image_results_real,echo=FALSE--------------------------------------------
# if (do_image) {
# mkfig(image(pred, main="Reconstruction", useRaster=TRUE),
# "volcano_results.png")
# }
## ----results_image,eval=TRUE,echo=FALSE---------------------------------------
usefig("volcano_results.png")
## ----fit.us.2-----------------------------------------------------------------
# vv0 <- VarCorr(fit.us)
# vv1 <- vv0$cond$group ## extract 'naked' V-C matrix
# n <- nrow(vv1)
# rpars <- getME(fit.us,"theta") ## extract V-C parameters
# ## first n parameters are log-std devs:
# all.equal(unname(diag(vv1)),exp(rpars[1:n])^2)
# ## now try correlation parameters:
# cpars <- rpars[-(1:n)]
# length(cpars)==n*(n-1)/2 ## the expected number
# cc <- diag(n)
# cc[upper.tri(cc)] <- cpars
# L <- crossprod(cc)
# D <- diag(1/sqrt(diag(L)))
# round(D %*% L %*% D,3)
# round(unname(attr(vv1,"correlation")),3)
## ----other_check--------------------------------------------------------------
# all.equal(c(cov2cor(vv1)),c(fit.us$obj$env$report(fit.us$fit$parfull)$corr[[1]]))
## ----fit.us.profile,cache=TRUE------------------------------------------------
# ## want $par, not $parfull: do NOT include conditional modes/'b' parameters
# ppar <- fit.us$fit$par
# length(ppar)
# range(which(names(ppar)=="theta")) ## the last n*(n+1)/2 parameters
# ## only 1 fixed effect parameter
# tt <- tmbprofile(fit.us$obj,2,trace=FALSE)
## ----fit.us.profile.plot_fake,eval=FALSE--------------------------------------
# confint(tt)
# plot(tt)
## ----fit.us.profile.plot_real,echo=FALSE--------------------------------------
# mkfig(plot(tt),"us_profile_plot.png")
## ----us_profile_image,eval=TRUE,echo=FALSE------------------------------------
usefig("us_profile_plot.png")
## ----fit.cs.profile,cache=TRUE------------------------------------------------
# ppar <- fit.cs$fit$par
# length(ppar)
# range(which(names(ppar)=="theta")) ## the last n*(n+1)/2 parameters
# ## only 1 fixed effect parameter, 1 dispersion parameter
# tt2 <- tmbprofile(fit.cs$obj,3,trace=FALSE)
## ----fit.cs.profile.plot_fake,eval=FALSE--------------------------------------
# plot(tt2)
## ----fit.cs.profile.plot_real,echo=FALSE--------------------------------------
# mkfig(plot(tt2),"cs_profile_plot.png")
## ----fit.cs.profile_image,echo=FALSE,eval=TRUE--------------------------------
usefig("cs_profile_plot.png")
## ----rr_ex, eval = FALSE------------------------------------------------------
# if (require(mvabund)) {
# data(spider)
# ## organize data into long format
# sppTot <- sort(colSums(spider$abund), decreasing = TRUE)
# tmp <- cbind(spider$abund, spider$x)
# tmp$id <- 1:nrow(tmp)
# spiderDat <- reshape(tmp,
# idvar = "id",
# timevar = "Species",
# times = colnames(spider$abund),
# varying = list(colnames(spider$abund)),
# v.names = "abund",
# direction = "long")
# ## fit rank-reduced models with varying dimension
# fit_list <- lapply(2:10,
# function(d) {
# fit.rr <- glmmTMB(abund ~ Species + rr(Species + 0|id, d = d),
# data = spiderDat)
# })
# ## compare fits via AIC
# aic_vec <- sapply(fit_list, AIC)
# aic_vec - min(aic_vec, na.rm = TRUE)
## ----mm_int, eval = TRUE------------------------------------------------------
model.matrix(~f, data.frame(f=factor(c("c", "s", "v"))))
## ----mm_noint, eval = TRUE----------------------------------------------------
model.matrix(~0+f, data.frame(f=factor(c("c", "s", "v"))))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.