autocorr: calculate an auto-correlation term

In gobbios/cfp: Christof's functions

calculate an auto-correlation term

Description

calculate an auto-correlation term

Usage

autocorr(
  modres,
  xdata,
  corrvar,
  contrfac,
  doplot = TRUE,
  plot_sds = c(0.001, 2),
  opti_sds = c(0.001, 1),
  plot_n = 51,
  stdzcorr = TRUE
)

Arguments

modres	the model result object from (g)lmer
xdata	the data.frame used to fit the model
corrvar	character, the column name in xdata for the variable relevant for the correlation (e.g. date)
contrfac	character, the column name in xdata for the control factor (e.g. ID or group)
doplot	logical, should a plot be produced
plot_sds	numeric of length 2, the SDs for visualization if required
opti_sds	numeric of length 2, the SDs for the optimization
plot_n	numeric, the number of steps for the visualization
stdzcorr	logical, should the correlation-variable corrvar be z-transformed

Details

This function is heavily based on code by Roger Mundry.

The difference between plotSDs and optiSDs is that the former is only used for the visualization, and generally should cover at least equal if not wider range than the latter

Value

a list (and a plot)

Examples

# generate model
set.seed(123)
n = 1000
xdata <- data.frame(pred1 = rnorm(n), pred2 = rnorm(n))
xdata$ID <- factor(sample(LETTERS[1:20], size = n, replace = TRUE))
re <- rnorm(20)
names(re) <- LETTERS[1:20]
xdata$resp <- re[as.numeric(xdata$ID)] + xdata$pred1 * 0.3 + xdata$pred2 * (-0.4) +
              diffinv(rnorm(n, sd = 1.5))[-n] + rnorm(n, sd = 0.5)
## Not run: 
library(lme4)
res <- lmer(resp ~ pred1 + pred2 + (1|ID), data = xdata)
summary(res)
# visual
autocorr(modres = res, xdata = xdata, corrvar = "pred1", contrfac = "ID", plot_sds = c(0.001, 1))
# narrow optimization range and suppress plot
xres <- autocorr(modres = res, xdata = xdata, corrvar = "pred1", contrfac = "ID",
                 opti_sds = c(0.05, 0.25), doplot = FALSE)
# add location of local minimum to previous plot
abline(v = xres$opti$minimum, col = "red")
xdata$ac <- xres$optiac
plot(xdata$ac, xdata$resp)
res2 <- update(res, .~. +ac)
summary(res2)

## End(Not run)

gobbios/cfp documentation built on April 11, 2022, 2:22 a.m.