# Options

# Surpress scientific notations, digits and lines

options(scipen=999, digits=10, max.print=99999999)


# Correlation

cor(Data.01.PPT)

cor(Data.01.PPT[,13:16])

library(Hmisc)

rcorr(as.matrix(Data.01.PPT[,13:16]))

library(corrgram)

corrgram(Data.01.PPT)


# Regression Analysis

library(car)

R00<-lm(ATT ~ PEU + PU, data=Data.01.PPT)

summary(R00)

vif(R00)

plot(R00)


# Standardized Coefficients

library(QuantPsyc)
lm.beta(R00)

# Alternatively

R00$coefficients[["PEU"]]*(sd(Data.01.PPT$PEU)/sd(Data.01.PPT$ATT))

R00$coefficients[["PU"]]*(sd(Data.01.PPT$PU)/sd(Data.01.PPT$ATT))

# Bootstrapping

bs<-function(formula, data, indices){
d<-data[indices,]
m<-lm(formula, data=d)
return(coef(m))
}

boot.R00<-boot(data=Data.01.PPT, statistic=bs, R=1000, formula=ATT ~ PEU + PU)
summary(boot.R00)

boot.ci(boot.R00, type="bca", index=1)


# Cook's D plot
# identify D values > 4/(n-k-1) 
cutoff <- 4/((nrow(Data.01.PPT)-length(R00$coefficients)-2)) 
plot(R00, which=4, cook.levels=cutoff)


# Equality of Coefficients

R01<-lm(ATT ~ I(PEU + PU), data=Data.01.PPT)
summary(R01)

anova(R01,R00)

# Alternative Using car Package

library(car)

linearHypothesis(R00, "PEU=PU")

linearHypothesis(R00, "PEU=2*PU")


# Dummy variables

A.01<-lm(sales ~ factor(promotion), data=Table16.2)
summary(A.01)

summary(lm(lm(sales ~ promotion, data=Table16.2)))

Table16.2$prom.1<-factor(Table16.2$promotion)
Table16.2$prom.2<-factor(Table16.2$promotion)

contrasts(Table16.2$prom.1)<-contr.treatment(3, base=3)
contrasts(Table16.2$prom.2)<-contr.sum(3)

summary(lm(sales ~ prom.1, data=Table16.2))

summary(lm(sales ~ prom.2, data=Table16.2))


# Dummy variables

Data.01.PPT$Age.Cat<-as.factor(sample.int(3, size=439, replace=TRUE, prob=c(0.3,0.3,0.4)))
table(Data.01.PPT$Age.Cat)

contrasts(Data.01.PPT$Age.Cat)<-contr.treatment(3)
contrasts(Data.01.PPT$Age.Cat)<-contr.sum(3)

# base = 1

R00.Age<-lm(ATT ~ PEU + PU + Age.Cat, data=Data.01.PPT)
summary(R00.Age)

anova(R00.Age, R00)


# Mediation

options(scipen=999)

library(mediation)

MM<-lm(M ~ X, data=Mediation)
MY<-lm(Y ~ X + M, data=Mediation)
MED1<-mediate(MM,MY, treat='X', mediator='M')
summary(MED1)

library(bda)

mediation.test(Mediation$M, Mediation$X, Mediation$Y)

library(QuantPsyc)

Data.Mediation<-Mediation
names(Data.Mediation)<-c('y','x','m')

proximal.med(Data.Mediation)

library(boot)
med.boot <- boot(Data.Mediation, proxInd.ef, R=1000)
sort(med.boot$t)[c(25,975)] #95% CI
plot(density(med.boot$t)) # Distribution of bootstapped indirect effect
summary(med.boot$t)

boot.ci(med.boot, type="bca", index=1)


library(psych)
MED.1<-mediate(Y ~ X + (M), data=Mediation, n.iter=0)
mediate.diagram(MED.1, digits=10)
print(MED.1, digits=10, short=FALSE)

# Moderation

summary(lm(Y ~ X*Z, data=Moderation))

library(QuantPsyc)

MOD1<-moderate.lm(X, Z, Y, Moderation, mc=FALSE)
summary(MOD1)

GMOD1<-sim.slopes(MOD1, Moderation$Z, zsd=2, mcz=FALSE)

graph.mod(GMOD1, X, Y, Moderation)


library(psych)
MOD2<-mediate(Y ~ X*Z, data=Moderation)
print(MOD2, digits=10)



# Surface Plot

MCX <- seq(-4,4, 0.1) # Mean-centered
MCZ <- seq(-4,4, 0.1) # Mean-centered
f<-function(MCX,MCZ) {3.69259+0.80104*MCX+0.43046*MCZ+0.07861*MCX*MCZ}
Y <- outer(MCX, MCZ, f)
#z[is.na(z)] <- 1
#op <- par(bg = "white")
persp(MCX, MCZ, Y, theta = 20, phi = 20, expand = 0.5, col = "orange")



# Rockchalk

library(rockchalk)

R.rc<-lm(Y ~ X*Z, data=Moderation)

PS.rc<-plotSlopes(R.rc, plotx="X", modx="Z", modxVals="std.dev")

TS.rc<-testSlopes(PS.rc)

TS.rc$hypotests

TS.rc$jn

print(TS.rc)

plot(TS.rc)


# jtools

library(jtools)

MOD3<-lm(Y ~ X*Z, data=Moderation)

johnson_neyman(model = MOD3, pred = X, modx = Z)


# Spiller

BMI$LargeQuantModelf<-factor(BMI$LargeQuantModel)

M.Spiller<-lm(Candies ~ BMI*LargeQuantModel, data=BMI)
summary(M.Spiller)

M.Spillerf<-lm(Candies ~ BMI*LargeQuantModelf, data=BMI)
summary(M.Spiller)

PS.Spiller<-plotSlopes(M.Spiller, plotx="LargeQuantModel", modx="BMI", modxVals=c(0,1), interval="conf")

print(PS.Spiller)

TS.Spiller<-testSlopes(PS.Spiller)

plot(TS.Spiller)



PS1.Spiller<-plotSlopes(M.Spillerf, plotx="BMI", modx="LargeQuantModelf", interval="conf")

print(PS1.Spiller)

TS.Spiller<-testSlopes(PS1.Spiller)




# jtools

library(jtools)

johnson_neyman(model = M.Spiller, pred = LargeQuantModel, modx = BMI)



# pwr Regression

install.packages("pwr")

library(pwr)

# R2=0.30
# u=3 (independent variables)
# n=v+u+1

pwr.f2.test(u=3, f2=0.30/(1-0.30), sig.level=0.05, power=0.80)