Latent variable의 group간 비교: R 사용(updated)
문정훈 교수님께서 궁금하게 여기시는 latent variable의 그룹 간 평균값을 비교하는 SEM에 대해 소개합니다.
두 가지 접근법이 가능합니다. 보다 자세한 내용은 Qureshi & Compeau (2009) 논문을 참고하시면 됩니다. 두 방법은 Parametric Approach (PA)와 Mean and Covariance Structure (MACS, 혹은 Covariance SEM Approach: CSA)입니다. Qureshi & Compeau (2009)는 정규성이 명확하고, 샘플 사이즈가 작을 때는 PA가 MACS보다 더 좋지만, 여하간 두 결과 모두 타당한 결론에 이른다고 주장합니다. 다만, 정규성 가정이 흔들리면 둘 다 믿기 어렵다는 결론에 이릅니다.
우선, {plspm}으로 PLS 모델을 구축한 다음 formative type의 latent variable에서 PA로 그룹 간 평균을 비교해보고, {lavaan}으로 LISREL 방식의 공분산구조모형을 구축하여 MACS 방법을 사용하는 방법을 보여드립니다.
저는 R을 주로 사용합니다. 아래 예제에서 R>은 R 콘솔을 의미합니다.
데이터 준비
데이터는 CSV 파일로 있어야 하고, 제일 첫 머리는 variable의 이름을 나타냅니다. 우선 데이터가 다음과 같이 Excel로 준비되어 있다면,
파일 이름이 정확한지 확인하고(예제는 testdata.csv),
R> sample_data <- read.csv("testdata.csv")
데이터를 확인하려면
R> View(sample_data)
R로 데이터를 잘 불러왔습니다.
Parametric Approach
이 방법은 PLS를 이용합니다. 구체적으로 다음과 같이 group 간의 t-통계량을 구하는 것이 목표입니다.
보시다시피, 그룹 간 t test합니다. 표준편차 S_p는 각 그룹의 bootstraping standard error의 결합으로 정의합니다.
따라서 t분포는 자유도가 n_1 + n_2 - 2 가 됩니다.
예제는 패키지 plspm의 college를 활용하겠습니다. 사례는 Sanchez (2013)을 참고했습니다.
R> install.packages("plspm")
R> library(plspm)
우선 college가 csv 파일로 있다 가정하고 데이터를 가져오겠습니다.
R> college <- read.csv("college_data.csv",row.names=F)
path 구조가 다음과 같다고 합시다.
path 구조를 다음과 같이 작성합니다.
R> HighSchool = c(0, 0, 0, 0)
R> Intro = c(1, 0, 0, 0)
R> Medium = c(1, 1, 0, 0)
R> Graduation = c(1, 1, 1, 0)
R> gpa_path = rbind(HighSchool, Intro, Medium, Graduation)
관측변수의 위치를 HighSchool, Intro, Medium, Graduation 별로 지정합니다. 즉 college 데이터 프레임에서 1,2,3번은 HighSchool, 4,6,7번은 Intro, 8,9,10,11은 Medium, 마지막으로 12번 변수는 Graduation의 관측변수입니다.
R> gpa_blocks = list(1:3, 4:7, 8:11, 12)
관측변수와 잠재변수의 모드(mode)를 전부 formative type으로 지정합시다. Formative type은 A, reflective type은 B입니다.
R> gpa_modes = c("A","A","A","A")
이제 pls 모형을 fitting 시키겠습니다.
R> gpa_pls = plspm(college, gpa_path, gpa_blocks, modes = gpa_modes,
boot.val = TRUE)
R> summary(gpa_pls)
결과의 일부를 보면...
이제 PA 방법을 통해 그룹 간 t-value 차이를 계산해보겠습니다.
> gpa_boot = plspm.groups(gpa_pls, college$Gender, method = "bootstrap")
> gpa_boot
GROUP COMPARISON IN PLS-PM FOR PATH COEFFICIENTS
Scale of Data: TRUE
Weighting Scheme: centroid
Selected method: bootstrap
Num of replicates: 100
$test
global group.FEMALE group.MALE diff.abs t.stat
HighSchool->Intro 0.4531 0.4278 0.5000 0.0722 0.9125
HighSchool->Medium -0.0460 -0.0679 0.0107 0.0786 1.0220
HighSchool->Graduation -0.0488 -0.0432 -0.0628 0.0197 0.4246
Intro->Medium 0.7996 0.8140 0.7775 0.0365 0.5613
Intro->Graduation 0.3438 0.2886 0.4370 0.1484 1.3239
Medium->Graduation 0.6439 0.6843 0.5805 0.1038 1.0440
deg.fr p.value sig.05
HighSchool->Intro 350 0.1811 no
HighSchool->Medium 350 0.1537 no
HighSchool->Graduation 350 0.3357 no
Intro->Medium 350 0.2875 no
Intro->Graduation 350 0.0932 no
Medium->Graduation 350 0.1486 no
Inner models in the following objects:
$global
$group1
$group2
회귀식을 활용하는 방법
Qureshi & Compeau (2009)의 논문을 보면 PA의 결과와 나중에 소개할 MACS 방법을 상호 비교하는 것이 핵심입니다만... PA는 그저 group 간 어떤 차이가 있나에 집중하지, latent variable의 mean 값을 계산해주지는 않습니다. 문정훈 교수님께서 원하시는 것이 이것은 아니지 않나 생각합니다.
PLS의 weight로 linear combination을 수행하여 latent variable의 값을 계산하고 이것으로 t-test를 하는 방법을 소개합니다. 관측 변수의 scale이 모두 동일해야 합니다.
앞서 college 데이터를 FEMALE과 MALE로 나누겠습니다.
{dplyr} 패키지는 data manipulation을 도와주기 때문에 사용하겠습니다. 없다면 설치합시다.
R> install.packages("dplyr")
R> library(dplyr)
데이터를 쪼갭니다.
R> college_f <- college %>% filter(Gender=="FEMALE")
R> college_m <- college %>% filter(Gender=="MALE")
각각을 PLS로 fitting 시키겠습니다.
R> gpa_pls_f = plspm(college_f, gpa_path, gpa_blocks, modes = gpa_modes,
boot.val = TRUE)
R> gpa_pls_m = plspm(college_m, gpa_path, gpa_blocks, modes = gpa_modes,
boot.val = TRUE)
Outer model의 결과값을 따로 저장합시다.
R> outer_model_f <- (summary(gpa_pls_f))$outer_model
R> outer_model_m <- (summary(gpa_pls_m))$outer_model
설명을 위해 열어보겠습니다.
R> outer_model_f
weight 열은 비율로 정의되지 않기 때문에 x_i/sum(x) 형태로 변환되어야 하겠습니다. 기초 과목의 성적에 남녀의 차이를 확인하려 한다면 block이 Intro인 4,5,6,7번 데이터를 대상으로 latent variable의 regression prediction을 계산합니다.
R> Intro_p_f <- college_f %>% select(4:7)
R> Intro_lv_f <- apply(Intro_p_f,1,function(x) as.vector( t(x) %*% matrix(outer_model_f[4:7,3]/sum(outer_model_f[4:7,3]))))
마찬가지로 MALE의 경우에,
R> Intro_p_m <- college_m %>% select(4:7)
R> Intro_lv_m <- apply(Intro_p_m,1,function(x) as.vector( t(x) %*% matrix(outer_model_m[4:7,3]/sum(outer_model_m[4:7,3]))))
Welch sample t-test로 그룹 간 차이를 보겠습니다.
R> t.test(Intro_lv_f,Intro_lv_m)
Welch Two Sample t-test
data: Intro_lv_f and Intro_lv_m
t = -0.8492, df = 295.75, p-value = 0.3965
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.15742020 0.06251693
sample estimates:
mean of x mean of y
3.245116 3.292568
표준편차를 구해봅시다.
R> sd(Intro_lv_f)
[1] 0.5231552
R> sd(Intro_lv_m)
[1] 0.5023077
마찬가지로 그룹 간 차이가 없습니다.
Mean and Covariance Structure 접근법
Basic Concept
간단히 개념들을 소개합니다(Beaujean, 2014).
- Invariance. The general question of invariance is one of whether or not, under different conditions of observing and studying phenomena, measurements yield measures of the same attributes.
- Mean and Covariance Structure (MACS) question. A question that often arises when examining group differences on a variable is if any observed differencs are due to differences in the groups themselves or differences in the instrument(s) used to measure the variable. One of the more useful ways of answering this question is by assessing invariance in latent variables. In MACS, the indicator variables' intercepts are labeled with a mu and the latent variable's mean is referred as alpha as shown in the following figure.
- Measurement Invariance (MI). In empirical research, comparisons of means or regression coefficients is often drawn from distinct population groups such as culture, gender, language spoken. Unless explicitly tested, these analysis automatically assumes the measurement of these outcome variables are equivalent across these groups. Measurement invariance can be tested and it is important to make sure that the variables used in the analysis are indeed comparable constructs across distinct groups.
- MI - Configural level. This is the most basic level, and just indicates that Latent Variable Model(LVM)s have the same structure in all the groups.
- MI - Weak Invariance level. This level adds that, for a given indicator, the loadings are the same across groups.
- MI - Strong Invariance level. This level of invariance adds the constraint that, for a given indicator variable, the intercepts are the same among groups.
- MI - Strict Invariance level. This level adds that, for a given indicator, the error variances are constrained to be equal across groups.
MI 확인법
- Fit값을 비교해서 너무 지나치게 떨어지면 non-invariance 의심
- Modification index를 활용하는 방법
R을 활용하면 더 간단하게 볼 수 있습니다. 조금 있다가 다시 설명하겠습니다.
사례
R을 사용하겠습니다.
필요한 패키지를 일단 설치를 합시다.
R> install.packages("lavaan")
R> install.packages("semTools")
R> library(lavaan)
R> library(semTools)
R> library(semTools)
샘플 데이터는 lavaan에 포함된 HolzingerSwineford1939 를 사용합시다.
R> data(HolzingerSwineford1939)
기본 모형을 다음과 같이 정의하겠습니다. 예는 visual이라는 latent construct에 대해 x1, x2, 그리고 x3라는 observation이 있다고 하겠습니다.
latent 구조를 정해주는 formula는 =~ 기호를 사용합니다.
R> model <- 'visual=~x1+x2+x3'
먼저 MI를 확인합시다.
R> semTools::measurementInvariance(model,data=HolzingerSwineford1939,group="sex")
Measurement invariance models:
Model 1 : fit.configural
Model 2 : fit.loadings
Model 3 : fit.intercepts
Model 4 : fit.means
Chi Square Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit.configural 0 2723.4 2790.1 0.0000
fit.loadings 2 2726.6 2785.9 7.2103 7.2103 2 0.027184 *
fit.intercepts 4 2725.1 2777.0 9.7076 2.4973 2 0.286885
fit.means 5 2732.0 2780.2 18.6492 8.9415 1 0.002788 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Fit measures:
cfi rmsea cfi.delta rmsea.delta
1 1.000 0.000 NA NA
2 0.952 0.132 0.048 0.132
3 0.948 0.097 0.005 0.034
4 0.875 0.135 0.073 0.037
결과를 보면 모형 간의 chi-square test와 주요 fit indices의 변화를 함께 보여줍니다.
결과를 보면 strong level까지는 MI가 유지됩니다.
평균 비교는 partial invariance가 유지되는 수준에서도 가능하다고 하기 때문에(Byrne et al., 1989), 계속 latent variable의 group간 평균을 비교해보겠습니다.
본 예제는 Dimitrov (2006)의 절차를 따라 합니다. 우선 두 그룹의 평균을 invariant하다 고정하는 경우(model1), 그리고 한 그룹의 평균을 기준으로(평균=0), 다른 평균의 상대차이(mean difference)를 보는 경우(model2)를 각각 생성한 다음, 두 모형의 chi-square 차이를 봐서 평균 차이의 significance를 확인합니다.
R> model.equl.means = cfa(model,data=HolzingerSwineford1939,group="sex",
group.equal=c("loadings","intercepts","residuals",
"lv.variances","means"))
R> model.unequl.means = cfa(model,data=HolzingerSwineford1939,group="sex",
group.equal=c("loadings","intercepts","residuals",
"lv.variances"))
Confirmatory Factor Analysis (CFA) 모형을 작성하는 함수인 cfa() 안에 모형과, 데이터, 그리고 그룹 결정 변수의 이름을 주고, group 동질성의 가능 수준을 적습니다. 우선, model.equal.means는 means까지 통제하지만, model.unequal.means는 그냥 둡니다.
group.equal 파라미터를 정리하면 다음과 같습니다.
[1] intercept 기울기 같게
[2] means 평균 같게
[3] residuals 잔차 구조 동일
[4] residual.covariances 공분산 구조 동일
[5] lv.variances 잠재변수의 잔차 구조 동일
[6] lv.covariances 잠재변수의 공분산 구조 동일
[7] regressions 회귀계수 동일
본 예에서 본 것처럼 cfa에서 covariance 구조만 허락하고 나머지는 다 묶었습니다.
우선,
R> summary(model.equl.means,fit.measures=T,standardized=T)
lavaan (0.5-18) converged normally after 29 iterations
Number of observations per group
1 146
2 155
Estimator ML
Minimum Function Test Statistic 26.562
Degrees of freedom 9
P-value (Chi-square) 0.002
Chi-square for each group:
1 12.057
2 14.505
Model test baseline model:
Minimum Function Test Statistic 115.121
Degrees of freedom 6
P-value 0.000
User model versus baseline model:
Comparative Fit Index (CFI) 0.839
Tucker-Lewis Index (TLI) 0.893
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1356.977
Loglikelihood unrestricted model (H1) -1343.696
Number of free parameters 9
Akaike (AIC) 2731.955
Bayesian (BIC) 2765.319
Sample-size adjusted Bayesian (BIC) 2736.776
Root Mean Square Error of Approximation:
RMSEA 0.114
90 Percent Confidence Interval 0.065 0.165
P-value RMSEA <= 0.05 0.019
Standardized Root Mean Square Residual:
SRMR 0.115
Parameter estimates:
Information Expected
Standard Errors Standard
Group 1 [1]:
Estimate Std.err Z-value P(>|z|) Std.lv Std.all
Latent variables:
visual =~
x1 1.000 0.724 0.621
x2 (.p2.) 0.778 0.141 5.532 0.000 0.563 0.479
x3 (.p3.) 1.107 0.214 5.173 0.000 0.801 0.710
Intercepts:
x1 (.p8.) 4.936 0.067 73.473 0.000 4.936 4.235
x2 (.p9.) 6.088 0.068 89.855 0.000 6.088 5.179
x3 (.10.) 2.250 0.065 34.579 0.000 2.250 1.993
visual 0.000 0.000 0.000
Variances:
x1 (.p4.) 0.835 0.118 0.835 0.614
x2 (.p5.) 1.065 0.105 1.065 0.771
x3 (.p6.) 0.633 0.129 0.633 0.496
visual (.p7.) 0.524 0.130 1.000 1.000
Group 2 [2]:
Estimate Std.err Z-value P(>|z|) Std.lv Std.all
Latent variables:
visual =~
x1 1.000 0.724 0.621
x2 (.p2.) 0.778 0.141 5.532 0.000 0.563 0.479
x3 (.p3.) 1.107 0.214 5.173 0.000 0.801 0.710
Intercepts:
x1 (.p8.) 4.936 0.067 73.473 0.000 4.936 4.235
x2 (.p9.) 6.088 0.068 89.855 0.000 6.088 5.179
x3 (.10.) 2.250 0.065 34.579 0.000 2.250 1.993
visual 0.000 0.000 0.000
Variances:
x1 (.p4.) 0.835 0.118 0.835 0.614
x2 (.p5.) 1.065 0.105 1.065 0.771
x3 (.p6.) 0.633 0.129 0.633 0.496
visual (.p7.) 0.524 0.130 1.000 1.000
결과에서 볼 수 있듯
visual의 mean은 0으로 고정입니다(이것이 null model).
이제 unequal mean를 가정한 것을 보면,
R> summary(model.unequl.means,fit.measures=T,standardized=T)
lavaan (0.5-18) converged normally after 31 iterations
Number of observations per group
1 146
2 155
Estimator ML
Minimum Function Test Statistic 17.078
Degrees of freedom 8
P-value (Chi-square) 0.029
Chi-square for each group:
1 6.483
2 10.594
Model test baseline model:
Minimum Function Test Statistic 115.121
Degrees of freedom 6
P-value 0.000
User model versus baseline model:
Comparative Fit Index (CFI) 0.917
Tucker-Lewis Index (TLI) 0.938
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1352.235
Loglikelihood unrestricted model (H1) -1343.696
Number of free parameters 10
Akaike (AIC) 2724.470
Bayesian (BIC) 2761.541
Sample-size adjusted Bayesian (BIC) 2729.827
Root Mean Square Error of Approximation:
RMSEA 0.087
90 Percent Confidence Interval 0.026 0.144
P-value RMSEA <= 0.05 0.128
Standardized Root Mean Square Residual:
SRMR 0.084
Parameter estimates:
Information Expected
Standard Errors Standard
Group 1 [1]:
Estimate Std.err Z-value P(>|z|) Std.lv Std.all
Latent variables:
visual =~
x1 1.000 0.678 0.587
x2 (.p2.) 0.804 0.142 5.645 0.000 0.545 0.466
x3 (.p3.) 1.198 0.222 5.392 0.000 0.812 0.729
Intercepts:
x1 (.p8.) 5.091 0.085 59.773 0.000 5.091 4.405
x2 (.p9.) 6.213 0.080 77.378 0.000 6.213 5.314
x3 (.10.) 2.437 0.088 27.547 0.000 2.437 2.186
visual 0.000 0.000 0.000
Variances:
x1 (.p4.) 0.876 0.111 0.876 0.656
x2 (.p5.) 1.070 0.104 1.070 0.783
x3 (.p6.) 0.582 0.129 0.582 0.469
visual (.p7.) 0.460 0.114 1.000 1.000
Group 2 [2]:
Estimate Std.err Z-value P(>|z|) Std.lv Std.all
Latent variables:
visual =~
x1 1.000 0.678 0.587
x2 (.p2.) 0.804 0.142 5.645 0.000 0.545 0.466
x3 (.p3.) 1.198 0.222 5.392 0.000 0.812 0.729
Intercepts:
x1 (.p8.) 5.091 0.085 59.773 0.000 5.091 4.405
x2 (.p9.) 6.213 0.080 77.378 0.000 6.213 5.314
x3 (.10.) 2.437 0.088 27.547 0.000 2.437 2.186
visual -0.302 0.103 -2.929 0.003 -0.445 -0.445
Variances:
x1 (.p4.) 0.876 0.111 0.876 0.656
x2 (.p5.) 1.070 0.104 1.070 0.783
x3 (.p6.) 0.582 0.129 0.582 0.469
visual (.p7.) 0.460 0.114 1.000 1.000
sex=1을 기준으로 sex=2의 차이를 보면 visual이 -0.302입니다(0.3정만큼 상대적으로 평균값이 더 작습니다). 이것이 유의미한지 알아보려면,
R> anova(model.equl.means,model.unequl.means)
Chi Square Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
model.unequl.means 8 2724.5 2761.5 17.078
model.equl.means 9 2731.9 2765.3 26.562 9.4844 1 0.002072 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
유의미한 p-value가 나오기 때문에 평균 차이가 유의미합니다.
참고문헌
- Beaujean, A. A. (2014). Latent variable modeling using R: A step-by-step guide. Routledge.
- Byrne, B. M., Shavelson, R. J., & Muthén, B. (1989). Testing for the equivalence of factor covariance and mean structures: The issue of partial measurement invariance. Psychological Bulletin, 105(3), 456-466.
- Dimitrov, D. M. (2006). Comparing groups on latent variables: A structural equation modeling approach. Work, 26(4), 429-436.
- Sanchez, G. (2013). PLS Path Modeling with R, Berkeley, http://gastonsanchez.com/PLS_Path_Modeling_with_R.pdf
- Qureshi, I., & Compeau, D. (2009). Assessing between-group differences in information systems research: a comparison of covariance-and component-based SEM. MIS Quarterly, 33(1), 197-214.
Appendix
변수의 가중치를 반영한 latent variable을 생성하는 R 코드를 남겨둡니다.파일: create_lv_by_pls.R
## ============================= ##
## CREATE A LATENT VARIABLE ##
## USING PARTIAL LEAST SQUARE ##
## TAEKYUNG KIM ##
## 2015.11. ##
## ============================= ##
## HOW TO RUN ##
# Let's assume that you have a data set including 3 observation variables.
# Suppose that those are indexed as 3,4,5.
# In order to calculate the latent variable of them,
# you should put data and index into "create_lv_pls" as parameters.
# Such that,
# R> lv_new <- create_lv_pls(your_data_here, 3:5)
if(!require(plspm)) {
install.packages("plspm")
library(plspm)
}
if(!require(dplyr)) {
install.packages("dplyr")
library(dplyr)
}
create_lv_pls <- function(data,index) {
options(warn=-1)
test_data0 <- data %>% select(index)
test_data <- cbind(test_data0,test_data0)
names(test_data) <- paste("ob",1:(length(index)*2))
X <- c(0,0)
Y <- c(1,0)
path = rbind(X,Y)
blocks = list(1:length(index),(1:length(index))+length(index))
modes = c("A","A")
pls.result = plspm(test_data,path,blocks,modes=modes,boot.val=TRUE)
result = summary(pls.result)
items_weight <- result$outer_model[1:length(index),3]
lv = apply(test_data0,1,function(x) as.vector(t(x) %*% matrix(items_weight/sum(items_weight))))
options(warn=0)
lv
}
이제 실행할 때, 같은 경로에 위의 코드를 가진 파일 create_by_pls.R을 두고
R> source("create_by_pls.R")
R> your_index = c(1,2,3) #만약 1,2,3번에 latent variable을 구성하는 관측변수라면...
R> new_latent_variable = create_lv_pls(your_data, your_index)
R> mean(new_latent_variable)
R> sd(new_latent_variable)
이상입니다.
댓글
댓글 쓰기