This is regarding , how to coerce a non tidy object into tidy universe. I am using a relatively new r package GLMMadaptive to fit a mixed model.
This is the code
Step1 : Creating the dataset
set.seed(1234)
n <- 200 # number of subjects
K <- 12 # number of measurements per subject
t_max <- 14 # maximum follow-up time
# construct a dataset with the design:
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.2, -0.25, 0.24, -0.05) # fixed effects coefficients
sigma <- 0.5 # errors' standard deviation
D11 <- 1.0 # variance of random intercepts
D22 <- 0.5 # variance of random slopes
# simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# simulate normal longitudinal data
DF$y <- rnorm(n * K, mean = eta_y, sd = sigma)
# assume that values below -4 are not observed, and set equal to -4
DF$ind <- as.numeric(DF$y < -4)
DF$y <- pmax(DF$y, -4)
Step2 : Run the model
foo <- mixed_model(cbind(y, ind) ~ sex * time,
random = ~ time | id,
data = DF,
family = censored.normal())
Step3 : Model results
summary(foo)
Call:
mixed_model(fixed = cbind(y, ind) ~ sex * time, random = ~time |
id, data = DF, family = censored.normal())
Fit statistics:
log.Lik AIC BIC
-2377.952 4771.904 4798.291
Random effects covariance matrix:
StdDev Corr
(Intercept) 0.9623
time 0.6341 0.1459
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -2.2468 0.1008 -22.2968 < 1e-04
sexfemale -0.0376 0.1445 -0.2599 0.79491
time 0.3231 0.0638 5.0628 < 1e-04
sexfemale:time -0.0913 0.0933 -0.9795 0.32734
log(residual std. dev.):
Estimate Std.Err
-0.6932 0.0175
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE
I am unable to pipe these model summary through tidy and I need help. Thanks in advance.
library(broom.mixed)
tidy(foo)
term estimate std.error statistic p.value conf.low conf.high
1 (Intercept) -2.24677968 0.10076694 -22.2967947 3.969327e-110 -2.4442792 -2.04928012
2 sexfemale -0.03757411 0.14454807 -0.2599419 7.949086e-01 -0.3208831 0.24573491
3 time 0.32312932 0.06382416 5.0628059 4.131305e-07 0.1980363 0.44822238
4 sexfemale:time -0.09134472 0.09325713 -0.9794932 3.273364e-01 -0.2741253 0.09143589
Doesn't appear the random effects extraction is implemented in broom.mixed for this model yet.