1 load packages

pacman::p_load(mlmRev, tidyverse, lme4)#pacman::p_load()同時取代install.packages()及library(),功能與p_load()一樣

1.1 data from mlmRev package

data(Hsb82) 

1.2 first 10 lines

head(Hsb82,10)
##    school minrty     sx    ses   mAch   meanses sector        cses
## 1    1224     No Female -1.528  5.876 -0.434383 Public -1.09361702
## 2    1224     No Female -0.588 19.708 -0.434383 Public -0.15361702
## 3    1224     No   Male -0.528 20.349 -0.434383 Public -0.09361702
## 4    1224     No   Male -0.668  8.781 -0.434383 Public -0.23361702
## 5    1224     No   Male -0.158 17.898 -0.434383 Public  0.27638298
## 6    1224     No   Male  0.022  4.583 -0.434383 Public  0.45638298
## 7    1224     No Female -0.618 -2.832 -0.434383 Public -0.18361702
## 8    1224     No   Male -0.998  0.523 -0.434383 Public -0.56361702
## 9    1224     No Female -0.888  1.527 -0.434383 Public -0.45361702
## 10   1224     No   Male -0.458 21.521 -0.434383 Public -0.02361702

##show structure of data

str(Hsb82)
## 'data.frame':    7185 obs. of  8 variables:
##  $ school : Ord.factor w/ 160 levels "8367"<"8854"<..: 59 59 59 59 59 59 59 59 59 59 ...
##  $ minrty : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sx     : Factor w/ 2 levels "Male","Female": 2 2 1 1 1 1 2 1 2 1 ...
##  $ ses    : num  -1.528 -0.588 -0.528 -0.668 -0.158 ...
##  $ mAch   : num  5.88 19.71 20.35 8.78 17.9 ...
##  $ meanses: num  -0.434 -0.434 -0.434 -0.434 -0.434 ...
##  $ sector : Factor w/ 2 levels "Public","Catholic": 1 1 1 1 1 1 1 1 1 1 ...
##  $ cses   : num  -1.0936 -0.1536 -0.0936 -0.2336 0.2764 ...

2 make a copy of data to the dta object

dta <- Hsb82 

3 variable names for human

names(dta) <- c("School", "Minority", "Gender", "SES", "Math",
                "Ave_SES", "Sector", "C_SES")

4 see the random number generator

set.seed(20160923) #隨機抽樣1/3的raw data,()內為起始值

5 take a sample of size 31 from total number of schools

n <- sample(length(levels(dta$School)), 31) 

sample(x, size)可進行隨機抽樣,從length(levels(dta$School))中取出放回,重複31次。

6

theme_set(theme_bw())

7

ggplot(data = filter(dta, School %in% levels(School)[n]), 
       aes(C_SES, Math, color = Sector)) + 
 geom_smooth(method = "lm", se = F, lwd = .6) +
 geom_point(size = rel(0.5)) +
 facet_wrap(~ School, nrow = 2) +
 labs(x = "Centered SES", y="Math achievement score")+
 theme(axis.text.x = element_text(size = 6))
## `geom_smooth()` using formula 'y ~ x'

filter(dta, School %in% levels(School)[n])

filter用來篩選要分析的觀察值(Row),觀察值包含(%in%)抽樣的31所學校的level

7.1

8 one regression line per school

