線性混合效應模型R語言應用+1

2022-01-31 一個小朋友w

線性混合效應模型R語言應用+1

數據的層次結構模型出發點是個體受環境的影響，從而在同一個環境中的個體，具有一定的相似性。其他類型的層次結構例如在縱向研究中的個體也一樣，每個個體為第二層層，個體的每次觀測為第一層。

層次結構最低的一級為Level 1，較高的一級為Level 2。例如，調查各個學校的學生，學生為Level 1，學校為Level 2。對這樣的數據建模通常被稱為多層次模型或者混合效應模型。

混合效應模型可以看作為具有層次結構的回歸模型，其中Level 1的參數是Level 2的方程。Level 1看起來和OLS很像，但是下標j表明其估計會隨著Level 2值的變化而不同，即研究對象的每個層級單位都有各自的平均水平

Level 2方程顯示了Level 1的參數如何受Level 2變量的影響：

第二個方程的解釋與此類似，但此方程擬合的是Level 2變量對

聯立後的形式如下：

第一部分為模型的固定效應部分（，第二部分為模型的隨機效應部分（u和r），其中

使用層次模型時需要注意的是：

3層以上的模型並非不可以，但社會科學中常用到的層次模型為2~3層。使用隨機截距，還是隨機截距和隨機斜率共同作為層級單位的自變量與因變量的關係的反應，應綜合考慮理論和數據作出決定。用到的包
# load libraries
library(nlme)
library(lme4)
library(lattice)
library(ggplot2)
library(dplyr)
我們使用nlme包的一份重複測量的數據集Orthodont，該數據集有4個變量
data1 <- as.data.frame(nlme::Orthodont)
data1$Subject <- factor(data1$Subject, ordered = FALSE)
head(data1)
str(data1)
'data.frame': 108 obs. of  4 variables:
 $ distance: num  26 25 29 31 21.5 22.5 23 26.5 23 22.5 ...
 $ age     : num  8 10 12 14 8 10 12 14 8 10 ...
 $ Subject : Factor w/ 27 levels "M16","M05","M02",..: 15 15 15 15 3 3 3 3 7 7 ...
 $ Sex     : Factor w/ 2 levels "Male","Female": 1 1 1 1 1 1 1 1 1 1 ...
每個個體age與distance的關係
data1 %>% 
  ggplot(aes(x=age, y=distance, 
             group=Subject, color=Subject, linetype=Sex)) + 
  geom_line(size=1) + theme_classic()

