library(tidyverse)
vec_prob <- c(6, 3, 2, 1.5, 1, 0.5, 0.1)
example_data <- crossing(
id = c(1:100),
day = 1:7
) %>%
group_by(id) %>%
mutate(sleep = sample(c("normal", "slow_wave", "REM"), size = 1)) %>%
mutate(ptsd = sample(c("ptsd", "trauma_exposed"), size = 1)) %>%
ungroup() %>%
rowwise() %>%
mutate(intrusion = sample(x = c(0,1:6),
size = 1,
prob = c(vec_prob / sum(vec_prob)))
) %>%
ungroup()
# I can't quite get the levels even in the sampling
count(example_data, sleep)
#> # A tibble: 3 x 2
#> sleep n
#> <chr> <int>
#> 1 normal 217
#> 2 REM 280
#> 3 slow_wave 203
count(example_data, ptsd)
#> # A tibble: 2 x 2
#> ptsd n
#> <chr> <int>
#> 1 ptsd 406
#> 2 trauma_exposed 294
# the intrusions happen more frequently for 0
ggplot(example_data,
aes(x = intrusion)) +
geom_bar()
ggplot(example_data,
aes(x = day,
y = jitter(intrusion),
group = factor(id))) +
geom_line()
# facetting by the conditions
ggplot(example_data,
aes(x = day,
y = jitter(intrusion),
group = factor(id))) +
geom_line() +
geom_point() +
facet_grid(ptsd~sleep)
library(MASS)
#>
#> Attaching package: 'MASS'
#> The following object is masked from 'package:dplyr':
#>
#> select
library(lme4)
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#>
#> expand, pack, unpack
glm_fit <-
glmer.nb(formula = intrusion ~ sleep*ptsd + (day | id),
data = example_data)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge with max|grad| = 0.0669895 (tol = 0.002, component 1)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see ?isSingular
# glm(formula = intrusion ~ sleep*ptsd + (day | id),
# family = quasipoisson(link = "log"),
summary(glm_fit)
#> Generalized linear mixed model fit by maximum likelihood (Laplace
#> Approximation) [glmerMod]
#> Family: Negative Binomial(1.3062) ( log )
#> Formula: intrusion ~ sleep * ptsd + (day | id)
#> Data: example_data
#>
#> AIC BIC logLik deviance df.resid
#> 2214.5 2260.0 -1097.2 2194.5 690
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -0.8337 -0.7965 -0.1759 0.5085 3.0826
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> id (Intercept) 0.021006 0.14493
#> day 0.002063 0.04542 -1.00
#> Number of obs: 700, groups: id, 100
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.24252 0.11762 2.062 0.0392 *
#> sleepREM -0.05326 0.15051 -0.354 0.7234
#> sleepslow_wave -0.01436 0.16327 -0.088 0.9299
#> ptsdtrauma_exposed 0.08494 0.16842 0.504 0.6140
#> sleepREM:ptsdtrauma_exposed -0.07178 0.22831 -0.314 0.7532
#> sleepslow_wave:ptsdtrauma_exposed 0.06522 0.24269 0.269 0.7881
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Correlation of Fixed Effects:
#> (Intr) slpREM slpsl_ ptsdt_ sREM:_
#> sleepREM -0.734
#> sleepslw_wv -0.706 0.529
#> ptsdtrm_xps -0.666 0.516 0.479
#> slpREM:pts_ 0.480 -0.659 -0.347 -0.736
#> slpslw_wv:_ 0.461 -0.357 -0.666 -0.694 0.511
#> convergence code: 0
#> boundary (singular) fit: see ?isSingular
plot(glm_fit)
Created on 2020-09-24 by the reprex package (v0.3.0)