User Tools

Site Tools


mediation_analysis
# 개인과제와 관련된 코드입니다. 
# Carseats 데이터 분석입니다. 새로 시작하면 
library(ISLR)
# 만약에 ISLR이 install도 안되어 
# 있으면 
# install.packages("ISLR")

?Carseats   # 데이터 설명 확인
str(Carseats)


# 이 중에서 CompPrice의 독립변인으로서의
# 역할에 문제점이 보이는 듯 하여 물어봅니다
# CompPrice와 다른 변인들을 모두 활용해도 
# 되겠지만 간단하게 보기 위해서 Price만을
# 독립변인으로 분석을 합니다

lm.c1 <- lm(Sales ~ Price + CompPrice, data = Carseats)
summary(lm.c1)

# output 을 살펴보면 R 제곱값이 
# 0.3578 (35.78%) 임을 알 수 있습니다. 
# 이는 우리가 흔히 쓰는 다이어그램에서
#
# http://commres.net/wiki/_detail/pasted/20201201-170048.png?id=multiple_regression_exercise
#
# b + c + d 에 해당하는 부분입니다. 즉 두 변인의
# 설명력을 보여주는 부분인 lm.c1에서의 
# R 제곱은 b + c + d / a + b + c + d 입니다. 

# 여기에서 아래를 수행하여 
# semipartial corrleation 값과 이를 
# 제곱한 값을 알아봅니다. 

library(ppcor)
attach(Carseats)
spcor.Price <- spcor.test(Sales, Price, CompPrice)
spcor.CompPrice <- spcor.test(Sales, CompPrice, Price)

spcor.Price
spcor.CompPrice

# 위 둘의 아웃풋에서 estimate값은 
# semipartial correlation 값이므로 이를 
# 제곱한 값은 각각 b와 d에 해당하는 
# 값입니다. 이를 더 자세히 이야기하면
# b = b / a + b + c + d
# d = d / a + b + c + d
# 라는 이야기입니다. 

b <- spcor.Price$estimate^2
d <- spcor.CompPrice$estimate^2 

# 각 b와 d를 출력해 봅니다. 
b
d

# 이 둘을 더해봅니다. 
# 이 값은 0.5135791 입니다
b + d

# 다시 아까 lm.c1의 R 제곱값은 
summary(lm.c1)$r.squared
# 0.3578332 입니다. 

# 다시 그림을 보면
# R 제곱은 b + c + d 에 해당하는 
# 것이라고 했고
# semi-partial 값을 각각 구하여
# 이를 제곱하여 더한 것이 b와 d입니다.
# 그런데 그림을 보면 b + d > b + c + d 가
# 됩니다.  그렇지만 그림대로라면 
# 이것은 말이 되지 않습니다. 그런데 
# 계산값을 이를 가르키고 있습니다. 
# 이 이유가 뭘까요?

# 여러분의 생각 후에 mediator 분석에 대해서
# 설명을 하도록 하겠습니다. 
#
#

cs.dat <- Carseats[, c(1,2,6)]
str(cs.dat)
round(cor(cs.dat),3)
ia.1 <- lm(Sales~CompPrice, data=cs.dat)
ia.2 <- lm(Sales~Price, data=cs.dat)
ia.3 <- lm(Sales~CompPrice+Price, data=cs.dat)

# Price 에 CompPrice를 더한 것에 따라서 
# r square 값의 증가분이 의미가 있는가를 
# 테스트 
anova(ia.2, ia.3)
# 증가분은 
summary(ia.3)$r.squared 
summary(ia.2)$r.squared
summary(ia.3)$r.squared - summary(ia.2)$r.squared

# mediated effect = CompPrice에 대한 Price의 설명력이
# 얼마나 Sales를 설명하는가?

ia.0 <- lm(CompPrice~Price, data=cs.dat)

# CompPrice에 대한 Price의 설명력
med.price.2.compprice <- ia.0$fitted.values - mean(cs.dat$CompPrice)
cor(med.price.2.compprice, cs.dat$Sales)
cor(cs.dat$Price, cs.dat$Sales)

# !!!! ????
# Are they the same?
cor(med.price.2.compprice, cs.dat$Price)
# examine the head of each
head(med.price.2.compprice)
head(cs.dat$Price)

ia.01 <- lm(Sales~med.price.2.compprice, data=cs.dat)
summary(ia.01)
summary(ia.2)