summary(m1b <- lm(Math ~ School / C_SES - 1, data = dta)) 
## 
## Call:
## lm(formula = Math ~ School/C_SES - 1, data = dta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.2662  -4.3241   0.1164   4.4434  18.3617 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## School8367        4.55279    1.61953   2.811 0.004950 ** 
## School8854        4.23978    1.07122   3.958 7.64e-05 ***
## School4458        5.81140    0.87465   6.644 3.28e-11 ***
## School5762        4.32486    0.99621   4.341 1.44e-05 ***
## School6990        5.97679    0.83237   7.180 7.68e-13 ***
## School5815        7.27136    1.21194   6.000 2.08e-09 ***
## School7172        8.06682    0.91354   8.830  < 2e-16 ***
## School4868       12.31018    1.03923  11.845  < 2e-16 ***
## School7341        9.79418    0.84853  11.543  < 2e-16 ***
## School1358       11.20623    1.10635  10.129  < 2e-16 ***
## School4383       11.46568    1.21194   9.461  < 2e-16 ***
## School2305       11.13776    0.74031  15.045  < 2e-16 ***
## School8800        7.33594    1.07122   6.848 8.13e-12 ***
## School3088        9.14585    0.97033   9.425  < 2e-16 ***
## School8775        9.46700    0.87465  10.824  < 2e-16 ***
## School7890        8.34110    0.84853   9.830  < 2e-16 ***
## School6144        8.54507    0.92410   9.247  < 2e-16 ***
## School6443        9.47553    1.10635   8.565  < 2e-16 ***
## School5192       10.40950    1.14518   9.090  < 2e-16 ***
## School6808        9.28611    0.91354  10.165  < 2e-16 ***
## School2818       13.87283    0.93504  14.837  < 2e-16 ***
## School9340       11.17855    1.12526   9.934  < 2e-16 ***
## School4523        8.35174    0.88390   9.449  < 2e-16 ***
## School6816       14.53824    0.81709  17.793  < 2e-16 ***
## School2277        9.29761    0.77587  11.983  < 2e-16 ***
## School8009       14.08472    0.88390  15.935  < 2e-16 ***
## School5783       13.15034    1.12526  11.686  < 2e-16 ***
## School3013       12.61083    0.83237  15.151  < 2e-16 ***
## School7101       11.84993    1.14518  10.348  < 2e-16 ***
## School4530        9.05570    0.76345  11.861  < 2e-16 ***
## School9021       14.69666    0.80976  18.149  < 2e-16 ***
## School4511       13.40903    0.79568  16.852  < 2e-16 ***
## School2639        6.61548    0.93504   7.075 1.64e-12 ***
## School3377        9.18671    0.90333  10.170  < 2e-16 ***
## School6578       11.99400    0.80976  14.812  < 2e-16 ***
## School9347       13.53875    0.80263  16.868  < 2e-16 ***
## School3705       10.33169    0.90333  11.437  < 2e-16 ***
## School3533       10.40904    0.87465  11.901  < 2e-16 ***
## School1296        7.63596    0.87465   8.730  < 2e-16 ***
## School4350       11.85542    1.05486  11.239  < 2e-16 ***
## School9397       10.35547    0.88390  11.716  < 2e-16 ***
## School4253        9.41286    0.79568  11.830  < 2e-16 ***
## School2655       12.34517    0.84033  14.691  < 2e-16 ***
## School7342       11.16641    0.79568  14.034  < 2e-16 ***
## School9292       10.27963    1.39020   7.394 1.59e-13 ***
## School3499       13.27653    0.98302  13.506  < 2e-16 ***
## School7364       14.17214    0.91354  15.513  < 2e-16 ***
## School8983       10.99200    0.84853  12.954  < 2e-16 ***
## School5650       14.27353    0.90333  15.801  < 2e-16 ***
## School2658       13.39616    0.90333  14.830  < 2e-16 ***
## School8188       12.74097    1.10635  11.516  < 2e-16 ***
## School4410       13.47298    0.94637  14.236  < 2e-16 ***
## School9508       13.57466    1.02428  13.253  < 2e-16 ***
## School8707       12.88394    0.87465  14.730  < 2e-16 ***
## School1499        7.66036    0.83237   9.203  < 2e-16 ***
## School8477       12.52224    0.99621  12.570  < 2e-16 ***
## School1288       13.51080    1.21194  11.148  < 2e-16 ***
## School6291       10.10731    1.02428   9.868  < 2e-16 ***
## School1224        9.71545    0.88390  10.992  < 2e-16 ***
## School4292       12.86435    0.75162  17.116  < 2e-16 ***
## School8857       15.29694    0.75747  20.195  < 2e-16 ***
## School3967       12.03508    0.84033  14.322  < 2e-16 ***
## School6415       11.86020    0.82462  14.383  < 2e-16 ***
## School1317       13.17769    0.87465  15.066  < 2e-16 ***
## School2629       14.90777    0.80263  18.574  < 2e-16 ***
## School4223       14.62262    0.90333  16.187  < 2e-16 ***
## School1462       10.49556    0.80263  13.076  < 2e-16 ***
## School9550       11.08914    1.12526   9.855  < 2e-16 ***
## School6464        7.09162    1.12526   6.302 3.12e-10 ***
## School4931       13.79081    0.79568  17.332  < 2e-16 ***
## School5937       16.77597    1.12526  14.908  < 2e-16 ***
## School7919       14.84997    0.99621  14.906  < 2e-16 ***
## School3716       10.36866    0.94637  10.956  < 2e-16 ***
## School1909       14.42332    1.14518  12.595  < 2e-16 ***
## School2651       11.08432    0.98302  11.276  < 2e-16 ***
## School2467       10.14752    0.84033  12.076  < 2e-16 ***
## School1374        9.72846    1.14518   8.495  < 2e-16 ***
## School6600       11.70389    0.80976  14.453  < 2e-16 ***
## School5667       13.77823    0.77587  17.758  < 2e-16 ***
## School5720       14.28230    0.83237  17.159  < 2e-16 ***
## School3498       16.39045    0.83237  19.691  < 2e-16 ***
## School3881       11.94922    0.94637  12.626  < 2e-16 ***
## School2995        9.54611    0.89346  10.684  < 2e-16 ***
## School5838       13.68961    1.08836  12.578  < 2e-16 ***
## School3688       14.65626    0.92410  15.860  < 2e-16 ***
## School9158        8.54517    0.83237  10.266  < 2e-16 ***
## School8946       10.37509    0.79568  13.039  < 2e-16 ***
## School7232       12.54263    0.84033  14.926  < 2e-16 ***
## School2917        7.97895    0.92410   8.634  < 2e-16 ***
## School6170       14.18105    1.32234  10.724  < 2e-16 ***
## School8165       16.45122    0.86567  19.004  < 2e-16 ***
## School9104       16.83211    0.81709  20.600  < 2e-16 ***
## School2030       12.07819    0.88390  13.665  < 2e-16 ***
## School8150       14.85236    0.91354  16.258  < 2e-16 ***
## School4042       14.31542    0.75747  18.899  < 2e-16 ***
## School8357       14.38185    1.16619  12.332  < 2e-16 ***
## School8531       13.52868    0.94637  14.295  < 2e-16 ***
## School6074       13.77909    0.80976  17.016  < 2e-16 ***
## School4420       13.87416    1.07122  12.952  < 2e-16 ***
## School1906       15.98317    0.83237  19.202  < 2e-16 ***
## School3992       14.64521    0.83237  17.595  < 2e-16 ***
## School3999       10.94404    0.89346  12.249  < 2e-16 ***
## School4173       12.72466    0.91354  13.929  < 2e-16 ***
## School4325       13.24000    0.83237  15.906  < 2e-16 ***
## School5761       11.13806    0.84033  13.254  < 2e-16 ***
## School6484       12.91240    1.02428  12.606  < 2e-16 ***
## School6897       15.09763    0.86567  17.440  < 2e-16 ***
## School7635       15.06553    0.84853  17.755  < 2e-16 ***
## School7734       10.55964    1.29194   8.173 3.54e-16 ***
## School8175       11.69809    1.05486  11.090  < 2e-16 ***
## School8874       12.05503    1.00995  11.936  < 2e-16 ***
## School9225       14.66733    1.00995  14.523  < 2e-16 ***
## School2458       13.98568    0.80263  17.425  < 2e-16 ***
## School3610       15.35495    0.75747  20.271  < 2e-16 ***
## School5640       13.16011    0.80263  16.396  < 2e-16 ***
## School3838       16.06281    0.82462  19.479  < 2e-16 ***
## School9359       15.27062    0.83237  18.346  < 2e-16 ***
## School2208       15.40467    0.78231  19.691  < 2e-16 ***
## School6089       15.56958    1.05486  14.760  < 2e-16 ***
## School1477       14.22847    0.76959  18.488  < 2e-16 ***
## School2768       10.88692    1.21194   8.983  < 2e-16 ***
## School3039       16.96386    1.32234  12.829  < 2e-16 ***
## School5819       12.13890    0.85697  14.165  < 2e-16 ***
## School6397       12.79610    0.78231  16.357  < 2e-16 ***
## School1308       16.25550    1.35500  11.997  < 2e-16 ***
## School1433       19.71914    1.02428  19.252  < 2e-16 ***
## School1436       18.11161    0.91354  19.826  < 2e-16 ***
## School1461       16.84264    1.05486  15.967  < 2e-16 ***
## School1637        7.02411    1.16619   6.023 1.80e-09 ***
## School1942       18.11090    1.12526  16.095  < 2e-16 ***
## School1946       12.90844    0.97033  13.303  < 2e-16 ***
## School2336       16.51770    0.88390  18.687  < 2e-16 ***
## School2526       17.05300    0.80263  21.246  < 2e-16 ***
## School2626       13.39661    0.98302  13.628  < 2e-16 ***
## School2755       16.47651    0.88390  18.641  < 2e-16 ***
## School2771       11.84411    0.81709  14.495  < 2e-16 ***
## School2990       18.44792    0.87465  21.092  < 2e-16 ***
## School3020       14.39527    0.78891  18.247  < 2e-16 ***
## School3152       13.20904    0.84033  15.719  < 2e-16 ***
## School3332       14.27816    0.98302  14.525  < 2e-16 ***
## School3351       11.46518    0.97033  11.816  < 2e-16 ***
## School3427       19.71559    0.86567  22.775  < 2e-16 ***
## School3657        9.52118    0.84853  11.221  < 2e-16 ***
## School4642       14.59903    0.77587  18.816  < 2e-16 ***
## School5404       15.41498    0.80263  19.206  < 2e-16 ***
## School5619       15.41624    0.74590  20.668  < 2e-16 ***
## School6366       15.65640    0.79568  19.677  < 2e-16 ***
## School6469       18.45572    0.80263  22.994  < 2e-16 ***
## School7011       13.81358    1.05486  13.095  < 2e-16 ***
## School7276       12.67940    0.83237  15.233  < 2e-16 ***
## School7332       14.63610    0.87465  16.734  < 2e-16 ***
## School7345       11.33855    0.80976  14.002  < 2e-16 ***
## School7688       18.42231    0.82462  22.340  < 2e-16 ***
## School7697       15.72178    1.07122  14.677  < 2e-16 ***
## School8193       16.23226    0.92410  17.565  < 2e-16 ***
## School8202       11.71243    1.02428  11.435  < 2e-16 ***
## School8627       10.88372    0.83237  13.076  < 2e-16 ***
## School8628       16.52838    0.77587  21.303  < 2e-16 ***
## School9198       19.09229    1.08836  17.542  < 2e-16 ***
## School9586       14.86369    0.78891  18.841  < 2e-16 ***
## School8367:C_SES  0.25037    2.22488   0.113 0.910403    
## School8854:C_SES  1.93884    1.35428   1.432 0.152292    
## School4458:C_SES  1.13184    1.36793   0.827 0.408034    
## School5762:C_SES -1.01410    1.95950  -0.518 0.604803    
## School6990:C_SES  0.94769    0.98169   0.965 0.334393    
## School5815:C_SES  3.01800    2.16485   1.394 0.163335    
## School7172:C_SES  0.99448    1.36612   0.728 0.466662    
## School4868:C_SES  1.28647    1.48571   0.866 0.386576    
## School7341:C_SES  1.70370    1.00528   1.695 0.090169 .  
## School1358:C_SES  5.06801    1.71168   2.961 0.003078 ** 
## School4383:C_SES  6.18019    2.11601   2.921 0.003504 ** 
## School2305:C_SES -0.78211    1.13763  -0.687 0.491795    
## School8800:C_SES  2.56813    1.42668   1.800 0.071893 .  
## School3088:C_SES  1.79134    1.23144   1.455 0.145807    
## School8775:C_SES  1.00146    1.32013   0.759 0.448114    
## School7890:C_SES -0.65597    1.44460  -0.454 0.649784    
## School6144:C_SES  2.77027    1.30158   2.128 0.033340 *  
## School6443:C_SES -0.74335    1.46027  -0.509 0.610734    
## School5192:C_SES  1.60349    1.46102   1.098 0.272452    
## School6808:C_SES  2.27610    1.15285   1.974 0.048384 *  
## School2818:C_SES  2.80241    1.52348   1.839 0.065888 .  
## School9340:C_SES  3.30950    2.00186   1.653 0.098334 .  
## School4523:C_SES  2.38079    1.24545   1.912 0.055971 .  
## School6816:C_SES  1.35272    1.16781   1.158 0.246767    
## School2277:C_SES -2.01503    1.17745  -1.711 0.087062 .  
## School8009:C_SES  1.55687    1.58675   0.981 0.326542    
## School5783:C_SES  3.11406    1.49402   2.084 0.037165 *  
## School3013:C_SES  3.83980    1.75094   2.193 0.028341 *  
## School7101:C_SES  1.29545    1.42842   0.907 0.364484    
## School4530:C_SES  1.64743    1.23926   1.329 0.183771    
## School9021:C_SES  2.52416    1.39478   1.810 0.070384 .  
## School4511:C_SES  0.04251    1.38066   0.031 0.975438    
## School2639:C_SES -0.63010    1.52971  -0.412 0.680417    
## School3377:C_SES -0.74685    1.31955  -0.566 0.571419    
## School6578:C_SES  2.39054    1.30725   1.829 0.067492 .  
## School9347:C_SES  2.68599    1.18051   2.275 0.022920 *  
## School3705:C_SES  1.15850    1.33049   0.871 0.383933    
## School3533:C_SES -0.31177    1.64469  -0.190 0.849658    
## School1296:C_SES  1.07596    1.36610   0.788 0.430948    
## School4350:C_SES  4.37192    1.59458   2.742 0.006127 ** 
## School9397:C_SES  2.44644    1.47583   1.658 0.097430 .  
## School4253:C_SES -0.39954    1.15329  -0.346 0.729023    
## School2655:C_SES  5.22704    1.35364   3.861 0.000114 ***
## School7342:C_SES  1.01246    1.42097   0.713 0.476173    
## School9292:C_SES  0.75868    2.74811   0.276 0.782501    
## School3499:C_SES  0.99238    1.58184   0.627 0.530445    
## School7364:C_SES  0.25950    1.82833   0.142 0.887139    
## School8983:C_SES  1.38594    1.34299   1.032 0.302119    
## School5650:C_SES  0.68062    1.17460   0.579 0.562308    
## School2658:C_SES  2.62990    1.42677   1.843 0.065335 .  
## School8188:C_SES  4.39759    1.65556   2.656 0.007920 ** 
## School4410:C_SES  2.76020    1.55700   1.773 0.076312 .  
## School9508:C_SES  3.95379    1.82178   2.170 0.030019 *  
## School8707:C_SES  3.39153    1.09903   3.086 0.002037 ** 
## School1499:C_SES  3.63473    1.16428   3.122 0.001805 ** 
## School8477:C_SES  3.81216    1.25548   3.036 0.002403 ** 
## School1288:C_SES  3.25545    1.84816   1.761 0.078205 .  
## School6291:C_SES  3.98087    1.69353   2.351 0.018770 *  
## School1224:C_SES  2.50858    1.42433   1.761 0.078243 .  
## School4292:C_SES -0.16061    1.16329  -0.138 0.890195    
## School8857:C_SES  0.80222    1.43452   0.559 0.576025    
## School3967:C_SES  3.31107    1.35069   2.451 0.014255 *  
## School6415:C_SES  3.53008    1.21524   2.905 0.003686 ** 
## School1317:C_SES  1.27391    1.58930   0.802 0.422837    
## School2629:C_SES  0.22235    1.14645   0.194 0.846225    
## School4223:C_SES  2.48659    1.50770   1.649 0.099140 .  
## School1462:C_SES -0.82881    1.39839  -0.593 0.553411    
## School9550:C_SES  3.89194    1.45938   2.667 0.007675 ** 
## School6464:C_SES  1.00349    1.82459   0.550 0.582348    
## School4931:C_SES  0.91185    1.17597   0.775 0.438133    
## School5937:C_SES  1.03962    2.00531   0.518 0.604175    
## School7919:C_SES  3.98937    1.88178   2.120 0.034042 *  
## School3716:C_SES  5.86379    1.18088   4.966 7.01e-07 ***
## School1909:C_SES  2.85479    1.63722   1.744 0.081260 .  
## School2651:C_SES  4.89906    1.55631   3.148 0.001652 ** 
## School2467:C_SES  3.13713    1.65830   1.892 0.058563 .  
## School1374:C_SES  3.85432    1.66836   2.310 0.020904 *  
## School6600:C_SES  4.70429    1.16489   4.038 5.44e-05 ***
## School5667:C_SES  3.52297    1.24071   2.839 0.004532 ** 
## School5720:C_SES  2.46631    1.26524   1.949 0.051302 .  
## School3498:C_SES -0.13109    1.45712  -0.090 0.928320    
## School3881:C_SES  2.39071    1.79271   1.334 0.182388    
## School2995:C_SES  1.43231    1.01804   1.407 0.159492    
## School5838:C_SES  1.85305    1.77992   1.041 0.297872    
## School3688:C_SES  1.53672    1.72035   0.893 0.371749    
## School9158:C_SES  3.86121    1.12532   3.431 0.000604 ***
## School8946:C_SES  1.69048    1.06749   1.584 0.113331    
## School7232:C_SES  5.00160    1.47738   3.385 0.000715 ***
## School2917:C_SES  1.13585    1.02982   1.103 0.270083    
## School6170:C_SES  4.81178    1.64962   2.917 0.003547 ** 
## School8165:C_SES  1.80224    1.38480   1.301 0.193148    
## School9104:C_SES  1.49398    1.83440   0.814 0.415430    
## School2030:C_SES  1.41198    1.31420   1.074 0.282681    
## School8150:C_SES -0.18571    1.51064  -0.123 0.902162    
## School4042:C_SES  1.69362    1.20501   1.405 0.159922    
## School8357:C_SES  2.67578    1.15717   2.312 0.020788 *  
## School8531:C_SES  3.31823    1.40287   2.365 0.018043 *  
## School6074:C_SES  1.52909    1.30292   1.174 0.240603    
## School4420:C_SES  2.95866    1.58878   1.862 0.062614 .  
## School1906:C_SES  2.14551    1.36955   1.567 0.117259    
## School3992:C_SES  0.53788    1.38802   0.388 0.698389    
## School3999:C_SES  3.56698    0.98836   3.609 0.000310 ***
## School4173:C_SES  3.36567    1.45956   2.306 0.021143 *  
## School4325:C_SES  2.75605    1.03273   2.669 0.007632 ** 
## School5761:C_SES  3.10801    1.19136   2.609 0.009106 ** 
## School6484:C_SES  0.60568    1.49351   0.406 0.685093    
## School6897:C_SES  3.58049    1.17477   3.048 0.002314 ** 
## School7635:C_SES  2.44847    1.26927   1.929 0.053767 .  
## School7734:C_SES  6.03523    1.48469   4.065 4.86e-05 ***
## School8175:C_SES  1.61237    1.54334   1.045 0.296186    
## School8874:C_SES  4.09630    1.43512   2.854 0.004326 ** 
## School9225:C_SES  2.88589    1.39642   2.067 0.038806 *  
## School2458:C_SES  2.95669    1.22988   2.404 0.016241 *  
## School3610:C_SES  2.95585    1.20868   2.446 0.014489 *  
## School5640:C_SES  3.82774    1.38890   2.756 0.005868 ** 
## School3838:C_SES  0.59790    1.23391   0.485 0.628005    
## School9359:C_SES -0.83348    1.46679  -0.568 0.569895    
## School2208:C_SES  2.63664    1.31898   1.999 0.045648 *  
## School6089:C_SES  1.69245    1.51926   1.114 0.265318    
## School1477:C_SES  1.23061    1.14871   1.071 0.284075    
## School2768:C_SES  3.51228    1.37517   2.554 0.010668 *  
## School3039:C_SES  2.95567    1.92045   1.539 0.123839    
## School5819:C_SES  1.97252    1.43955   1.370 0.170658    
## School6397:C_SES  2.75901    1.03849   2.657 0.007908 ** 
## School1308:C_SES  0.12602    2.89741   0.043 0.965308    
## School1433:C_SES  1.85429    1.74798   1.061 0.288809    
## School1436:C_SES  1.60056    1.61653   0.990 0.322150    
## School1461:C_SES  6.26650    1.54684   4.051 5.15e-05 ***
## School1637:C_SES  3.11681    1.56230   1.995 0.046081 *  
## School1942:C_SES  0.08938    2.17762   0.041 0.967260    
## School1946:C_SES  3.58583    1.40776   2.547 0.010881 *  
## School2336:C_SES  1.90497    1.39610   1.364 0.172457    
## School2526:C_SES  0.15950    1.32609   0.120 0.904264    
## School2626:C_SES  4.09968    1.77861   2.305 0.021197 *  
## School2755:C_SES  0.56050    1.41945   0.395 0.692950    
## School2771:C_SES  4.26819    1.60528   2.659 0.007859 ** 
## School2990:C_SES  1.32454    1.32327   1.001 0.316883    
## School3020:C_SES  1.65368    1.21797   1.358 0.174595    
## School3152:C_SES  2.76825    1.00158   2.764 0.005727 ** 
## School3332:C_SES  2.03095    1.71772   1.182 0.237106    
## School3351:C_SES  2.45504    1.16971   2.099 0.035867 *  
## School3427:C_SES -0.48817    1.31793  -0.370 0.711094    
## School3657:C_SES  3.73591    1.10047   3.395 0.000691 ***
## School4642:C_SES  3.27239    1.17344   2.789 0.005306 ** 
## School5404:C_SES  1.21423    1.60712   0.756 0.449955    
## School5619:C_SES  5.25753    1.25841   4.178 2.98e-05 ***
## School6366:C_SES  1.51752    1.21527   1.249 0.211812    
## School6469:C_SES  1.75529    1.29272   1.358 0.174562    
## School7011:C_SES  5.07465    1.72585   2.940 0.003289 ** 
## School7276:C_SES  3.77336    1.06358   3.548 0.000391 ***
## School7332:C_SES  2.46320    1.32282   1.862 0.062634 .  
## School7345:C_SES  4.21192    0.98954   4.256 2.10e-05 ***
## School7688:C_SES  0.11634    1.47469   0.079 0.937119    
## School7697:C_SES  3.13622    1.77439   1.767 0.077190 .  
## School8193:C_SES  2.33521    1.62482   1.437 0.150704    
## School8202:C_SES  3.70590    1.47606   2.511 0.012073 *  
## School8627:C_SES  1.86956    1.18737   1.575 0.115408    
## School8628:C_SES  1.23139    1.34892   0.913 0.361339    
## School9198:C_SES  2.61055    1.69902   1.537 0.124462    
## School9586:C_SES  1.67208    1.33730   1.250 0.211217    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.06 on 6865 degrees of freedom
## Multiple R-squared:  0.8328, Adjusted R-squared:  0.825 
## F-statistic: 106.8 on 320 and 6865 DF,  p-value: < 2.2e-16

