########################################
#Cohen, J., Cohen, A., West, S., & Aiken, L. (2003).
#Interactions Among Continuous Variables.
#In Applied multiple regression/correlation analysis for the behavioral sciences(3rd ed., pp.255--301)
#Mahwah, NJ: Erlbaum.
#######################################

##(a) データの読み込み
dat1 <- read.table("http://blue.zero.jp/yokumura/R/reg/C07E01DT.txt",header=F)
names(dat1) <- c("id","x","z","y")


##(b) センタリング
dat2 <- dat1
names(dat2) <- c("x", "z", "xz", "y")
dat2$x      <- dat1$x - mean(dat1$x)
dat2$z      <- dat1$z - mean(dat1$z)
dat2$xz     <- dat2$x * dat2$z


##(c) Summary Statistics for centered x and z and cross-product xz (Table 7.4.1 A)
digits <- 2
round(mean(dat2), digits)             #Mean
round(sd(dat2),   digits)             #SD
round(cor(dat2),  digits)             #Correlation


##(d) Centered regression equations (Table 7.4.1 B)
reg1 <- lm(y~x+z, data=dat2)
reg2 <- lm(y~x*z, data=dat2)
summary(reg1)
summary(reg2)
anova(reg1, reg2)           #F_gain(1, 241 = 12.08) p < .01


##(e) Regression of endurance on age at three levels of exercise (Figure 7.4.1 A)
z.sd <- sd(dat2$z)
coef.high <- c(coef(reg2)[["z"]] * z.sd + coef(reg2)[["(Intercept)"]],
               coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * z.sd)
coef.mean <- c(coef(reg2)[["z"]] * 0 + coef(reg2)[["(Intercept)"]],
               coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * 0)
coef.low  <- c(coef(reg2)[["z"]] * -z.sd + coef(reg2)[["(Intercept)"]],
               coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * -z.sd)

x.sd <- sd(dat2$x)
plot(y~x, data=dat2, xlim=c(-x.sd*1.2, x.sd*1.2), type="n", xaxt="n",
                     xlab="Age centered", ylab="Endurance (minutes)")
axis(side=1, at=c(-x.sd, 0, x.sd), labels=round(c(-x.sd, 0, x.sd), 2))

abline(coef.high, col="red")
abline(coef.mean)
abline(coef.low,  col="blue")

graphics.off()


##(f) Covariance matrix of regression coefficients in the centered regression equation (Table 7.4.1 C)
Sb <- summary(reg2)$sigma^2 * summary(reg2)$cov.unscaled[-1,-1]
round(Sb, 7)


##(g) Analysis of simple regression equations (Table 7.4.1 D)
sig.level <- .05

#
SEb.at.z.high <- sqrt(Sb[1,1] + 2 * z.sd * Sb[1,3] + z.sd^2 * Sb[3,3])                 #equation 7.4.2
t.b.at.z.high <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * z.sd) / SEb.at.z.high      #equation 7.4.4
p.high        <- pt(abs(t.b.at.z.high), summary(reg2)$df[2], lower.tail=F)

me <- qt(1-sig.level/2, summary(reg2)$df[2]) * SEb.at.z.high      #equation 7.4.6
lower.high <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * z.sd) - me
upper.high <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * z.sd) + me


#
SEb.at.z.mean <- sqrt(Sb[1,1] + 2 * 0 * Sb[1,3] + 0^2 * Sb[3,3])                 #equation 7.4.2
t.b.at.z.mean <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * 0) / SEb.at.z.mean      #equation 7.4.4
p.mean        <- pt(abs(t.b.at.z.mean), summary(reg2)$df[2], lower.tail=F)

me <- qt(1-sig.level/2, summary(reg2)$df[2]) * SEb.at.z.mean      #equation 7.4.6
lower.mean <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * 0) - me
upper.mean <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * 0) + me


#
SEb.at.z.low <- sqrt(Sb[1,1] + 2 * -z.sd * Sb[1,3] + z.sd^2 * Sb[3,3])                 #equation 7.4.2
t.b.at.z.low <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * -z.sd) / SEb.at.z.low      #equation 7.4.4
p.low        <- pt(abs(t.b.at.z.low), summary(reg2)$df[2], lower.tail=F)

me <- qt(1-sig.level/2, summary(reg2)$df[2]) * SEb.at.z.low      #equation 7.4.6
lower.low <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * -z.sd) - me
upper.low <- (coef(reg2)[["x"]] + coef(reg2)[["x:z"]] * -z.sd) + me


#
output <- matrix(0, ncol=6, nrow=3)
output <- data.frame(output)
colnames(output) <- c("value", "simple", "std", "t.value", "p.value", "conf.int")
rownames(output) <- c("low", "mean", "high")

output$value   <- round(c(-z.sd, 0, z.sd),2)
output$simple  <- c(paste(round(coef.low[2], 3), "x", " + ", round(coef.low[1], 2), sep=""),
                   paste(round(coef.mean[2], 3), "x", " + ", round(coef.mean[1], 2), sep=""),
                   paste(round(coef.high[2], 3), "x", " + ", round(coef.high[1], 2), sep=""))
output$std     <- round(c(SEb.at.z.low, SEb.at.z.mean, SEb.at.z.high),3)
output$t.value <- round(c(t.b.at.z.low, t.b.at.z.mean, t.b.at.z.high),2)
output$p.value <- round(c(p.low, p.mean, p.high), 3)
output$conf.int <- c(paste(round(c(lower.low, upper.low), 2), collapse="/"),
                    paste(round(c(lower.mean, upper.mean), 2), collapse="/"),
                    paste(round(c(lower.high, upper.high), 2), collapse="/"))
output


##(h) Standardized solution for the regression (Table 7.5.1 B)
dat3 <- dat2
dat3$x      <- scale(dat1$x)
dat3$z      <- scale(dat1$z)
dat3$xz     <- dat3$x * dat3$z
dat3$y      <- scale(dat1$y)

round(mean(dat3), digits)             #Mean
round(sd(dat3),   digits)             #SD
round(cor(dat3),  digits)             #Correlation

summary(lm(y~x*z, data=dat3))