
############################
dvar <- read.csv("http://commres.net/wiki/_media/elemapi2_.csv", sep="\t", fileEncoding="UTF-8-BOM")
dvar
mod <- lm(api00 ~ ell + acs_k3 + avg_ed + meals, data=dvar)
summary(mod)
anova(mod)

mod01 <- lm(api00 ~ acs_k3 + ell + avg_ed + meals, data=dvar)
summary(mod01)
anova(mod01)

mod02 <- lm(api00 ~ ell + avg_ed + acs_k3 + meals, data=dvar)
summary(mod02)
anova(mod02)

mod03 <- lm(api00 ~ ell + avg_ed + meals + acs_k3, data=dvar)
summary(mod03)
anova(mod03)

#################################################
#setwd("user location") #Working directory
set.seed(123) #Standardizes the numbers generated by rnorm; see Chapter 5
n <- 100 #Number of participants; graduate students
x <- rnorm(N, 175, 7) #IV; hours since dawn
m <- 0.7*X + rnorm(N, 0, 5) #Suspected mediator; coffee consumption 
y <- 0.4*M + rnorm(N, 0, 5) #DV; wakefulness
md <- data.frame(x, m, y)

lm.yx <- lm(y~x,data=md)
summary(lm.yx)

lm.mx <- lm(m~x, data=md)
summary(lm.mx)

lm.ymx <- lm(y~m+x, data=md)
summary(lm.ymx)

#####################################
# x와 m이 상관관계가 높고 
# 이론적으로 혹은 논리적으로 
# 인과관계가 분명할 때,
# x의 (y로 가는) 영향력이
# m으로 옮겨가서 영향력을 
# 간접적으로 주는 것으로 해석할 
# 수 있다. 이 때 
# mediation effect를 계산한다
#####################################

lm.xm <- lm(x~m, data = md)
summary(lm.xm)
lm.xm.res <- lm.xm$residuals
lm.yxmres <- lm(y~lm.xm.res, data=md)
summary(lm.yxmres)

lm.reg <- lm.xm$fitted.values
lm.yreg <- lm(y~lm.reg,data=md)
summary(lm.yreg)

lm.reg2 <- lm.mx$fitted.values
lm.yreg2 <- lm(y~lm.reg2,data=md)
summary(lm.yreg2)

cor(lm.reg, lm.reg2)

#####################################
attach(md)
lma <- lm(x~m, data=md)
lmb <- lm(m~x, data=md)
summary(lma)
summary(lmb)



reg <- lma$fitted.values
yreg <- lm(y~reg, data = md)
summary(yreg)
res <- lma$residuals
yres <- lm(y~res, data=md)
summary(yres)

ss.x <- var(x)*(length(x)-1)
ss.m <- var(m)*(length(m)-1)
var(x)
ss.x
ss.reg <- sum(reg^2)
rsa <- ss.reg/ss.x
rsa

reg2 <- lmb$fitted.values
var(m)
ss.m
ss.reg2 <- sum(reg2^2)
rsb <- ss.reg2/ss.m
rsb


###########################
options(digits = 4)

HSGPA <- c(3.0, 3.2, 2.8, 2.5, 3.2, 3.8, 3.9, 3.8, 3.5, 3.1)
FGPA <-  c(2.8, 3.0, 2.8, 2.2, 3.3, 3.3, 3.5, 3.7, 3.4, 2.9)
SATV <-  c(500, 550, 450, 400, 600, 650, 700, 550, 650, 550)
GREV <-  c(600, 670, 540, 800, 750, 820, 830, 670, 690, 600)
##GREV <- c(510, 670, 440, 800, 750, 420, 830, 470, 690, 600)

scholar <- data.frame(HSGPA, FGPA, SATV, GREV) # collect into a data frame
describe(scholar) # provides descrptive information about each variable

attach(scholar)

m00 <- lm(FGPA ~ SATV, data=scholar)
summary(m00)
m01 <- lm(HSGPA ~ SATV, data=scholar)
summary(m01)
m03 <- lm(FGPA ~ HSGPA + SATV, data=scholar)
summary(m03)

pcor.test(FGPA, SATV, HSGPA)
pcor.test(FGPA, HSGPA, SATV)

# proof of other variables controlled
###########################

a <- spcor.test(FGPA, SATV, HSGPA)$estimate^2
b <- spcor.test(FGPA, HSGPA, SATV)$estimate^2
c <- cor(FGPA, SATV)^2
d <- cor(FGPA, HSGPA)^2
c-a
d-b
# common affected area proof
###########################