lm(Math ~ -1 + School, data = dta)可得到每一所學校與數學成績關係的迴歸線。

lm(Math ~ School / C_SES - 1, data = dta)則比上一條式子多數學和C_SES的交互作用項(School / C_SES),語法與lm(Math ~ -1 + School/C_SES, data = dta)相同。

好奇-1是無截距項嗎?試著拿掉-1,看會不會出現截距項,結果無法執行。

Error in h(simpleError(msg, call)) : error in evaluating the argument ‘object’ in selecting a method for function ‘summary’: orthogonal polynomials cannot be represented accurately enough for 159 degrees of freedom

9 collect regression estimates from individual regression by school

m1b_betas <- m1b %>% 
             coef %>% #取得模型的係數
             matrix(ncol = 2) %>% #將y #將資料放入矩陣中,ncol=2是「直行數=2」。
             as.data.frame %>% #as.data.frame將矩陣轉換為資料框架
             rename( b0 = V1, b1 = V2) #rename將舊檔案更換新檔名,b0重新命名為V1,b1重新命名為V2

Q:不懂為何matrix(ncol = 2)中不需要設定nrow=160

10

colMeans(m1b_betas) #同時計算column的平均
##        b0        b1 
## 12.620755  2.201641

Q:rename(“oldname.txt”, “newname.txt”),前面放舊檔,後面放新檔。然而rename( b0 = V1, b1 = V2)得到的結果和自己預期的剛好相反….好難理解>"<

