{r setup, include=FALSE} rm(list = ls()) library(lme4) library(arm) knitr::opts_chunk$set(echo = TRUE)

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r cars} #Read data

lmm.data <- read.table(“http://bayes.acs.unt.edu:8083/BayesContent/class/Jon/R_SC/Module9/lmm.data.txt”, header = TRUE, sep = “,”, na.strings = “NA”, dec = “.”, strip.white = TRUE)

lmm.data %>% head(n=6)

#Nun-multilevel model

lm.ols <- lm(data = lmm.data, extro ~ open + agree + social)

lm.ols %>% summary()

lm.ols %>% display()

#According P-value and R-square results, this model is not well fit at all. #Try to see the other measures such as AIC

lm.glm <- glm(formula = extro ~ open + agree + social, data = lmm.data)

lm.glm %>% summary()

#AIC result is very poor!

#Run another glm

lm.glm2 <- glm(formula = extro ~ open + agree + social + school, data = lmm.data) lm.glm2 %>% summary()

AIC(lm.glm2)

anova(lm.glm, lm.glm2, test = “F”)

#Above are the fixed effect model, we can try to test the interaction between school and class variable as follows:

lm.glm.interac <- glm(formula = extro ~ open + agree + social + school:class, data = lmm.data)

lm.glm.interac %>% summary()

#While using those approaches above we can not define which is the best. It’s time to be mix-effected models being useful. Using syntax (1|school) tells function to fit a linear model with a varying-intercept group effect using the school variable

lm.mix <- lmer(extro ~ open + agree + social + (1|school), data = lmm.data)

lm.mix %>% summary()

#Or fit the model with multiple groups effect terms

#Crossed effect thể hiện hiệu ứng ngẫu nhiên của các yếu tố giao nhau, hay không có cơ sở để lập luận rằng chúng lồng ghép với nhau. Ví dụ để khảo sát mức độ tương hợp giữa biến school và biến class, crossed effect = (1|school) + (1|class) chứ không phải là nested effect = (1|school:class) vì mỗi level của class là như nhau cho cả 3 level

lm.mix1 <- lmer(extro ~ open + agree + social + (1|school) + (1|class), data = lmm.data)

lm.mix1 %>% summary()

#Random effects được chia làm 2 loại chính: Nested Effect và Crossed Effect #Nested Effect thể hiện một yếu tố (F) lồng trong yếu tố khác (G), hiệu ứng của G bị giới hạn bên trọng một cấp bậc nhất định của G, ví dụ G1 nhưng không hiện diện ở cấp bậc khác như G2, G3,… Ví dụ (1|school/class) sẽ đánh giá sự thay đổi của hệ số các trường học, trong đó các yếu tố lớp học được lồng ghép bên trong yếu tố trường học

lm.mix2 <- lmer(formula = extro ~ open + agree + social + (1 | school/class), data = lmm.data)

summary(lm.mix2)

#Multilevel model hay Hierachical model phổ biến hơn trong khoa học xã hội và kinh tế. Mô hình này có dạng cấu trúc đa cấp cũng tương tự như hiệu ứng lồng ghép nested như: #Điểm thi đại học ~ Giới tính + (1|Thành phố) + (1|Thành phố:Trường) + (1:Thành phố:Trường:lớp)

#Fit a varying slope with lmer##################### #What if we want to explore the effect of different student level indicators as they vary accross classrooms. Instead of fitting unique models by schools (or school/class) we can fit a varying slope model. Here we modify our random effect term to include variables before the grouping terms (1 + open|school/class) tells R to fit a varying slope and varying intercept

#lm.mix3 <- lmer(formula = extro ~ open + agree + social + (1 + open|school/class), data = lmm.data)

#lm.mix3 %>% summary()

```

Including Plots

You can also embed plots, for example:

{r pressure, echo=FALSE} plot(pressure)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.