Trong bài hôm nay, mình sẽ trình bày cách sử dụng vòng lặp (For loop) trong R để thi hành những quy trình mang tính tổ hợp và lặp lại.

Nhằm minh họa, chúng ta sẽ sử dụng dataset về kết quả kì thi PISA năm 2013 của Việt Nam, bộ số liệu này từng được giới thiệu trong Workshop của thầy Nguyễn Văn Tuấn, cũng như các bài thực hành về Data wrangling của bạn Nguyễn Chí Dũng trước đây.

df=read.csv("PISA DATA (INTERNATIONAL).CSV")%>%as_tibble()%>%subset(.,CNT == "Viet Nam")

df=df[,c(2:5,26:36)]%>%separate(STRATUM,into = c("Region", "School","Zone"), sep = "/")%>%separate(Region,into=c("Code","Region"),sep=":")%>%.[,-c(1,6)]%>%dplyr::rename(.,Gender=ST04Q01)%>%.[,c(4,1:3,5:16)]
df$PARED=as.numeric(df$PARED)

Mục tiêu nghiên cứu giả định của chúng ta là : Đánh giá hiệu ứng của mỗi yếu tố trong số 5 yếu tố là Vùng miền (cụ thể là miền Bắc), Trường công lập, Vùng nông thôn (so với thành thị), Giới tính Nam (so với Nữ) và trình độ học vấn của cha mẹ (PARED) đối với điểm của từng môn trong số 10 môn thi (MATH, MACC, MACQ, MACS, MACU, MAPE, MAPF, MAPI, READ, SCIE) của học sinh Việt Nam, trong đó có xét đến hiệu ứng ngẫu nhiên của Trường học (School ID)

Bài toán này yêu cầu chúng ta phải dựng 50 mô hình tương ứng với 50 cặp tổ hợp giữa Outcome và Predictor, nội dung mô hình có dạng :

Outcome ~ Predictor + (1|SCHOOLID)

Đây là 1 mixed linear model và ta có thể làm việc này dễ dàng trong package lme4, tuy nhiên không ai đủ kiên nhẫn ngồi viết code cho 50 mô hình và xử lý kết quả với 50 output objects khác nhau. Chính lúc này ta cần sử dụng vòng lặp (For loop).

Nguyên tắc của vòng lặp FOR trong R rất đơn giản, nó gồm 3 yếu tố

  1. Một hay nhiều object nhận kết quả, có thể là vector hay dataframe

  2. Điểm bắt đầu và kết thúc

  3. Nhiệm vụ cho mỗi lượt thi hành

Trong nghiên cứu này, chúng ta có tới 2 nhóm cần tổ hợp : Outcome (n=10) và Predictors (n=5), ta sẽ tạo ra danh sách, điểm bắt đầu và kết thúc của chúng:

outcome_start=7
outcome_end= 16
n_outcome=outcome_end-outcome_start+1

outcomelist=rep(NA, n_outcome)

pred_start=2
pred_end=6
n_pred=pred_end-pred_start+1

predictorlist=rep(NA,n_pred)

Sau đó, chúng ta sẽ nghĩ đến chuyện cần có bao nhiêu kết quả sẽ được khai thác trong vòng lặp ? Để minh họa, mình sẽ trích xuất tất cả những gì có thể từ 1 mô hình Mixed model, bao gồm cả những chỉ số đánh giá độ chính xác dự báo như rmse, mse, R2, BIC, cũng như khoảng tin cậy 97.5% của hệ số hồi quy và P value cho factor cần khảo sát.

Mỗi chỉ số sẽ được chuẩn bị sẵn 1 object để nhận kết quả từ vòng lặp. Như vậy ta có 12 objects cần chuẩn bị sẵn (10 kết quả + 1 danh sách Outcome và 1 danh sách predictor).

predbeta=rep(NA, n_pred)
predse = rep(NA, n_pred)
predpvalue=rep(NA, n_pred)
modrmse=rep(NA, n_pred)
modmse=rep(NA, n_pred)
modR2=rep(NA, n_pred)
modmae=rep(NA, n_pred)
LIC=rep(NA, n_pred)
UIC=rep(NA, n_pred)
modAIC=rep(NA, n_pred)
modBIC=rep(NA, n_pred)

