Mixed Models
Mixed model 은 오랜동안 통계학 분야에서 사용되어 왔다. 예를들어, 표준 ANOVA 방법은 mixed model 의 특별한 케이스라고 볼 수 있다. 최근에는, mixed model 은 다양하게 응용되었고 확장되어, 다양한 범위의 데이터 상황들에 적용할 수 있게 되었다. 일반화 선형모형 이상의 툴셋을 확장하는데에 첫 걸음으로 볼 수 있다.
용어
일반인들에게, mixed model 을 둘러싼 용어들은 조금 헷갈릴 수 있다. 이러한 종류의 모형과 연관된 용어들은 다음과 같다:
- Variance components
- Random intercepts and slopes
- Random effects
- Random coefficients
- Varying coefficients
- Intercepts- and/or slopes-as-outcomes
- Hierarchical linear models
- Multilevel models (implies multiple levels of hierarchically clustered data)
- Growth curve models (possibly Latent GCM)
- Mixed effects models
이 모두는 mixed model 의 다른 종류를 말한다. 어떤 용어들은 역사가 깊고, 어떤 것들은 특수 분야에서 자주 사용되고, 어떤 것들은 특정 데이터 구조를 뜻하고, 어떤 것들은 특수한 케이스들이다. Mixed effects 혹은 mixed model 은 일반적으로 fixed 와 random effect들이 함께 있는 것을 의미한다. 일반적으로 나는 ‘mixed models’ 혹은 ’random effects models’를 선호하는데, 더 간단한 용어들이고, 특정한 구조를 의미하지 않고, 후자는 다른 용어들이 사용될 때 많은 사람들이 생각하지 않았을 확장적용될 수 있다1. Mixed effects 에 관해서는, fixed effects 는 linear regression model 의 main effect 로서는 약하지만 고루한 용어, 즉, mixed model 에서 random 이 아닌 파트이고, 어떤 문맥에서 population average effect 를 의미한다. Random effect 에 대해서 많은 정의를 듣게 되겠지만, 어떻게 정의되던지 간에 한 관측단위에 특정한 것이다. 이 문서에서는 관측단위가 어떤 그룹화 factor 의 레벨인 경우를 다루지만, 이 경우는 단지 한 예일 뿐이다.
Kinds of Clustering
데이터에 하나 이상의 군집이 있을 수 있고 이 클러스터들은 hierachical 해서 다른 클러스터에 포함되어(nested) 있을 수 있다. 학생들이 여러번 응시한 수학능력시험(SAT) 가 그런 예이다 (학생들 내에 반복된 관측값들이 포함되어 있고, 학생들은 학교 안에 포함되었고, 학교는 학군 내에 포함되어 있다). 포함하지 않는 구조도 있다. 참가자들이 같은 업무를 수행하는 반응시간 실험이 예이다. 관측값들은 개인 내에 포함되어 있지만, 관측값들은 또한 업무유형에 따라 클러스터를 이룬다. 이러한 시나리오를 구분짓기 위해 겹침(nested) 과 crossed 라는 용어가 사용되기도 한다. 또한, 클러스터는 균형을 이룰 수도 있고 아닐 수도 있다. 실험에서는 균형 클러스터가 많겠지만, 아닌 경우도 있다. 예를 들어, 지리적 단위가 클러스터이거나 관측값이 사람인 경우이다.
아래에서 우리는 이러한 데이터 상황들에서 mixed effect model 을 살펴볼 것이다. 일반적으로 우리는 클러스터링이 모델보다는 데이터의 특성이라는 접근법을 취할 것이다. 그러나, 중요한 점은 mixed model 은 유연해서 다양한 데이터 상황을 다룰 수 있다는 점이다.
Random Intercepts Model
지금부터 mixed model 중 가장 간단하고2 가장 자주 일어나는 케이스인 random effect 의 그룹화/군집 구조가 하나인 케이스를 살펴볼 것이다. 이를 random intercepts model 이라고 부르는데, 이유는 곧 알게 될 것이다.
예: 학점 (GPA)
지금부터 우리는 대학 학점을 예측하는 요소를 조사할 것이다. 200 명의 학생을 여섯 번 평가했다 (처음 3년간 모든 학기).
관측값들은 학생 내에서 클러스터되어 있다.
취업상황, 성별, 고등학교 학점과 같은 변수들도 있다. 이들 중 일부는 라벨되어 있고 수치형 형태로 되어 있다.
자세한 내용은 appendix 을 참고하라.
일반적인 회귀 모형
기본모형을 여러 방법으로 표현할 수 있다. 첫번째로는 표준 회귀부터 시작해서 방향을 잡아본다.
\[\mathscr{gpa} = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathscr{occasion} + \epsilon\]
절편과 시간효과를 위한 계수들 (\(b\)) 이 있다. 오차 (\(\epsilon\)) 는 평균 0 이고 어떤 편차 \(\sigma\) 를 가지는 정규분포를 이룬다고 가정한다 .
\[\epsilon \sim \mathscr{N}(0, \sigma)\] 또 다른 방법인, \(\mathrm{gpa}\) 데이터 생성과정을 강조하는 모델은 다음과 같이 쓸 수 있다.
\[\mathscr{gpa} \sim \mathscr{N}(\mu, \sigma)\] \[\mu = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathscr{occasion}\] GPA 와 \(\mu\) 변수는 각 관측값을 나타내는 첨자를 생략했지만, 한 시점에 한 개인에 관한 모델이라고 생각할 수 있다.
The Mixed Model
첫모델
이제 우리는 각 학생들마다 고유한 효과를 포함하는 mixed model 로 설명하는 방법을 살펴볼 것이다. 한 학생에 대해 다음 모형을 살펴보자.3 이는 학생에 특수한 효과, 즉, 주어진 학생이 GPA 에서의 변인은 분산 원인으로 볼 수 있다.
\[\mathscr{gpa} = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathscr{occasion} + (\mathrm{effect}_{\mathscr{student}} + \epsilon)\] 학생 효과에 대해 (관용적으로) 다음을 가정할 수 있다.
\[\mathrm{effect}_{\mathrm{student}} \sim \mathscr{N}(0, \tau)\]
따라서 학생효과는 랜덤, 구체적으로는 평균 0 과 측정된 표준편차 (\(\tau\)) 를 갖는 정규분포를 이룬다. 다른 말로 하면, mixed model 과 표준 회귀 사이의 개념적 차이는 학생효과 뿐인데, 평균으로 보면 효과가 없지만, 학생마다 평균 \(\tau\) 만큼 차이가 난다.
다시 정리하면, 오차의 추가 소스보다 모델 계수들이 강조된다.
\[\mathscr{gpa} = (b_{\mathrm{intercept}} + \mathrm{effect}_{\mathscr{student}}) + b_{\mathrm{occ}}\cdot \mathscr{occasion} + \epsilon\] 더 간단하게 하면:
\[\mathscr{gpa} = b_{\mathrm{int\_student}} + b_{\mathrm{occ}}\cdot \mathscr{occasion} + \epsilon\]
우리는 학생고유의 절편을 갖게 되는데, 구체적으로는 학생 각자가 전체 절편에 고유 효과가 더해지게 되어, 각 사람마다 다른 절편을 갖게 된다.
\[b_{\mathrm{int\_student}} \sim \mathscr{N}(b_{\mathrm{intercept}}, \tau)\] 이제 절편이 평균이 전체(overall) 절편인 정규분포를 이루는 것을 보았다. 이를 종종 랜덤 절편(random intercepts) 모델이라고 부른다.
Multi-level model
mixed model 을 표현하는 다른 방법은 multilevel 모델링 논문들에서 볼 수 있다. 두 개의 파트 회귀 모형으로 표현하면 더 명시적인데, 한 파트는 관측값 레벨과 다른 파트는 학생 레벨이다.
\[\mathrm{gpa} = b_{\mathrm{int\_student}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\]
\[b_{\mathrm{int\_student}} = b_{\mathrm{intercept}} + \mathrm{effect}_{\mathrm{student}}\] 하지만, 두번째 레벨부분을 첫번째에 ‘치환하면’ 이전 모델과 동일하게 된다.
우리는 학기마다 학생-특수한 효과를 갖지 않는다는 것을 주목하라. 이러한 배경에서, 학기는 fixed effect 만 갖고, 랜덤 파트는 없다고 말한다. 항상 이래야만 하는 것은 아닌데, 나중에 볼 것이다.
적용
시각화
미리 살펴보는 것은 늘 도움이 된다. GPA vs. 학기 플롯을 그려서 출발점과 변화에 있어서 차이를 간략히 보자.
모든 학생 path들은 연한 선으로 표현되었고 10개 샘플은 진한색으로 표현되어있다. 나중에 볼 회귀모형으로 측정한 전체 트렌드는 빨간색으로 나타냈다. 두 가지가 눈에 띈다. 하나는 처음에 학생들 사이의 변동이 크다는 점이다. 두번째는 일반적으로 GPA 트렌드가 우상향하지만, 개별 학생마다 궤적에 차이를 보인다는 점이다.
표준 회귀
시작해보자. 우선, 회귀적합 결과를 살펴보는데, 시간 indicator 를 covariate 로 두는데 수치형으로 할 것이다. 여기서 요약객체들 중 보기 편한 버전을 사용한 결과이다.
load('data/gpa.RData')
gpa_lm = lm(gpa ~ occasion, data = gpa)
## summary(gpa_lm)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.599 | 0.018 | 145.7 | 0 |
| occasion | 0.106 | 0.006 | 18.04 | 0 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1200 | 0.3487 | 0.2136 | 0.2129 |
이는 절편이 의미하는 시작시점에서, 즉 학기가 0 일 때 평균 학점이 2.6 임을 보여준다. 또한, 학기가 진행됨에 따라, GPA 가 약 0.11 점 만큼 증가함을 알 수 있다. 우리가 군집을 무시한다는 점을 제외하고는 큰 문제가 생기지 않는다. 이렇게 하는 것의 단점은 표준오차가 편향되고 따라서 이를 기반으로 한 통계적 유의성에 관한 주장은 잘못되기 쉽다는 점이다. 하지만 더 중요한 점은 바로 우리의 관심사인 학생 효과를 탐색할 수 없다는 점이다.
Regression by cluster
다른 방법은 각 학생마다 독립된 회귀를 실행하는 것이다. 하지만, 여기에는 많은 단점들이 있다 - 그룹이 많고, 각 군집 내에 요약할 만한 데이터가 매우 적을 때 (이 경우에 해당한다) 쉽게 요약되지 않고, 그리고 모델들이 over-contextualized 된다는 점인데, 이는 학생들이 공통으로 가지고 있는 것을 무시함을 의미한다. 우리는 뒤에서 이러한 방법을 mixed model 과 비교할 것이다.
Running a mixed model
다음으로 학생 효과가 있는 mixed model 을 수행한다. 이 모델은 R 에서 lme4 패키지로 쉽게 수행할 수 있다. 다음의 코드는 lme4 으로 회귀에 사용한 것과 똑같이 보일 수 있지만, 그룹을 특정하는 추가요소, 즉 학생 효과가 있다. (1|student) 는 1 로 표시하는 절편이 학생마다 변하도록 허용한다는 것을 의미한다. 이러한 mixed model 은 회귀와 같은 결과를 얻겠지만, 이에 대해서는 뒤에서 이야기할 것이다.
library(lme4)
gpa_mixed = lmer(gpa ~ occasion + (1 | student), data = gpa)
## summary(gpa_mixed)| term | value | se | t | p_value | lower_2.5 | upper_97.5 |
|---|---|---|---|---|---|---|
| Intercept | 2.599 | 0.022 | 119.800 | 0 | 2.557 | 2.642 |
| occasion | 0.106 | 0.004 | 26.096 | 0 | 0.098 | 0.114 |
| group | effect | variance | sd |
|---|---|---|---|
| student | Intercept | 0.064 | 0.252 |
| Residual | 0.058 | 0.241 |
절편과 학기에 관한 계수들 즉 우리 문맥에서 fixed 효과들은, 표준 회귀에서 보았던 것과 같고4, 해석도 그렇다는 것을 확인할 수 있다. 반면 표준오차는 다르다. 통계적 유의성이 같은 한, 결론은 같겠지만 말이다. 절편의 표준오차가 증가했다는 것에 주목하라. 개념적으로 사람마다 랜덤 절편을 허용해서 개인에 대한 정보를 얻게 해주면서도 전에 과소추정한5, 전체 평균에 관한 불확실성을 알 수 있게 된다.
lme4 가 계수와 표준오차는 출력했지만, p-value 를 주지 않는다는 것을 알아차렸는가? 여기에는 몇몇 이유가 있는데, mixed model 로 다루는 것이 첫번째로는 클러스터마다 샘플 크기가 다른 것들 \(N_c\) 인데, 하나의 관측값이 있는 경우도 있다. 두번째로는 N 개의 총 관측값을 다루는데, 이는 레퍼런스 분포, 분모의 자유도와 최상의 해법을 근사하는 방법에 있어서 상황이 이상해지게 된다. 자동으로 p-value 를 제공하는 프로그램들은, p-value 를 계산하는 여러 방법 중 어떤 방법을 사용했는지 말해주지 않는다. 더군다나, 이러한 근사법은 어떤 시나리오에서는 매우 안 좋거나, 그 상황에 적절하지 않은 가정을 한다6.
하지만, 신뢰 구간을 구하는 것은 더 직관적인데, 이는 다음과 같이7 lme4 로 구할 수 있다.
confint(gpa_mixed)| group | effect | variance | sd | sd_2.5 | sd_97.5 | var_prop | |
|---|---|---|---|---|---|---|---|
| sd_(Intercept)|student | student | Intercept | 0.064 | 0.252 | 0.225 | 0.282 | 0.523 |
| sigma | Residual | 0.058 | 0.241 | 0.231 | 0.252 | 0.477 |
Variance components
표준 회귀의 결과와 비교하여 새로운 점은 학생효과의 표준편차/분산 추정값 (앞의 공식에서 \(\tau/\tau^2\)) 이다. 이는 GPA 가 학생마다 평균적으로 얼마나 변하는지를 알려준다. 다른말로 하면, 시점에 기반하여 예측을 한 뒤라도, 각 학생들은 고유한 변동을 가지고 있는데, (표준 편차라는 이름을 가지고 있는) 이 값은 학생들 사이의 평균 편차의 추정값이다. 학생에 기반하여 점수가 바뀌는 것은 학기가 바뀌는 것에 기반한 변화량의 두 배 이상인 것을 주목하라.
분산 출력값을 해석하는 다른 방법은 전체 분산에 비례한 학생 분산의 백분율인데 0.064 / 0.122 = 52% 와 같다. 이를 intraclass correlation 라고 부르는데, within 클러스터 상관계수의 추정값이며, 뒤에서 더 살펴보겠다.
랜덤효과의 추정값
모델을 실행한 후, 학생 효과8의 값을 구할 수 있다. 앞에서부터 다섯 명의 학생에 대해 두가지 방법에 대해 보여주는데, 하나는 랜덤효과로 다른 하나는 랜덤 절편이다. (즉, 절편 + 랜덤 효과)
ranef(gpa_mixed)$student %>% head(5)
# showing mixedup::extract_random_effects(gpa_mixed)| group_var | effect | group | value | se | lower_2.5 | upper_97.5 |
|---|---|---|---|---|---|---|
| student | Intercept | 1 | -0.071 | 0.092 | -0.251 | 0.109 |
| student | Intercept | 2 | -0.216 | 0.092 | -0.395 | -0.036 |
| student | Intercept | 3 | 0.088 | 0.092 | -0.091 | 0.268 |
| student | Intercept | 4 | -0.187 | 0.092 | -0.366 | -0.007 |
| student | Intercept | 5 | 0.030 | 0.092 | -0.149 | 0.210 |
coef(gpa_mixed)$student %>% head(5)| group_var | effect | group | value | se | lower_2.5 | upper_97.5 |
|---|---|---|---|---|---|---|
| student | Intercept | 1 | 2.528 | 0.095 | 2.343 | 2.713 |
| student | Intercept | 2 | 2.383 | 0.095 | 2.198 | 2.568 |
| student | Intercept | 3 | 2.687 | 0.095 | 2.502 | 2.872 |
| student | Intercept | 4 | 2.412 | 0.095 | 2.227 | 2.597 |
| student | Intercept | 5 | 2.629 | 0.095 | 2.444 | 2.814 |
학기가 모든 학생들마다 변하지 않고 일정하다고, 즉, fixed effect 라고 가정했었다.
우리는 이러한 효과들에 매우 관심이 있고, 이에 관한 불확실성에 관해 알고 싶다. lme4 에서 bootstrapping 을 통해 할 수 있는데, 구체적으로는 lme4 안에 있는 bootMer 함수로 수행할 수 있다. 일부 독자에게는 약간 복잡한 작업일 수 있다. merTools 패키지의 predictInterval 함수9는 쉬운 방법을 제공한다. 바로 이들이 제공하는 플롯을 보자.
library(merTools)
predictInterval(gpa_mixed) # for various model predictions, possibly with new data
REsim(gpa_mixed) # mean, median and sd of the random effect estimates
plotREsim(REsim(gpa_mixed)) # plot the interval estimates다음의 플롯은 각 학생마다의 랜덤 효과의 추정값과 구간 추정값이다 (코드10의 마지막 라인이 만든 플롯의 수정버전). 랜덤 효과는 수평선으로 나타낸, 평균 0 인 정규분표를 이루고 있음을 기억하라. 0을 포함하지 않는 구간은 굵게 나타내었다.
Prediction
표준 예측과 클러스터-기반 예측을 비교해보자. R에 있는 대부분의 모델들에서와 같이, 모델객체에 predict 함수를 사용할 수 있다.
predict(gpa_mixed, re.form = NA) %>% head() 1 2 3 4 5 6
2.599214 2.705529 2.811843 2.918157 3.024471 3.130786 위 코드에서 우리는 랜덤 효과를 사용하지 않는다고 re.form = NA 설정했었고, 따라서 관측값들에 대한 우리의 예측값들은 표준 선형 모형에서 얻는 것과 거의 비슷하다.
predict_no_re = predict(gpa_mixed, re.form = NA)
predict_lm = predict(gpa_lm)하지만 각 사람들은 고유의 절편을 갖기 때문에, 이 정보를 포함했을 때 예측값들이 얼마나 달라지는지 살펴보자.
predict_with_re = predict(gpa_mixed)주어진 학생에 대해, 학생 효과의 추정값에 따라, 그 학생은, 모든 학생들의 절편 추정값보다 크거나 작은 값에서 시작할 것이다. 다음은 앞에서 부터 두 명의 학생에 대한 비조건부(unconditional) 예측 vs. 랜덤 절편을 포함한 조건부 예측의 시각화이다.
우리는 mixed model 의 예측값은 다른 절편을 갖기 때문에 옮겨졌음을 볼 수 있다. 이러한 변동의 의미는 이 두 학생이 성적이 상대적으로 좋지 않게 시작했다는 것이다.
Cluster Level Covariates
Mixed model 을 multilevel 모델로서의 해석을 살펴보자.
\[\mathrm{gpa} = b_{\mathrm{int\_student}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\]
\[b_{\mathrm{int\_student}} = b_{\mathrm{intercept}} + \mathrm{effect}_{\mathrm{student}}\] 예를 들어 성별과 같은 학생레벨의 covariate 을 모델에 추가한다면, 다음과 같이 된다.
\[b_{\mathrm{int\_student}} = b_{\mathrm{intercept}} + b_{sex}\cdot \mathrm{sex} + \mathrm{effect}_{\mathrm{student}}\] 집어 넣은 다음 우리는 전과 같은 모델에 설명변수들이 추가된 형태가 된다.
\[\mathrm{gpa} = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathrm{occasion}+ b_{sex}\cdot \mathrm{sex} + (\mathrm{effect}_{\mathscr{student}} + \epsilon)\] 따라서, 클러스터 레벨 covariate 을 추가한다고 해서, 모델11에 이상한 효과가 더해지지 않는다. 이것들을 설명 변수들에 포함시키기만 하면 된다. 관측값 레벨 변수들의 평균이나 다른 요약값으로 클러스터 레벨 covariate 을 만들 수도 있다. 클러스터가 지리적 단위를 대표하고 관측값들이 사람들인 경우 이러한 상황이 매우 자주 생긴다. 예를 들어, 소득을 개인레벨 covairate 으로 갖고, 그 지역의 평균 경제적 수준을 대표하기 위해 중간값을 사용할 수 있다.
Mixed Model 기초사항의 요약
Mixed model 을 이용하면 데이터 안의 군집을 설명할 수 있다. 이러한 이유가 전부라고 하더라도, 데이터의 구조를 무시했을 때 보다 상대적으로 더 정확한 추론을 할 수 있을 것이다. 하지만 그 이상인 것들이 있다. 우리는 목적 변수의 변동 원인을 더 잘 이해한다. 모델의 파라미터의 그룹고유의 추정값을 얻을 수 있는데, 그룹들이 서로 얼마만큼 정확히 다른지를 이해할 수 있다. 더군다나, 그 결과로 그룹 고유한 예측을 하게 해주어, 그룹 때문에 발생하는 변동이 상당하다는 전제하에 훨씬 더 정확한 예측을 할 수 있다. 요약하면, 가장 간단한 상황에서조차 mixed model 로 훨씬 많은 것을 얻을 수 있다.
연습문제
수면
이 예시에서 lme4 패키지의 수면연구 (sleep study) 데이터를 사용할 것이다. 다음에서 기술되어 있다.
The average reaction time per day for subjects in a sleep deprivation study. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. The observations represent the average reaction time (in milliseconds) on a series of tests given each day to each subject.
패키지를 로딩한 뒤 데이터를 다음과 같이 불러 온다. 첫 몇 관측값은 다음과 같다.
library(lme4)
data("sleepstudy")| Reaction | Days | Subject |
|---|---|---|
| 249.5600 | 0 | 308 |
| 258.7047 | 1 | 308 |
| 250.8006 | 2 | 308 |
| 321.4398 | 3 | 308 |
| 356.8519 | 4 | 308 |
| 414.6901 | 5 | 308 |
Run a regression with Reaction as the target variable and Days as the predictor.
Run a mixed model with a random intercept for Subject.
Interpret the variance components and fixed effects.
클러스터 레벨 공변량(covariate) 추가하기
클러스터 레벨 공변량, sex, or high school GPA (highgpa) 를 각각 혹은 둘 다 추가하여 GPA data 로 mixed model 을 다시 실행하라. 결과의 모든 측면에 대해 해석하라.
모델에 클러스터 레벨 공변량을 추가한 뒤 학생 분산에 어떤 일이 일어났는가?
mixed model 시뮬레이트하기
The following represents a simple way to simulate a random intercepts model. Note each object what each object is, and make sure the code make sense to you. Then run it.
set.seed(1234) # this will allow you to exactly duplicate your result
Ngroups = 100
NperGroup = 3
N = Ngroups * NperGroup
groups = factor(rep(1:Ngroups, each = NperGroup))
u = rnorm(Ngroups, sd = .5)
e = rnorm(N, sd = .25)
x = rnorm(N)
y = 2 + .5 * x + u[groups] + e
d = data.frame(x, y, groups)Which of the above represent the fixed and random effects? Now run the following.
model = lmer(y ~ x + (1|groups), data=d)
summary(model)
confint(model)
library(ggplot2)
ggplot(aes(x, y), data=d) +
geom_point()Do the results seem in keeping with what you expect?
In what follows we’ll change various aspects of the data, then rerun the model after each change, then summarize and get confidence intervals as before. For each note specifically at least one thing that changed in the results.
- First calculate or simply eyeball the intraclass correlation coefficient \(\frac{\textrm{random effect variance}}{\textrm{residual + random effect variance}}\). In addition, create a density plot of the random effects as follows.
re = ranef(model)$groups
qplot(x = re, geom = 'density', xlim = c(-3, 3))- Change the random effect variance/sd and/or the residual variance/sd and note your new estimate of the ICC, and plot the random effect as before.
- Reset the values to the original. Change Ngroups to 50. What differences do you see in the confidence interval estimates?
- Set the Ngroups back to 100. Now change NperGroup to 10, and note again the how the CI is different from the base condition.
I actually like Richly Parameterized Linear Models, or Structured Additive Regression Models. Both are a mouthful, but at least the latter reduces to STARs.↩︎
Actually, the simplest model would have no covariates at all, just variance components, with no correlations among the random effects. Such a model can be interesting to look at while exploring your data, but would probably never suffice on its own to tell the story you desire to.↩︎
Note that I leave out the observation level subscript to keep things clean. I find that multilevel style notation quickly becomes unwieldy, and don’t wish to reproduce it. It also tends to add confusion to a lot of applied researchers starting out with mixed models.↩︎
This will not always be the case, e.g. with unbalanced data, but they should be fairly close.↩︎
The standard error for our time covariate went down due to our estimate of \(\sigma\) being lower for this model, and there being no additional variance due to cluster membership.↩︎
Note that many common modeling situations involve a fuzzy p setting, but especially penalized regression approaches such as mixed, additive, ridge regression models etc. Rather than be a bad thing, this usually is a sign you’re doing something interesting, or handling complexity in an appropriate way.↩︎
See
?confint.merModfor details and options. The output you see is based on my wrappermixedup::extract_vc.↩︎These are sometimes referred to as BLUPs or EBLUPs, which stands for (empirical) best linear unbiased prediction. However, they are only BLUP for linear mixed effects models. As such you will also see them referred to as conditional mode. Furthermore, in the Bayesian context, the effects are actually estimated as additional model parameters, rather than estimated/predicted after the fact.↩︎
Note that while predictionInterval does not quite incorporate all sources of uncertainty as does bootMer, it’s actually feasible for larger data sets, and on par with the Bayesian results (e.g. with rstanarm).↩︎
Note that the default plot from merTools is confusingly labeled for single random effect, because it unnecessarily adds a facet. You’ll understand it better by looking the plot in the discussion of crossed random effects later. However, the one displayed is from my own package, visibly.↩︎
This is why the multilevel depiction is sub-par, and leads many to confusion at times. You have a target variable and predictor variables based on theory. Whether they are cluster level variables or if there are interactions doesn’t have anything to do with the data structure as much as it does the theoretical motivations. However, if you choose to depict the model in multilevel fashion, the final model must adhere to the ‘plugged in’ result. So if, for example, you posit a cluster level variable for a random slope, you must include the implied interaction of the cluster level and observation level covariates.↩︎