Jak symulować dane, aby wykazać mieszane efekty z R (lme4)?

Jako odpowiednik tego postu pracowałem nad symulacją danych ze zmiennymi ciągłymi, nadając się do skorelowanych przechwyceń i nachyleń.

Chociaż na stronie i poza nią są świetne posty na ten temat , miałem trudności ze znalezieniem od początku do końca przykładu z symulowanymi danymi, który byłby podobny do prostego scenariusza z prawdziwego życia.

Pytanie brzmi więc, jak symulować te dane i „testować” je lmer. Dla wielu nic nowego, ale być może przydatne dla wielu innych, którzy szukają zrozumienia modeli mieszanych.

Jeśli wolisz format artykułu na blogu, Hierarchiczne modele liniowe i lmer to artykuł, który napisałem, który zawiera symulację z losowymi nachyleniami i przechwytywaniami. Oto kod symulacyjny, którego użyłem:

rm(list = ls())
set.seed(2345)

N <- 30
unit.df <- data.frame(unit = c(1:N), a = rnorm(N))

head(unit.df, 3)
unit.df <-  within(unit.df, {
  E.alpha.given.a <-  1 - 0.15 * a
  E.beta.given.a <-  3 + 0.3 * a
})
head(unit.df, 3)

library(mvtnorm)
q = 0.2
r = 0.9
s = 0.5
cov.matrix <- matrix(c(q^2, r * q * s, r * q * s, s^2), nrow = 2,
                     byrow = TRUE)
random.effects <- rmvnorm(N, mean = c(0, 0), sigma = cov.matrix)
unit.df$alpha <- unit.df$E.alpha.given.a + random.effects[, 1]
unit.df$beta <- unit.df$E.beta.given.a + random.effects[, 2]
head(unit.df, 3)

J <- 30
M = J * N  #Total number of observations
x.grid = seq(-4, 4, by = 8/J)[0:30]

within.unit.df <-  data.frame(unit = sort(rep(c(1:N), J)), j = rep(c(1:J),
                              N), x =rep(x.grid, N))
flat.df = merge(unit.df, within.unit.df)

flat.df <-  within(flat.df, y <-  alpha + x * beta + 0.75 * rnorm(n = M))
simple.df <-  flat.df[, c("unit", "a", "x", "y")]
head(simple.df, 3)

library(lme4)
my.lmer <-  lmer(y ~ x + (1 + x | unit), data = simple.df)
cat("AIC =", AIC(my.lmer))
my.lmer <-  lmer(y ~ x + a + x * a + (1 + x | unit), data = simple.df)
summary(my.lmer)

Dane są całkowicie fikcyjne, a kod, którego użyłem do ich wygenerowania, można znaleźć tutaj .

Chodzi o to, że wykonujemy pomiary na glucose concentrationsgrupie 30 athletespo zakończeniu 15 racesw stosunku do stężenia makijażu amino acid A( AAA) we krwi tych sportowców.

Model to: lmer(glucose ~ AAA + (1 + AAA | athletes)

Istnieje stałe nachylenie efektu (stężenie glukozy ~ aminokwasu A); Jednak stoki również różnią się między różnymi sportowców z mean = 0i sd = 0.5, zaś przechwytuje dla poszczególnych zawodników są szerzyć wokół efektów losowych 0z sd = 0.2. Ponadto istnieje korelacja między przechwytywaniem i nachyleniem 0,8 w obrębie tego samego sportowca.

Te losowe efekty są dodawane do wybranych intercept = 1dla ustalonych efektów, oraz slope = 2.

alpha + AAA * beta + 0.75 * rnorm(observations)1 + random effects changes in the intercept$+$AAA $*$2 + random effect changes in slopes for each athlete$+$ noise$\epsilon$sd = 0.75

Dane wyglądają więc tak:

      athletes races      AAA   glucose
    1        1     1  51.79364 104.26708
    2        1     2  49.94477 101.72392
    3        1     3  45.29675  92.49860
    4        1     4  49.42087 100.53029
    5        1     5  45.92516  92.54637
    6        1     6  51.21132 103.97573
    ...

Nierealistyczne poziomy glukozy, ale wciąż ...

Podsumowanie zwraca:

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 athletes (Intercept) 0.006045 0.07775      
          AAA         0.204471 0.45218  1.00
 Residual             0.545651 0.73868      
Number of obs: 450, groups:  athletes, 30

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   1.31146    0.35845 401.90000   3.659 0.000287 ***
AAA           1.93785    0.08286  29.00000  23.386  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

1Zamiast tego występuje korelacja efektów losowych 0.8. sd = 2Na przypadkowej zmienności przechwytuje interpretowany jest jako 0.07775. Standardowe odchylenie 0.5dla losowych zmian stoków wśród sportowców oblicza się jako 0.45218. Hałas ustawiony ze standardowym odchyleniem 0.75został zwrócony jako 0.73868.

Przechwytywanie ustalonych efektów miało być 1, i dostaliśmy 1.31146. Dla zbocza miało to być 2, a oszacowanie było 1.93785.

Dość blisko!