Vòng lặp của chúng ta là sự tổ hợp của 1 vòng lặp lớn (10 lần cho Outcome) và 1 vòng lặp thứ cấp (cho 5 predictors).

number=1 #Điểm khởi đầu cho vị trí kết quả

library(lme4)

for (i in outcome_start:outcome_end) {
  outcome = colnames(df)[i]
  for (j in pred_start:pred_end){
    predictor = colnames(df)[j]
    
    # Mô hình Mixed linear
    
    model <- lmer(get(outcome) ~ get(predictor) + (1|SCHOOLID),
                  na.action = na.exclude,
                  data=df)
    
    # Bóc tách kết quả
    
    Vcov <- vcov(model, useScale = FALSE)
    beta <- fixef(model)
    se <- sqrt(diag(Vcov))
    zval <- beta / se
    pval <- 2 * pnorm(abs(zval), lower.tail = FALSE)
    
    l=confint(model)%>%.[4]
    u=confint(model)%>%.[8]
    AIC=AIC(model)
    BIC=BIC(model)
    
    newpred=predict(model,newdata=df)
    true=df[i]
    error=(true-newpred)
    rmse=sqrt(mean(error^2))
    mse=mean((newpred-true)^2)
    R2=1-(sum((true-newpred)^2)/sum((true-mean(true[[1]]))^2))
    
    # Gán kết quả
    predbeta[number] = as.numeric(beta[2])
    predse[number] = as.numeric(se[2])
    predpvalue[number] = as.numeric(pval[2])
    predictorlist[number] = predictor
    modrmse[number]= rmse
    modmse[number]=mse
    modR2[number]=R2
    LIC[number]=l
    UIC[number]=u
    modAIC[number]=AIC
    modBIC[number]=BIC
    outcomelist[number]=outcome
    
    number = number + 1
  }
}

Quy trình bên trong vòng lặp mới nhìn có vẻ phức tạp, nhưng có thể diễn giải theo trình tự như sau:

giá trị [i] chỉ thứ tự Outcome từ vị trí cột 7 đến 16 trong dataframe, [j] chỉ thứ tự Predictor từ vị trí ,cột 2 đến 6 trong dataframe và number chỉ vị trí trong vector object kết quả. Khi vòng lặp chạy, number tăng dần từ 1 đến 50, i tăng dần từ 7 – 16 và j tăng dần từ 2-6

Trước tiên, dựng một mô hình mixed linear có nội dung Outcome ~ Predictor + (1|SCHOOLID)

Sau đó bắt đầu khai thác kết quả mô hình, bao gồm :

Cuối cùng, gán các giá trị này cho từng object tương ứng đã chuẩn bị sẵn từ trước.

Sau khi vòng lặp chạy đủ 50 bước, mỗi object đã nhận đủ 50 kết quả, ta chỉ còn việc ghép nối chúng lại với nhau để có 1 dataframe kết quả.

CompiledDF= data.frame(outcomelist,predictorlist, predbeta,predse,LIC,UIC,predpvalue,modR2,modrmse,modmse,modAIC,modBIC)

CompiledDF=CompiledDF%>%rename(.,Exam=outcomelist,Factor=predictorlist,Coef=predbeta,SE=predse,LL=LIC,UL=UIC,Pvalue=predpvalue,R2=modR2,RMSE=modrmse,MSE=modmse,AIC=modAIC,BIC=modBIC)

Từ dataframe này ta có thể bắt đầu trình diễn kết quả bằng đồ thị như sau:

Hình 1 Với các bạn ưa chuộng giá trị p, đây là đồ thị so sánh giá trị p của 5 factor cho 10 outcome:

