# Homework 13 key

# 1. read in, describe by condition combos
d = dfReadDat('HW13_Data.dat', SubID=NULL)
varDescribeBy(d$Y, list(d$X, d$G))

# 2. example
# I'm just going to keep the variables as-is for the key

# 3. long to wide, renames
dW = reshape(d, timevar = 'X', idvar = c('SubID','G'), direction='wide')
258 / 3 # = 86, so this makes sense.

# 4. overall effect
dW$YM = rowMeans(dW[,c('Y.1','Y.2','Y.3')])
m1 = lm(YM ~ G, data=dW)
modelSummary(m1)
modelEffectSizes(m1)
# the group difference in significant. Those with G = 2 have Y scores 3.3 points higher, on average

# 5. pattern of means hypothesis
# 3 > 2 and 1
m2 = lm(Y.3 - (Y.2+Y.1)/2 ~ G, data=dW)
modelSummary(m2)
# This pattern of means is consistent with the hypothesized relationship, t(84) = 2.63, p = .01.
# This pattern is moderated by G, t(84) = 3.45, p < .001.

# 6. pairwise comparisons
# 3 > 2 > 1. I'll do all three comparisons
m3 = lm(Y.3 - Y.2 ~ G, data=dW)
modelSummary(m3) # G moderates the difference between Y3 and Y2. there is no overall difference
m4 = lm(Y.2 - Y.1 ~ G, data=dW)
modelSummary(m4) # G does not moderate the difference between Y2 and Y1. There is an overall difference though
m5 = lm(Y.3 - Y.1 ~ G, data=dW)
modelSummary(m5) # G moderates the difference between Y3 and Y1 AND there is an overall difference.

# CORRECTIONS
ps = c(0,0,0)

ps[1] = modelSummary(m3)$coefficients[2,4]  
ps[2] = modelSummary(m4)$coefficients[2,4]  
ps[3] = modelSummary(m5)$coefficients[2,4]  

ps = sort(ps)
ps

p.adjust(ps, "holm")
round(p.adjust(ps, "holm"), digits = 5)
round(ps, digits=5)

# our conclusions don't change


# 7. graph
d$G = as.factor(d$G)
mY1 = lm(Y.1 ~ G, data=dW)
mY2 = lm(Y.2 ~ G, data=dW)
mY3 = lm(Y.3 ~ G, data=dW)


pX1 = expand.grid(G = c(1,2), X = 1)
pX2 = expand.grid(G = c(1,2), X = 2)
pX3 = expand.grid(G = c(1,2), X = 3)

pY1 = modelPredictions(mY1, pX1)
pY2 = modelPredictions(mY2, pX2)
pY3 = modelPredictions(mY3, pX3)

pYs = rbind(pY1, pY2, pY3)

library(ggplot2)
library(cowplot)
bp1 = ggplot(data=pYs, aes(x=as.factor(X), y = Predicted, fill=as.factor(G))) +
  geom_bar(mapping = aes(), stat='identity', width=.5, data=pYs, position_dodge()) +
  geom_errorbar(width=.25, aes(ymin = CILo, ymax = CIHi), stat='identity', data=pYs, position_dodge(.5)) +
  coord_cartesian(ylim = c(3.8, 80), expand=T)
bp1