11

apply(m1b_betas, 2, sd) #第二個引數:1為列,2為行,第三個引數sd為所要應用的函數
##       b0       b1 
## 3.117651 1.630998

12

cor(m1b_betas) #計算相關係數
##            b0         b1
## b0 1.00000000 0.01346686
## b1 0.01346686 1.00000000

13 scatter plot of regression coefficient estimates

ggplot(m1b_betas, aes(b0, b1)) + #將160條回歸線的截距估計值b0放x軸,斜率估計值b1放Y軸,畫出散佈圖。
 geom_point(pch = 20) +
 stat_smooth(method = "lm", se = F) +
 stat_ellipse() +
 annotate("text", x = 5, y = 5, label = "r = 0.013") + 
 labs(x = "Intercept estimates", y = "Slope estimates")
## `geom_smooth()` using formula 'y ~ x'

14 random intercept only

summary(m3a <- lmer(Math ~ C_SES + (1 | School), data = dta))#C_SES是固定效果,(1 | School)是隨機截距項的R語法,意思是假設每個學校的截距不同,其中"1"代表截距項。
## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ C_SES + (1 | School)
##    Data: dta
## 
## REML criterion at convergence: 46724
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0969 -0.7322  0.0194  0.7572  2.9147 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  School   (Intercept)  8.672   2.945   
##  Residual             37.010   6.084   
## Number of obs: 7185, groups:  School, 160
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  12.6361     0.2445   51.68
## C_SES         2.1912     0.1087   20.17
## 
## Correlation of Fixed Effects:
##       (Intr)
## C_SES 0.000