my_theme <- function(base_size = 10, base_family = "sans"){
  theme_minimal(base_size = base_size, base_family = base_family) +
    theme(
      axis.text = element_text(size = 10),
      axis.text.x = element_text(angle = 0, vjust = 0.5, hjust = 0.5),
      axis.title = element_text(size = 12),
      panel.grid.major = element_line(color = "grey"),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = "#faefff"),
      strip.background = element_rect(fill = "#400156", color = "#400156", size =0.5),
      strip.text = element_text(face = "bold", size = 10, color = "white"),
      legend.position = "bottom",
      legend.justification = "center",
      legend.background = element_blank(),
      panel.border = element_rect(color = "grey30", fill = NA, size = 0.5)
    )
}
theme_set(my_theme())

mycolors=c("#ff0180","#9809bc", "#0bc1c4" , "#088969", "#a2c40b")

CompiledDF%>%ggplot(aes(x=Factor,y=Pvalue,fill=Factor))+geom_point(size=3,shape=21,color="black",show.legend = F)+facet_wrap(~Exam,ncol=3,scales="free")+coord_flip()+scale_y_continuous(limits=c(0,0.07))+geom_hline(yintercept=0.05)+scale_fill_manual(values=mycolors)+scale_color_manual(values=mycolors)

Thí dụ với môn Khoa học (Scie), ta có thể diễn giải: hiệu ứng của Miền Bắc không có ý nghĩa (p> 0.05), Trường công lập, Giới tính Nam, Khu vực thành thị vs nông thôn, và trình độ học vấn cha mẹ có hiệu ứng ý nghĩa với điểm môn Khoa học (p<0.05).

Hình 2

Biểu đồ tiếp theo so sánh giá trị hệ số hồi quy giữa 5 yếu tố:

Thí dụ cho môn Đọc (Read), có thể diễn giải như sau: Trình độ học vấn của cha mẹ không gây hiệu ứng đáng kể cho điểm Đọc, Nam sinh có điểm đọc thấp hơn Nữ sinh, Học sinh miền Bắc không khác biệt đáng kể so với điểm đọc trung bình.

CompiledDF%>%ggplot(aes(x=Factor,y=Coef,fill=Factor))+geom_bar(stat="identity",show.legend = F,width=0.5)+facet_wrap(~Exam,ncol=3,scales="free")+coord_flip()+geom_hline(yintercept=0,color="blue4",linetype = "dashed")+scale_y_continuous(limits=c(-30,130))+scale_fill_manual(values=mycolors)+scale_color_manual(values=mycolors)

Hình 3

Tương tự, đây là biểu đồ trình bày khoảng tin cậy của hệ số hồi quy cho 5 yếu tố và trong 10 môn thi.

Nếu CI chứa giá trị 0, ta có thể nói là hiệu ứng của yếu tố đó không có ý nghĩa thống kê.

CompiledDF%>%ggplot()+geom_errorbar(stat="identity",position="identity",aes(x=reorder(Factor,Coef),ymin=LL,ymax=UL,y=Coef,color=Factor),size=1,width=.5,show.legend = F)+geom_hline(yintercept = 0,color="red",linetype="dashed")+geom_point(stat="identity",position="identity",shape=21,size=2,aes(x=reorder(Factor,Coef),y=Coef,fill=Factor),color="black",show.legend = F)+coord_flip()+scale_y_continuous(limits=c(-100,550))+facet_wrap(~Exam,ncol=3)+scale_fill_manual(values=mycolors)+scale_color_manual(values=mycolors)

Hình 4

Tiếp theo, ta so sánh giá trị hệ số R2 của 50 mô hình:

