Day 5. Mô hình tăng trưởng tuyến tính và mô hình đa tầng

Việc 1: Đánh giá tiến bộ của học sinh

wil = read.csv("C:\\Thach\\VN trips\\VN trip 5 (Apr 2023)\\Datasets\\willett.csv", header = T)

1.1 Vẽ biểu đồ bánh tằm đánh giá tiến bộ của học sinh

library(ggplot2)
ggplot(data = wil, aes(x = time, y = y, group = id, col = id)) + geom_line() + theme(legend.position = "none")

1.2 Thực hiện mô hình tuyến tính

m.12 = lm(y ~ time, data = wil)
summary(m.12)
## 
## Call:
## lm(formula = y ~ time, data = wil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -88.374 -25.584   1.186  28.926  64.746 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  164.374      5.035   32.65   <2e-16 ***
## time          26.960      2.691   10.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.6 on 138 degrees of freedom
## Multiple R-squared:  0.421,  Adjusted R-squared:  0.4168 
## F-statistic: 100.3 on 1 and 138 DF,  p-value: < 2.2e-16

1.3 Thưc hiện mô hình tuyến tính cho từng học sinh

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
wil %>% group_by(id) %>% mutate(beta = lm(y ~ time)$coefficients[2]) %>% mutate(alpha = lm(y ~ time)$coefficients[1]) %>% summarize(mbeta=mean(beta), malpha=mean(alpha))
## # A tibble: 35 × 3
##       id mbeta malpha
##    <int> <dbl>  <dbl>
##  1     1  34.2  197. 
##  2     2  28.5  218  
##  3     3  47.9  151. 
##  4     4  21.9  206. 
##  5     5   9.9  200. 
##  6     6  29.6  166. 
##  7     7  32.6  160. 
##  8     8  37.9   92.9
##  9     9  22.7  139. 
## 10    10  26.8  220. 
## # ℹ 25 more rows

1.4 Thực hiện mô hình ảnh hưởng hỗn hợp

library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
m.14 = lmer(y ~ time + (1 + time | id), data = wil)
summary(m.14)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ time + (1 + time | id)
##    Data: wil
## 
## REML criterion at convergence: 1266.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.27205 -0.60035  0.02632  0.60730  1.58439 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 1198.8   34.62         
##           time         132.4   11.51    -0.45
##  Residual              159.5   12.63         
## Number of obs: 140, groups:  id, 35
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  164.374      6.119   26.86
## time          26.960      2.167   12.44
## 
## Correlation of Fixed Effects:
##      (Intr)
## time -0.489

Việc 2: Điểm năm học trước có ảnh hưởng như thế nào đến tiến bộ của học sinh?

2.1 Thực hiện mô hình tuyến tính

wil$ccovar = wil$covar - mean(wil$covar)
m.21 = lm(y ~ time + ccovar, data = wil)
summary(m.21)
## 
## Call:
## lm(formula = y ~ time + ccovar, data = wil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.522 -24.591   3.165  28.033  70.018 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 164.3743     4.9646  33.109   <2e-16 ***
## time         26.9600     2.6537  10.159   <2e-16 ***
## ccovar        0.5357     0.2410   2.223   0.0278 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.11 on 137 degrees of freedom
## Multiple R-squared:  0.4412, Adjusted R-squared:  0.433 
## F-statistic: 54.08 on 2 and 137 DF,  p-value: < 2.2e-16

2.2 Thực hiện mô hình tuyến tính có yếu tố tương tác

m.22 = lm(y ~ time + ccovar + time:ccovar, data = wil)
summary(m.22)
## 
## Call:
## lm(formula = y ~ time + ccovar + time:ccovar, data = wil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -88.767 -26.040   3.676  26.224  73.433 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 164.3743     4.9089  33.485   <2e-16 ***
## time         26.9600     2.6239  10.275   <2e-16 ***
## ccovar       -0.1136     0.3987  -0.285   0.7762    
## time:ccovar   0.4329     0.2131   2.031   0.0442 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 34.71 on 136 degrees of freedom
## Multiple R-squared:  0.4576, Adjusted R-squared:  0.4457 
## F-statistic: 38.25 on 3 and 136 DF,  p-value: < 2.2e-16

2.3 Thực hiện mô hình ảnh hưởng hỗn hợp

m.23 = lmer(y ~ time + ccovar + time:ccovar + ( 1 + time | id), data = wil)
summary(m.23)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ time + ccovar + time:ccovar + (1 + time | id)
##    Data: wil
## 
## REML criterion at convergence: 1260.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.24817 -0.61873  0.00428  0.61472  1.55688 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 1236.4   35.16         
##           time         107.2   10.36    -0.49
##  Residual              159.5   12.63         
## Number of obs: 140, groups:  id, 35
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 164.3743     6.2061  26.486
## time         26.9600     1.9939  13.521
## ccovar       -0.1136     0.5040  -0.225
## time:ccovar   0.4329     0.1619   2.673
## 
## Correlation of Fixed Effects:
##             (Intr) time   ccovar
## time        -0.522              
## ccovar       0.000  0.000       
## time:ccovar  0.000  0.000 -0.522

Việc 3: Yếu tố ảnh hưởng đến điểm số học sinh

pop = read.csv("C:\\Thach\\VN trips\\VN trip 5 (Apr 2023)\\Datasets\\popular.csv")

