讀進資料

並將 ID 視為類別變項(factor)

dta <- read.csv('Chemistry100.csv',header=T)
head(dta)
##   schID stuID gender sciscore
## 1     1     1      M    6.875
## 2     1     2      M    5.875
## 3     1     3      M    8.000
## 4     1     4      F    6.625
## 5     1     5      M    6.125
## 6     1     6      M    7.250
dta$schID <- as.factor(dta$schID)

載入 lme4 套件,用來分析多層次資料

library(pacman)
## Warning: 套件 'pacman' 是用 R 版本 4.2.3 來建造的
p_load(lme4,lmerTest,ggplot2)

分開做 100 次迴歸

rstsep <- lmList(sciscore ~ gender | schID, data=dta)
IntSlope <- as.data.frame(coef(rstsep))
names(IntSlope) <- c('Intercept','genderSlope')
head(IntSlope)
##   Intercept genderSlope
## 1  6.180778  0.43907937
## 2  6.656429 -0.63672269
## 3  5.747071  0.05865584
## 4  7.500000 -0.30084000
## 5  6.309667 -0.65900000
## 6  6.739806 -0.37417398

注意一下第一列

200 個參數的迴歸

rstfix <- lm(sciscore~schID+gender+schID:gender,data=dta)
summary(rstfix)
## 
## Call:
## lm(formula = sciscore ~ schID + gender + schID:gender, data = dta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.13630 -0.49724  0.03323  0.52237  2.09105 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.180778   0.256983  24.051  < 2e-16 ***
## schID2            0.475651   0.388521   1.224 0.220937    
## schID3           -0.433706   0.329385  -1.317 0.188020    
## schID4            1.319222   0.339956   3.881 0.000106 ***
## schID5            0.128889   0.406325   0.317 0.751106    
## schID6            0.559028   0.287315   1.946 0.051771 .  
## schID7            0.619008   0.329385   1.879 0.060289 .  
## schID8            0.415347   0.374613   1.109 0.267620    
## schID9           -0.087248   0.317809  -0.275 0.783693    
## schID10          -0.563340   0.321228  -1.754 0.079568 .  
## schID11          -0.673778   0.299691  -2.248 0.024623 *  
## schID12          -0.817965   0.321228  -2.546 0.010928 *  
## schID13          -0.364298   0.299691  -1.216 0.224227    
## schID14          -0.028000   0.363428  -0.077 0.938593    
## schID15           1.009222   0.354226   2.849 0.004410 ** 
## schID16          -0.025678   0.354226  -0.072 0.942216    
## schID17           0.770949   0.346515   2.225 0.026154 *  
## schID18          -0.185473   0.303120  -0.612 0.540656    
## schID19           0.058012   0.311965   0.186 0.852490    
## schID20           0.381815   0.296738   1.287 0.198281    
## schID21           0.430918   0.303120   1.422 0.155229    
## schID22           0.576710   0.283790   2.032 0.042211 *  
## schID23           0.050197   0.284428   0.176 0.859923    
## schID24          -0.025009   0.334305  -0.075 0.940372    
## schID25           0.066139   0.339956   0.195 0.845756    
## schID26           0.066739   0.269663   0.247 0.804543    
## schID27           0.242299   0.334305   0.725 0.468633    
## schID28          -0.119421   0.329385  -0.363 0.716958    
## schID29          -0.143013   0.317809  -0.450 0.652740    
## schID30           0.110131   0.346515   0.318 0.750636    
## schID31           0.205915   0.334305   0.616 0.537969    
## schID32           0.214589   0.293005   0.732 0.463990    
## schID33           0.728131   0.346515   2.101 0.035686 *  
## schID34           0.139431   0.301339   0.463 0.643606    
## schID35           0.348068   0.334305   1.041 0.297869    
## schID36           0.130125   0.283790   0.459 0.646604    
## schID37           0.452089   0.325060   1.391 0.164379    
## schID38           0.384847   0.374613   1.027 0.304341    
## schID39           0.282495   0.346515   0.815 0.414987    
## schID40          -0.088153   0.374613  -0.235 0.813977    
## schID41           0.007722   0.354226   0.022 0.982609    
## schID42           0.355131   0.346515   1.025 0.305498    
## schID43           0.025822   0.354226   0.073 0.941892    
## schID44          -0.013349   0.329385  -0.041 0.967675    
## schID45           0.242222   0.374613   0.647 0.517938    
## schID46           0.337556   0.325060   1.038 0.299137    
## schID47           0.859422   0.354226   2.426 0.015308 *  
## schID48           0.283722   0.374613   0.757 0.448877    
## schID49          -0.120278   0.339956  -0.354 0.723507    
## schID50           0.231275   0.311965   0.741 0.458532    
## schID51           0.388530   0.298161   1.303 0.192631    
## schID52           0.348340   0.317809   1.096 0.273125    
## schID53           0.320064   0.311965   1.026 0.304980    
## schID54           0.090972   0.339956   0.268 0.789023    
## schID55          -0.298325   0.277947  -1.073 0.283203    
## schID56          -0.456444   0.406325  -1.123 0.261367    
## schID57           0.979097   0.301339   3.249 0.001168 ** 
## schID58           0.752222   0.363428   2.070 0.038545 *  
## schID59           0.263508   0.329385   0.800 0.423765    
## schID60           0.059639   0.339956   0.175 0.860751    
## schID61           0.102937   0.388521   0.265 0.791068    
## schID62          -0.006178   0.354226  -0.017 0.986086    
## schID63           0.522556   0.339956   1.537 0.124353    
## schID64           0.510294   0.329385   1.549 0.121416    
## schID65           0.541972   0.374613   1.447 0.148056    
## schID66           0.247389   0.339956   0.728 0.466841    
## schID67           0.489677   0.346515   1.413 0.157702    
## schID68           0.852556   0.363428   2.346 0.019038 *  
## schID69           0.336163   0.317809   1.058 0.290241    
## schID70           0.444222   0.374613   1.186 0.235776    
## schID71           0.594222   0.430014   1.382 0.167101    
## schID72           0.045657   0.303120   0.151 0.880281    
## schID73           0.077079   0.307153   0.251 0.801869    
## schID74           0.652556   0.339956   1.920 0.054999 .  
## schID75           0.346915   0.334305   1.038 0.299473    
## schID76          -0.005778   0.354226  -0.016 0.986987    
## schID77           0.379607   0.334305   1.136 0.256240    
## schID78           0.273738   0.291913   0.938 0.348442    
## schID79           0.243060   0.286538   0.848 0.396348    
## schID80          -0.140159   0.307153  -0.456 0.648191    
## schID81           0.649667   0.363428   1.788 0.073926 .  
## schID82          -0.236194   0.339956  -0.695 0.487240    
## schID83          -0.346232   0.346515  -0.999 0.317775    
## schID84           0.129222   0.339956   0.380 0.703883    
## schID85          -0.266515   0.311965  -0.854 0.392992    
## schID86           0.510660   0.321228   1.590 0.111991    
## schID87          -0.067361   0.339956  -0.198 0.842942    
## schID88           0.043365   0.307153   0.141 0.887733    
## schID89           0.535722   0.321228   1.668 0.095459 .  
## schID90           0.289778   0.363428   0.797 0.425305    
## schID91           0.237530   0.334305   0.711 0.477430    
## schID92           0.706972   0.339956   2.080 0.037635 *  
## schID93           0.174741   0.296738   0.589 0.555985    
## schID94           0.737472   0.339956   2.169 0.030126 *  
## schID95          -0.028992   0.329385  -0.088 0.929867    
## schID96           0.460422   0.354226   1.300 0.193756    
## schID97           0.869056   0.339956   2.556 0.010619 *  
## schID98           0.783508   0.388521   2.017 0.043809 *  
## schID99           0.081583   0.287315   0.284 0.776465    
## schID100          0.510313   0.305052   1.673 0.094441 .  
## genderM           0.439079   0.329385   1.333 0.182609    
## schID2:genderM   -1.075802   0.477876  -2.251 0.024434 *  
## schID3:genderM   -0.380424   0.452749  -0.840 0.400824    
## schID4:genderM   -0.739919   0.426379  -1.735 0.082766 .  
## schID5:genderM   -1.098079   0.497170  -2.209 0.027263 *  
## schID6:genderM   -0.813253   0.395331  -2.057 0.039746 *  
## schID7:genderM   -0.526365   0.447748  -1.176 0.239842    
## schID8:genderM   -1.426051   0.478027  -2.983 0.002872 ** 
## schID9:genderM   -0.321819   0.418018  -0.770 0.441429    
## schID10:genderM  -0.728987   0.424976  -1.715 0.086367 .  
## schID11:genderM  -0.161437   0.417999  -0.386 0.699363    
## schID12:genderM  -0.869225   0.412804  -2.106 0.035305 *  
## schID13:genderM  -0.828559   0.399105  -2.076 0.037963 *  
## schID14:genderM  -1.168274   0.473354  -2.468 0.013632 *  
## schID15:genderM  -1.382979   0.444577  -3.111 0.001881 ** 
## schID16:genderM  -0.797329   0.444577  -1.793 0.072986 .  
## schID17:genderM  -0.584070   0.440238  -1.327 0.184690    
## schID18:genderM  -0.668191   0.396480  -1.685 0.092018 .  
## schID19:genderM  -0.902226   0.426885  -2.114 0.034628 *  
## schID20:genderM  -0.753731   0.377044  -1.999 0.045680 *  
## schID21:genderM  -0.906346   0.420465  -2.156 0.031184 *  
## schID22:genderM  -1.101276   0.384390  -2.865 0.004195 ** 
## schID23:genderM  -0.897721   0.369542  -2.429 0.015179 *  
## schID24:genderM  -1.025582   0.440271  -2.329 0.019893 *  
## schID25:genderM  -1.088150   0.451381  -2.411 0.015973 *  
## schID26:genderM  -0.629647   0.348104  -1.809 0.070569 .  
## schID27:genderM  -0.896862   0.434945  -2.062 0.039280 *  
## schID28:genderM  -0.451068   0.426885  -1.057 0.290744    
## schID29:genderM  -1.101844   0.444398  -2.479 0.013207 *  
## schID30:genderM  -1.067488   0.471130  -2.266 0.023524 *  
## schID31:genderM  -1.587688   0.451381  -3.517 0.000441 ***
## schID32:genderM  -1.037495   0.377893  -2.745 0.006073 ** 
## schID33:genderM  -0.868211   0.442207  -1.963 0.049684 *  
## schID34:genderM  -0.963636   0.398875  -2.416 0.015748 *  
## schID35:genderM  -0.894396   0.434945  -2.056 0.039823 *  
## schID36:genderM  -0.893416   0.366340  -2.439 0.014787 *  
## schID37:genderM  -0.705846   0.455581  -1.549 0.121393    
## schID38:genderM  -0.881093   0.464553  -1.897 0.057957 .  
## schID39:genderM  -0.835424   0.452749  -1.845 0.065089 .  
## schID40:genderM  -0.536441   0.462679  -1.159 0.246362    
## schID41:genderM  -1.032125   0.471130  -2.191 0.028536 *  
## schID42:genderM  -0.773930   0.444398  -1.742 0.081680 .  
## schID43:genderM  -0.398537   0.442982  -0.900 0.368359    
## schID44:genderM  -0.466294   0.439776  -1.060 0.289082    
## schID45:genderM  -0.671541   0.478027  -1.405 0.160164    
## schID46:genderM  -0.608713   0.421707  -1.443 0.148983    
## schID47:genderM  -1.304126   0.462223  -2.821 0.004808 ** 
## schID48:genderM  -0.468579   0.478027  -0.980 0.327037    
## schID49:genderM  -0.914496   0.455581  -2.007 0.044793 *  
## schID50:genderM  -0.571372   0.404414  -1.413 0.157792    
## schID51:genderM  -0.827923   0.390617  -2.120 0.034116 *  
## schID52:genderM  -1.346157   0.408939  -3.292 0.001005 ** 
## schID53:genderM  -1.050614   0.430693  -2.439 0.014763 *  
## schID54:genderM  -0.574529   0.433293  -1.326 0.184939    
## schID55:genderM  -0.695967   0.364183  -1.911 0.056082 .  
## schID56:genderM  -0.389886   0.488709  -0.798 0.425047    
## schID57:genderM  -0.811954   0.389658  -2.084 0.037254 *  
## schID58:genderM  -0.752913   0.473354  -1.591 0.111792    
## schID59:genderM  -0.607365   0.443473  -1.370 0.170911    
## schID60:genderM  -0.856223   0.460497  -1.859 0.063062 .  
## schID61:genderM  -0.868579   0.485651  -1.788 0.073785 .  
## schID62:genderM  -0.468367   0.452855  -1.034 0.301089    
## schID63:genderM  -0.824231   0.460497  -1.790 0.073561 .  
## schID64:genderM  -0.932951   0.425049  -2.195 0.028234 *  
## schID65:genderM  -0.879123   0.466639  -1.884 0.059655 .  
## schID66:genderM  -0.672913   0.473354  -1.422 0.155236    
## schID67:genderM  -1.373677   0.452749  -3.034 0.002430 ** 
## schID68:genderM  -0.859623   0.453670  -1.895 0.058199 .  
## schID69:genderM  -0.712754   0.427879  -1.666 0.095846 .  
## schID70:genderM  -1.410233   0.478027  -2.950 0.003198 ** 
## schID71:genderM  -0.485785   0.512180  -0.948 0.342958    
## schID72:genderM  -0.755969   0.384505  -1.966 0.049368 *  
## schID73:genderM  -1.506937   0.414439  -3.636 0.000281 ***
## schID74:genderM  -1.127975   0.441783  -2.553 0.010715 *  
## schID75:genderM  -0.800301   0.434945  -1.840 0.065852 .  
## schID76:genderM  -0.841352   0.471130  -1.786 0.074215 .  
## schID77:genderM  -0.918082   0.414362  -2.216 0.026779 *  
## schID78:genderM  -0.843295   0.373570  -2.257 0.024045 *  
## schID79:genderM  -0.770765   0.370782  -2.079 0.037712 *  
## schID80:genderM  -0.546763   0.394931  -1.384 0.166308    
## schID81:genderM  -1.532857   0.462773  -3.312 0.000935 ***
## schID82:genderM  -1.338218   0.473354  -2.827 0.004724 ** 
## schID83:genderM  -1.512425   0.471130  -3.210 0.001338 ** 
## schID84:genderM  -0.962767   0.441783  -2.179 0.029378 *  
## schID85:genderM  -0.712161   0.408403  -1.744 0.081287 .  
## schID86:genderM  -1.487384   0.430425  -3.456 0.000556 ***
## schID87:genderM  -0.960829   0.455581  -2.109 0.035014 *  
## schID88:genderM  -0.714572   0.408063  -1.751 0.080011 .  
## schID89:genderM  -0.871198   0.417067  -2.089 0.036792 *  
## schID90:genderM  -0.477635   0.462773  -1.032 0.302089    
## schID91:genderM  -1.546054   0.469313  -3.294 0.000997 ***
## schID92:genderM  -0.531488   0.415356  -1.280 0.200773    
## schID93:genderM  -0.902646   0.398510  -2.265 0.023571 *  
## schID94:genderM  -1.017329   0.466327  -2.182 0.029207 *  
## schID95:genderM  -0.546683   0.452749  -1.207 0.227330    
## schID96:genderM  -1.044221   0.450436  -2.318 0.020493 *  
## schID97:genderM  -0.549463   0.433293  -1.268 0.204843    
## schID98:genderM  -1.083626   0.468236  -2.314 0.020710 *  
## schID99:genderM  -0.355544   0.381444  -0.932 0.351349    
## schID100:genderM -1.550863   0.425712  -3.643 0.000273 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7709 on 3500 degrees of freedom
## Multiple R-squared:  0.2405, Adjusted R-squared:  0.1973 
## F-statistic:  5.57 on 199 and 3500 DF,  p-value: < 2.2e-16