Kết quả này cho thấy 1 điều, đó là toàn bộ các mô hình đơn biến đều có độ chính xác rất thấp (R2 dao động trong khoảng 0.4 -0.5, điều này có nghĩa là không có yếu tố độc lập nào cho phép giải thích được hoàn toàn sự biến thiên của điểm số của các môn thi.

CompiledDF%>%ggplot(aes(x=Factor,y=R2,fill=Factor))+geom_point(size=2,shape=21,color="black",show.legend = F)+facet_wrap(~Exam,ncol=3,scales="free")+coord_flip()+geom_hline(yintercept=0.5,color="blue",linetype="dashed")+scale_fill_manual(values=mycolors)+scale_color_manual(values=mycolors)

Hình 5

Đây là 1 heatmap, mô tả cùng 1 thông tin về giá trị R2 của 50 mô hình: Màu sắc cho ta lợi thế cảm nhận trực quan, theo đó ta có thể diễn giải nhanh hơn, thí dụ:

3 môn MACS, SCIE, MAPF có độ chính xác dự báo Thấp nhất trong số 10 môn thi.

Giới tính là yếu tố có độ chính xác dự báo tương đối cao so với 4 yếu tố còn lại.

library(viridis)

CompiledDF%>%ggplot(aes(x=reorder(Exam,R2),y=Factor,fill=R2,label=round(R2,2)))+geom_tile(color="black")+geom_text(color = "white")+scale_fill_viridis(option="A",begin=0.2,end=0.8)+scale_x_discrete("Exams")

Hình 6

Tương tự, 1 heatmap khác trình bày giá trị p.

Đa số các yếu tố có ý nghĩa thống kê (p<0.05), ngoại trừ Miền Bắc (level 1 của Region)

CompiledDF%>%ggplot(aes(x=reorder(Exam,Pvalue),y=Factor,fill=Pvalue,label=round(Pvalue,2)))+geom_tile(color="black")+geom_text(color = "white")+scale_fill_viridis(option="D",begin=0.7,end=0.22)+scale_x_discrete("Exams")

Hình 7

Độ chính xác của mô hình còn có thể được so sánh bằng RMSE. Sự phân phối của màu sắc cho thấy sai số dự báo có vẻ như không phụ thuộc vào predictor mà liên quan đến Outcome (hay bản chất của môn thi). Theo heat map này, môn Đọc (Read) có độ chính xác dự báo cao nhất.

CompiledDF%>%ggplot(aes(x=reorder(Exam,RMSE),y=Factor,fill=RMSE,label=round(RMSE,1)))+geom_tile(color="black")+geom_text(color = "white")+scale_fill_viridis(option="A",begin=0.8,end=0.2)+scale_x_discrete("Exams")

Hình 8

Một heatmap khác trình bày độ mạnh của hiệu ứng (hệ số hồi quy) cho 5 yếu tố ở 10 môn thi:

Trình độ học vấn của cha mẹ có vẻ là tương quan nghịch với điểm thi (tuy nhiên hiệu ứng này rất yếu và không có ý nghĩa).

Trường công lập có hiệu ứng mạnh, làm tăng điểm thi tất cả các môn (tuy nhiên p value khá thấp).

Nữ chỉ giỏi hơn Nam ở môn Read, nhưng Nam hơn Nữ ở tất cả các môn còn lại

CompiledDF%>%ggplot(aes(x=reorder(Exam,Coef),y=Factor,fill=scale(CompiledDF$Coef),label=round(Coef,2)))+geom_tile(color="black")+geom_text(color = "white")+scale_fill_viridis(option="A",begin=0.8,end=0,name="Fixed effect")+scale_x_discrete("Exams")

Tổng kết

Trong bài này chúng ta đã học được rất nhiều chiêu thức hữu dụng, bao gồm:

Mở rộng

Vòng lặp trong bài sử dụng 1 mô hình mixed linear, tuy nhiên bạn có thể dùng nó cho bất cứ loại phân tích hồi quy nào khác.

Trong bài, tác giả đưa ra 1 trường hợp tổng quát, tổ hợp nhiều outcome với nhiều predictors, tuy nhiên bạn có thể cải biên lại vòng lặp để áp dụng cho những trường hợp khác, thí dụ dựng 1 mô hình duy nhất cho nhiều outcome, hoặc ngược lại: dựng nhiều mô hình đơn biến cho cùng 1 outcome.

Cả 2 vấn đề này khá thường gặp trên thực tế. Một ngày nào đó bạn sẽ phải giải quyết dữ liệu chứa hang trăm, thậm chí hàng ngàn outcome hoặc predictor (nghiên cứu genomic hoặc biomarker) bạn sẽ thấy vòng lặp vô cùng quan trọng.

Tạm biệt các bạn