3.1 Mô tả bằng biểu đồ

a. Phân bố điểm của học sinh

ggplot(data = pop, aes(popular)) + geom_histogram(fill = "blue", col = "white")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

b. Điểm học sinh nam có khác học sinh nữ không?

ggplot(data = pop, aes(x = sex, y = popular, col = sex)) + geom_boxplot() + geom_jitter(alpha = 0.1)

c.1 Ngoại cảnh ảnh hưởng thế nào đến điểm học sinh?

ggplot(data = pop, aes(x = factor(extrov), y = popular, col = factor(extrov))) + geom_boxplot() + geom_jitter(alpha = 0.1) 

c.2 Ảnh hưởng của ngoại cảnh có khác biệt giữa nam và nữ không?

ggplot(data = pop, aes(x=factor(extrov), y = popular, col = sex)) + geom_boxplot() + geom_jitter(alpha = 0.1)

3.2 Dữ liệu này có gì đặc biệt?

3.3 Điểm của học sinh có khác nhau giữa các lớp không?

m.33 = lmer(popular ~ 1 + (1 | class), data = pop)
summary(m.33)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + (1 | class)
##    Data: pop
## 
## REML criterion at convergence: 6330.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5655 -0.6975  0.0020  0.6758  3.3175 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  class    (Intercept) 0.7021   0.8379  
##  Residual             1.2218   1.1053  
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.07786    0.08739    58.1

3.4 Điểm học sinh nam có khác nữ không?

a. Không khác biệt giữa các lớp

m.34a = lmer(popular ~ 1 + Csex + (1 | class), data = pop)
summary(m.34a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + Csex + (1 | class)
##    Data: pop
## 
## REML criterion at convergence: 5564.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6890 -0.6415  0.0040  0.6788  3.1586 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  class    (Intercept) 0.4839   0.6957  
##  Residual             0.8308   0.9115  
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.06951    0.07252   69.91
## Csex         1.35010    0.04403   30.66
## 
## Correlation of Fixed Effects:
##      (Intr)
## Csex -0.004

b. Có khác biệt giữa các lớp

m.34b = lmer(popular ~ 1 + Csex + (1 + Csex | class), data = pop)
summary(m.34b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + Csex + (1 + Csex | class)
##    Data: pop
## 
## REML criterion at convergence: 5559.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5848 -0.6405  0.0062  0.6664  3.1420 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  class    (Intercept) 0.48299  0.6950        
##           Csex        0.05461  0.2337   -0.30
##  Residual             0.81896  0.9050        
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.07291    0.07254   69.93
## Csex         1.35218    0.05030   26.88
## 
## Correlation of Fixed Effects:
##      (Intr)
## Csex -0.144

3.5 Ngoại cảnh ảnh hưởng như thế nào đến điểm số?

a. Không có khác biệt giữa các lớp

m.35a = lmer(popular ~ 1 + Cextrav + (1 | class), data = pop)
summary(m.35a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + Cextrav + (1 | class)
##    Data: pop
## 
## REML criterion at convergence: 5832.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0644 -0.7267  0.0165  0.7088  3.3587 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  class    (Intercept) 0.8406   0.9168  
##  Residual             0.9304   0.9646  
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.07824    0.09421   53.90
## Cextrav      0.48631    0.02015   24.13
## 
## Correlation of Fixed Effects:
##         (Intr)
## Cextrav 0.000

b. Có khác biệt giữa các lớp

m.35b = lmer(popular ~ 1 + Cextrav + (1 + Cextrav | class), data = pop)
summary(m.35b)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + Cextrav + (1 + Cextrav | class)
##    Data: pop
## 
## REML criterion at convergence: 5779.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1961 -0.7291  0.0146  0.6816  3.2217 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  class    (Intercept) 0.89178  0.9443        
##           Cextrav     0.02599  0.1612   -0.88
##  Residual             0.89492  0.9460        
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.03127    0.09702   51.86
## Cextrav      0.49286    0.02546   19.36
## 
## Correlation of Fixed Effects:
##         (Intr)
## Cextrav -0.552

3.6 Ngoại cảnh và giới tính của học sinh có ảnh hưởng như thế nào đến điểm số

m.36 = lmer(popular ~ 1 + Cextrav + Csex + (1 + Cextrav + Csex | class), data = pop)
## boundary (singular) fit: see help('isSingular')
summary(m.36)
## Linear mixed model fit by REML ['lmerMod']
## Formula: popular ~ 1 + Cextrav + Csex + (1 + Cextrav + Csex | class)
##    Data: pop
## 
## REML criterion at convergence: 4870.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.01902 -0.64955 -0.01055  0.67099  3.11761 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  class    (Intercept) 0.673941 0.82094             
##           Cextrav     0.029847 0.17276  -0.73      
##           Csex        0.005365 0.07325  -0.65 -0.03
##  Residual             0.552883 0.74356             
## Number of obs: 2000, groups:  class, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.02051    0.08411   59.69
## Cextrav      0.44300    0.02343   18.91
## Csex         1.24483    0.03728   33.39
## 
## Correlation of Fixed Effects:
##         (Intr) Cextrv
## Cextrav -0.533       
## Csex    -0.127 -0.065
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Việc 4: Ghi lại các lệnh và chia sẻ trên tài khoản RPubs của bạn