# cor between residuals and regressions CompPrice 
cor(ia.0$residuals, cs.dat$Sales)
cor(med.price.2.compprice, cs.dat$Sales)
cor.test(ia.0$residuals, cs.dat$Sales)
cor.test(med.price.2.compprice, cs.dat$Sales)
cor(ia.0$residuals+med.price.2.compprice, cs.dat$Sales)
cor.test(ia.0$residuals+med.price.2.compprice, cs.dat$Sales)

ia.02 <- lm(Sales~ia.0$residuals, data=cs.dat)
summary(ia.02)

summary(ia.01)$r.squared
summary(ia.02)$r.squared
summary(ia.01)$r.squared + summary(ia.02)$r.squared
summary(ia.3)


#########################################
## data file: PlannedBehavior.csv
#########################################
df <- read.csv("http://commres.net/wiki/_media/r/plannedbehavior.csv")
head(df)
str(df)
# attitude
# norms
# control
# intention
# behavior

pb.m01 <- lm(intention~attitude+norms+control, data=df)
summary(pb.m01)
reg.intention <- pb.m01$fitted.values - mean(df$intention)
head(reg.intention)
pb.m0x <- lm(behavior~reg.intention, data=df)
summary(pb.m0x)
pb.m0y <- lm(behavior~intention, data = df)
summary(pb.m0y)
pb.m0z <- lm(behavior~pb.m01$residuals, data=df)
summary(pb.m0z)

summary(pb.m0x)$r.squared # regression part의 설명력
summary(pb.m0z)$r.squared # residuals part의 설명력
# 둘을 더한 아래값을 
summary(pb.m0x)$r.squared + summary(pb.m0z)$r.squared
# behavior~intention regression r squared 값과 비교
# intention 전체가 behavior에 (variance에) 기여한 부분
summary(pb.m0y)$r.squared 
# 0.1983717 ~= 0.1982297

############################
# 2nd 
############################
# using only attitude

lm.pb021 <- lm(intention~attitude, data=df)
lm.pb022 <- lm(behavior~attitude, data=df)
lm.pb023 <- lm(behavior~attitude+intention, data=df)

reg.intention <- summary(lm.pb021)$fitted.values - mean(df$intention)
res.intention <- summary(lm.pb021)$residuals

# what if the other way around?
lm.pb024 <- lm(attitude~intention, data=df)
summary(lm.pb024)
summary(lm.pb021)
# note that r square values are the same
summary(lm.pb024)$r.squared
summary(lm.pb021)$r.squared
summary(lm.pb024)$r.squared == summary(lm.pb021)$r.squared

# the effect of residuals of attitude
res.attitude <- summary(lm.pb024)$residuals
reg.attitude <- summary(lm.pb024)$fitted.values - mean(df$attitude)

lm.pb025 <- lm(behavior~res.attitude, data=df)
summary(lm.pb025)
# for the reference the whole part of attitude
# explains behavior significantly
summary(lm.pb022)
# summary(lm(behavior~attitude, data=df))
# 그러나 위의 분석으로는 
# attitude는 intention을 설명하는 부분을 제외하고는 
# behavior 설명에 기여하는 부분이 작다
summary(lm.pb023)

###########################
###########################
###########################
# Mediation Analysis using lavaan 
# in R (path analysis framework, SEM)
###########################
###########################
# install,packages("lavann")
library(lavaan)

specmod <- "
    # path c
    behavior ~ c*attitude

    # path a
    intention ~ a*attitude

    # path b 
    behavior ~ b*intention
    
    # indirect effect (a*b): Sobel test (Delta Method)
    ab := a*b
"
# Fit/estimate the model
fitmod <- sem(specmod, data=df)

# summarize the model
summary(fitmod, fit.measures=TRUE, rsquare=T)
lm.ba.01 <- lm(behavior~attitude+intention, data=df)
lm.ba.02 <- lm(behavior~intention, data=df)
lm.ba.03 <- lm(intention~attitude, data=df)
lm.ba.04 <- lm(attitude~intention, data=df)
lm.ba.05 <- lm(behavior~attitude, data=df)

reg.att <- lm.ba.03$fitted.values - mean(df$intention)
res.att <- summary(lm.ba.03)$residuals
lm.ba.021 <- lm(behavior~reg.att, data=df)
summary(lm.ba.021)
lm.ba.022 <- lm(behavior~res.att, data=df)
summary(lm.ba.022)

summary(lm.ba.02)$r.squared 
summary(lm.ba.021)$r.squared 
summary(lm.ba.022)$r.squared
summary(lm.ba.021)$r.squared + summary(lm.ba.022)$r.squared

summary(lm.ba.021)
summary(lm.ba.05)
cor(reg.att, df$attitude)
mediation_analysis.txt · Last modified: 2021/12/09 09:35 by hkimscil