零模型
mod1 <- lmer(distance ~ (1|Subject), data=data1, REML=F)
summary(mod1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: distance ~ (1 | Subject)
##    Data: data1
## 
##      AIC      BIC   logLik deviance df.resid 
##    521.5    529.5   -257.7    515.5      105 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2391 -0.5248 -0.1103  0.4827  2.7734 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 3.567    1.889   
##  Residual             4.930    2.220   
## Number of obs: 108, groups:  Subject, 27
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  24.0231     0.4216   56.98
# 檢驗隨機截距是否顯著
rand(mod1)
ANOVA-like table for random-effects: Single term deletions

Model:
distance ~ (1 | Subject)
              npar  logLik    AIC   LRT Df Pr(>Chisq)    
<none>           3 -257.75 521.49                        
(1 | Subject)    2 -268.79 541.58 22.09  1  2.601e-06 ***
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1
可視化模型結果：
data1 %>% 
    # save predicted values
    mutate(pred_dist = fitted(mod1)) %>% 
    # graph
    ggplot(aes(x=age, y=pred_dist, group=Subject, color=Subject)) + theme_classic() +
    geom_line(size=1) 
隨機截距+固定斜率模型
在零模型上加一個一級的自變量age
mod2 <- lmer(distance ~ age + (1|Subject), data=data1, REML=F)
summary(mod2)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: distance ~ age + (1 | Subject)
   Data: data1

     AIC      BIC   logLik deviance df.resid 
   451.4    462.1   -221.7    443.4      104 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6870 -0.5386 -0.0123  0.4910  3.7470 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 4.294    2.072   
 Residual             2.024    1.423   
Number of obs: 108, groups:  Subject, 27

Fixed effects:
            Estimate Std. Error t value
(Intercept) 16.76111    0.79456   21.09
age          0.66019    0.06122   10.78

Correlation of Fixed Effects:
    (Intr)
age -0.848
可視化模型結果：
data1 %>% 
    # save predicted values
    mutate(pred_dist = fitted(mod2)) %>% 
    # graph
    ggplot(aes(x=age, y=pred_dist, group=Subject, color=Subject)) + theme_classic() +
    geom_line(size=1) 
隨機截距和隨機斜率模型
允許年齡在每個人的身上對於因變量的關係是變化的
mod3 <- lmer(distance ~ age + (1|Subject) + (0+age|Subject), data=data1, REML=F)
summary(mod3)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: distance ~ age + (1 | Subject) + (0 + age | Subject)
   Data: data1

     AIC      BIC   logLik deviance df.resid 
   449.7    463.1   -219.9    439.7      103 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7542 -0.5056  0.0181  0.5216  3.8017 

Random effects:
 Groups    Name        Variance Std.Dev.
 Subject   (Intercept) 1.82570  1.3512  
 Subject.1 age         0.02141  0.1463  
 Residual              1.85944  1.3636  
Number of obs: 108, groups:  Subject, 27

Fixed effects:
            Estimate Std. Error t value
(Intercept) 16.76111    0.70816   23.67
age          0.66019    0.06509   10.14

Correlation of Fixed Effects:
    (Intr)
age -0.822
可視化模型結果
data1 %>% 
    # save predicted values
    mutate(pred_dist = fitted(mod3)) %>% 
    # graph
    ggplot(aes(x=age, y=pred_dist, group=Subject, color=Subject)) + theme_classic() +
    geom_line(size=1) 
隨機截距和隨機斜率（允許兩者相關）
允許隨機截距和隨機斜率相關
mod4 <- lmer(distance ~ age + (1+age|Subject), data=data1, REML=F)
summary(mod4)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: distance ~ age + (1 + age | Subject)
##    Data: data1
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.2    467.3   -219.6    439.2      102 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3060 -0.4874  0.0076  0.4822  3.9228 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Subject  (Intercept) 4.81409  2.1941        
##           age         0.04619  0.2149   -0.58
##  Residual             1.71620  1.3100        
## Number of obs: 108, groups:  Subject, 27
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 16.76111    0.76075  22.032
## age          0.66019    0.06992   9.442
## 
## Correlation of Fixed Effects:
##     (Intr)
## age -0.848
可視化模型結果
data1 %>% 
    # save predicted values
    mutate(pred_dist = fitted(mod4)) %>% 
    # graph
    ggplot(aes(x=age, y=pred_dist, group=Subject, color=Subject)) + theme_classic() +
    geom_line(size=1) 
模型診斷
殘差診斷
# 殘差與擬合值散點圖  檢查非線性關係、離群值、異方差
data1 %>% 
  # data
  mutate(pred_dist = fitted(mod4)) %>% 
  mutate(L1_resid = residuals(mod4, type = "response")) %>% 
  # graph
  ggplot(aes(x=pred_dist, y=L1_resid)) + theme_classic() +
  geom_point(size=3, shape=1) + 
  geom_hline(yintercept=0)
隨機截距  每個個體的差異還挺大的
dotplot(ranef(mod4, condVar = T), strip = T, scales=list(relation='free'))$Subject
模型比較
選擇mod2
anova(mod1,mod2,mod3,mod4)
Data: data1
Models:
mod1: distance ~ (1 | Subject)
mod2: distance ~ age + (1 | Subject)
mod3: distance ~ age + (1 | Subject) + (0 + age | Subject)
mod4: distance ~ age + (1 + age | Subject)
     npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
mod1    3 521.49 529.54 -257.75   515.49                          
mod2    4 451.39 462.12 -221.69   443.39 72.1016  1    < 2e-16 ***
mod3    5 449.74 463.15 -219.87   439.74  3.6513  1    0.05603 .  
mod4    6 451.21 467.30 -219.61   439.21  0.5267  1    0.46801    
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1
- END -

相關焦點

  • 使用R語言進行混合線性模型(mixed linear model) 分析代碼及詳解
    1.混合線性模型簡介混合線性模型，又名多層線性模型（Hierarchical linear model）。
  • R語言|線性混合模型到底是什麼?該如何在自己的研究中進行模型搭建?
    想寫這個主題已經很長時間了，但是一直沒法理解線性混合模型，特別是它適用的場景，所以就一直在等找到合適的資料，然後理解清楚了再寫。
  • 線性混合模型系列二:模型假定
    這裡：觀測值：體重weight固定因子：Chang隨機因子：Sire殘差：e如果公牛之間沒有親緣關係，那麼他們之間的親緣關係矩陣是單位矩陣，它的分布為這裡b為固定因子的效應值，加入固定因子有多個，場，年，季，性別等等，那麼b 可以分解為：[b1, b2, b3,...]
  • 生態學模擬對廣義線性混合模型GLMM進行功率(功效、效能、效力)分析power analysis環境監測數據
    它包括用於 (i) 對給定模型和設計進行功效分析的工具；(ii) 計算功效曲線以評估功效和樣本量之間的權衡。本文提供了一個教程，使用具有混合效果的計數數據的簡單示例（具有代表環境監測數據的結構）。介紹假設檢驗的功效定義為假設原假設為假，檢驗拒絕原假設的概率。換句話說，如果一個效應是真實的，那麼分析判斷該效應具有統計顯著性的概率是多少？
  • 常用混合線性模型演示:lme4和asreml
    這裡使用sleepstudy數據集，看一下免費的R包lme4和付費包asreml如何處理不同的混合線性模型，以加深對混合線性模型的理解。1. 隨機斜率，相同截距通用的混線性模型，包括固定因子和隨機因子。育種中常用的混線性模型，一般是針對於隨機因子關係矩陣，而其它領域的一般是針對於不同截距的定義。
  • 直播丨基於R語言的結構方程模型分析及應用
    、結構方程模型R實現的流程、潛變量和組成變量分析、混合效應模型及貝葉斯方法在結構方程模型中的應用等主要環節。課程目標（1）了解結構方程模型的基本原理（2）掌握利用R實現結構方程模型的步驟和流程； （3）掌握利用R實現潛變量和組成變量分析的原理和流程；（4）掌握混合效應模型和貝葉斯方法結構方程分析中的應用；（5）通過理論知識學習與上機實踐操作，讓學員具備根據自己的研究問題構建結構方程模型的能力
  • R語言 | 一個多元線性回歸在R中的實現示例
    從技術上來說，簡單多項式回歸也可以視為多元線性回歸的特例，例如二次回歸有兩個預測變量（X1=X和X2=X2），三次回歸有三個預測變量（X1=X、X2=X2和X3=X3）。 本篇繼續通過一個示例，展示在R語言中擬合多元線性回歸的方法過程。
  • R語言LME4混合效應模型研究教師的受歡迎程度
    我們第一個模型是截距模型。## 通過REML進行線性混合模型擬合。>如果查看匯總輸出，則在「隨機效應」下可以看到，類別層0.7021上的殘差和第一層（學生層）上的殘差為1.2218。它檢查如果刪除了某種隨機效應（稱為似然比檢驗），則模型是否變得明顯更差，如果不是這種情況，則隨機效應不顯著。
  • 一元(多元)線性回歸分析之R語言實現
    上篇介紹了《一元（多元）線性回歸分析之Excel實現》，本篇來探討一下回歸分析在R語言中的實現，我們將從更專業的角度對模型進行一些解讀。
  • r語言 檢驗p值 - CSDN
    輸入1： rdata = matrix(rnorm(1000* 6, 0, 3), 6) rvar = apply(rdata, 2, var) mean(rvar)結果1： [1] 8輸入2： var(rvar)結果2： [1] 32=2*81/5輸入3：
  • R語言ggplot2統計作圖
    收錄於話題 #r語言1個
  • 如何用潛類別混合效應模型(Latent Class Mixed Model ,LCMM)分析老年痴呆年齡數據
    高斯數據示例在此示例中，我們研究了認知標記的二次軌跡，即在老年人樣本（納入時年齡 65 歲及以上）中進行預先標準化（具有高斯分布）並對簡易智能量表評分 ( MMSE )進行了長達 15 年的跟蹤研究，可根據教育水平進行調整。儘管可以考慮任何回歸，但模型在此處不考慮交互作用。數據集子樣本這是來自原始前瞻性研究 的 500 名受試者的子樣本。
  • R語言——交叉驗證法計算線性回歸模型擬合優度的第三種R方
    本來打算這周繼續更新一篇LaTex的小短文，但是貌似我已經很久沒有寫R語言相關的東西了。想來想去，今天就寫一篇和R語言有關的，畢竟不能忘記初心呀！凡是學過計量的同學，哪怕只記得一點點皮毛，對於R方和調整R方也應該是再熟悉不過了。R方和調整R方是判斷回歸模型擬合度的最為方便簡單的指標，一般來說，數值越大代表模型的擬合度越好。
  • 線性概率模型、probit模型和logit模型的聯繫與區別
    如果被解釋變量是離散的，而非連續值，則需要運用"離散選擇模型"（discrete choice
  • R語言 | 回歸分析(四)
    為了將非線性的結果轉變為線性內容，我們可以採用優勢比進行轉化。R: The R Project for Statistical Computinghttps://www.r-project.org/RStudio:https://rstudio.com/今天是我們最後一期R語言入門系列。
  • 機器學習 · R語言包大全(共99個包)
    有很多R語言包都可以實現機器學習相關的思想和方法。()函數可以進行參數分類模型；擁有強大的用於可視化二叉樹的功能 vcrpart：樹狀變係數模型LogicReg：Logic回歸樹maptree：樹狀圖的可視化工具REEMtree：縱向數據隨機效應樹模型RPMM：混合效應模型的分類算法partykit：用於樹模型的預測和可視化evtree：結合partykit
  • 線性模型(一)普通線性回歸到廣義線性模型
    同時提醒讀者避免只從字面理解「線性」帶來誤會，即線性模型只能解決線性問題。本章將線性模型定位和表述為在數學表達式上具有線性的表示方式的數學模型，包含普通線性回歸模型和廣義線性模型（線性支持向量機本章不進行講述）。
  • 一般線性模型與廣義線性模型
    你所接觸到的統計分析方法都是基於已有或者正在探索中的統計模型；醫學統計學，大致應用在醫學科研、藥物臨床試驗兩個領域，兩者有交叉，而後者更加注重實用性；臨床試驗中的統計分析方法基本上都來自於經過實證的數據分析模型。
  • 重複測量數據分析系列:廣義線性混合模型(GLMM)
    有多個統計模型可以實現重複測量數據的分析：【1】一般線性模型中的重複測量方差分析，可以採用一元方差分析和多元方差分析。重複測量方差分析要求還是比較苛刻的，要求多元正態性、組間方差-協方差矩陣相等（Box』M檢驗），數據上也不能有缺失值。一元差分析雖然考慮了個體隨機效應，但要求方差協方差結構滿足複合對稱或者球性假設。
  • R語言中COX模型構建
    該模型以生存結局和生存時間為應變量，可同時分析眾多因素對生存期的影響，能分析帶有截尾生存時間的資料，且不要求估計資料的生存分布類型。由於上述優良性質，該模型自問世以來，在醫學隨訪研究中得到廣泛的應用，是迄今生存分析中應用最多的多因素分析方法 [引自百度百科]。今天我們介紹下在R語言中COX模型是如何實現又是如何來評價準確性的。