注意一下截距與斜率

混合效果迴歸

rstmix <- lmer(sciscore~1+gender+(1+gender|schID),data=dta)
summary(rstmix)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sciscore ~ 1 + gender + (1 + gender | schID)
##    Data: dta
## 
## REML criterion at convergence: 8810
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9698 -0.6464  0.0519  0.6983  2.6584 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  schID    (Intercept) 0.10330  0.3214       
##           genderM     0.01994  0.1412   0.20
##  Residual             0.59449  0.7710       
## Number of obs: 3700, groups:  schID, 100
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  6.40435    0.03795 97.92954  168.76   <2e-16 ***
## genderM     -0.39534    0.03020 76.38063  -13.09   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## genderM -0.262

效果不等於100個迴歸的直接綜合

summary(IntSlope)
##    Intercept      genderSlope     
##  Min.   :5.363   Min.   :-1.1486  
##  1st Qu.:6.179   1st Qu.:-0.6001  
##  Median :6.423   Median :-0.4107  
##  Mean   :6.416   Mean   :-0.4211  
##  3rd Qu.:6.660   3rd Qu.:-0.2195  
##  Max.   :7.500   Max.   : 0.4391
var(IntSlope)
##               Intercept genderSlope
## Intercept    0.14012705 -0.02958089
## genderSlope -0.02958089  0.10575264

給第 50 個學校的男學生,看看預測結果

newdta1 <- data.frame(schID='50',gender='M') 

predict(rstfix,newdta1)
##       1 
## 6.27976
predict(rstmix,newdta1,re.form=~(1+gender|schID))
##        1 
## 6.191668

以第 50 個迴歸式預測

IntSlope[50,1] + IntSlope[50,2]*1
## [1] 6.27976

給第 101 個學校的男學生,看看預測結果

newdta2 <- data.frame(schID='101',gender='M') 

predict(rstfix,newdta2)
## Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels): 因子 schID 裡出現了新的層次 101
predict(rstmix,newdta2,allow.new.levels = TRUE,re.form=~(1+gender|schID))
##        1 
## 6.009013