Nothing
## ----setup, include=FALSE-----------------------------------------------------
old <- options(digits = 3)
knitr::opts_chunk$set(echo = TRUE)
## -----------------------------------------------------------------------------
library("plyr")
library("mvtnorm")
make_data_generator <- function(resid_var = 1,
ranef_covar = diag(c(1, 1)), n = 100
) {
ni <- nrow(ranef_covar)
generate_data <- function() {
# sample data set under mixed effects model with random slope/intercepts
simulated_data <- rdply(n, {
b <- t(rmvnorm(n = 1, sigma = ranef_covar))
epsilon <- rnorm(n = length(b), mean = 0, sd = sqrt(resid_var))
b + epsilon
})
data.frame(
subject = rep(1:n, each = ni),
age = rep(1:ni, n),
simulated_data)
}
}
## -----------------------------------------------------------------------------
resid_var <- 0
resid_var <- 1
set.seed(77711)
covar <- matrix(c(1, 0.7, 0.5, 0.3,
0.7, 1, 0.8, 0.5,
0.5, 0.8, 1, 0.6,
0.3, 0.5, 0.6, 1), nrow = 4)
gen_dat <- make_data_generator(n = 10000,
ranef_covar = covar,
resid_var = resid_var)
data <- gen_dat()
head(data)
## -----------------------------------------------------------------------------
library("tidyr")
library("dplyr")
d <- as_tibble(data[,-3])
broad <- t(spread(d, subject, X1))[-1,]
cor(broad)
## -----------------------------------------------------------------------------
library("brokenstick")
knots <- 1:3
boundary <- c(1, 4)
fit <- brokenstick(X1 ~ age | subject, data,
knots = knots, boundary = boundary,
method = "lmer")
omega <- get_omega(fit, hide = "no")
beta <- coef(fit, hide = "no")
sigma2 <- fit$sigma2
round(beta, 2)
round(sigma2, 4)
# correlation random effects
round(covar, 3)
round(omega, 2)
# covariances measured data
round(omega + diag(sigma2, 4), 3)
round(cov(broad), 3)
# convert to time-to-time correlation matrix
round(cov2cor(omega + diag(sigma2, 4)), 3)
round(cor(broad), 3)
z <- predict(fit, x = "knots", include_data = FALSE, shape = "wide")[, -1]
# off-diagonal elements of covariance of broken stick estimates approach correlation
# not enough variance in the diagonal because of smoothing
cov(z)
# correlations of broken stick estimates are inflated because of smoothing
cor(z)
## -----------------------------------------------------------------------------
cov <- get_omega(fit)
chat <- cov + diag(fit$sigma2, nrow(cov))
r <- cov2cor(chat)
r
## ----echo = FALSE-------------------------------------------------------------
options(old)
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.