15

m3a <- lme4::lmer(Math ~ C_SES + (1 | School), data=dta)
sjPlot::tab_model(m3a, show.p=FALSE, show.r2=FALSE)
  Math
Predictors Estimates CI
(Intercept) 12.64 12.16 – 13.12
C_SES 2.19 1.98 – 2.40
Random Effects
σ2 37.01
τ00 School 8.67
ICC 0.19
N School 160
Observations 7185

16 random intercept and slope - correlated

summary(m3b <- lmer(Math ~ C_SES + (C_SES | School), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ C_SES + (C_SES | School)
##    Data: dta
## 
## REML criterion at convergence: 46714.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.09680 -0.73193  0.01855  0.75386  2.89924 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  School   (Intercept)  8.681   2.9464       
##           C_SES        0.694   0.8331   0.02
##  Residual             36.700   6.0581       
## Number of obs: 7185, groups:  School, 160
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  12.6362     0.2445   51.68
## C_SES         2.1932     0.1283   17.10
## 
## Correlation of Fixed Effects:
##       (Intr)
## C_SES 0.009

C_SES + (C_SES | School)的意義是「相關的隨機截距與斜率」

17 random intercept and slope - uncorrelated

summary(m3c <- lmer(Math ~ C_SES + (1 | School) + 
                           (C_SES - 1 | School), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Math ~ C_SES + (1 | School) + (C_SES - 1 | School)
##    Data: dta
## 
## REML criterion at convergence: 46714.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.09607 -0.73095  0.01694  0.75391  2.90167 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  School   (Intercept)  8.6807  2.9463  
##  School.1 C_SES        0.6941  0.8331  
##  Residual             36.7002  6.0581  
## Number of obs: 7185, groups:  School, 160
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  12.6360     0.2445   51.68
## C_SES         2.1926     0.1283   17.09
## 
## Correlation of Fixed Effects:
##       (Intr)
## C_SES 0.000

18

#summary(m3c <- lme(Math ~ C_SES, 
 #         random = list(School = pdDiag(~ C_SES)), data = dta))

19 default residual plot

plot(m3c, xlab = "Fitted values", ylab = "Pearson residuals", 
     pch = 20, cex = .5, type = c("p", "g"))

20 model comparisons

anova(m3a, m3b, m3c)
## refitting model(s) with ML (instead of REML)
## Data: dta
## Models:
## m3a: Math ~ C_SES + (1 | School)
## m3c: Math ~ C_SES + (1 | School) + (C_SES - 1 | School)
## m3b: Math ~ C_SES + (C_SES | School)
##     npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)   
## m3a    4 46728 46756 -23360    46720                        
## m3c    5 46721 46755 -23356    46711 9.4196  1   0.002147 **
## m3b    6 46723 46764 -23356    46711 0.0134  1   0.907954   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

21 CIs for model parameters

confint(m3c, method="boot")
## Computing bootstrap confidence intervals ...
## 
## 3 message(s): boundary (singular) fit: see ?isSingular
## 2 warning(s): Model failed to converge with max|grad| = 0.00255544 (tol = 0.002, component 1) (and others)
##                  2.5 %    97.5 %
## .sig01       2.6100481  3.269883
## .sig02       0.3848003  1.141918
## .sigma       5.9425454  6.150370
## (Intercept) 12.1304752 13.088948
## C_SES        1.9525498  2.439487

22 The end