#Dataset for the analysis
D_M <- c(16, 13, 12, 15, 11, 12, 14, 13)
D_S <- c(16, 10, 11, 12, 6, 8, 14, 12)
Im_M <- c(14, 16, 17, 15, 13, 17, 13, 16)
Im_S <- c(13, 7, 3, 10, 4, 2, 4, 9)
I_M <- c(15, 15, 12, 14, 13, 11, 11, 12)
I_S <- c(15, 10, 11, 7, 5, 12, 6, 8)
D <- rep(c("Densitization","Implosion","Insight"), c(16, 16, 16))
S <- rep(rep(c("Mild","Severe"),c(8,8)),3)
#1a. What are the factors in this study?
#The factors are therapy and phobia.
#1b.How many levels are associated with each factor?
##Five levels in total. phobia therapies with 3 levels namely;Desensitization, Implosion and Insight) and 2 levels for phobia (mild or severe).
##1c.How many treatments are in this study?
# Six treatments in total.
## 2. What type of treatment structure is present in this study?
#This is a factorial treatment structure, because it contains a unique combination of factor levels
#3.Construct an interaction plot of the means associated with unique factor level combinations. What do you conclude from this display.
myDat <- data.frame(fear=c(D_M,D_S,Im_M,Im_S,I_M,I_S),phobia=S, therapy=D)
with(myDat,interaction.plot(x.factor=D,trace.factor=S,response =fear))

#The plot display shows that there is an interaction between the two factors.
# 4. Write the model for this analysis. Describe all components of this model.
# (fear)ij = µ + therapy +phobia + (therapy*phobia)ij + ??ijk
# i = 1, ., a (Factor A levels- Therap)
# j = 1, ., b (Factor B levels- Phobia)
# k = 1, ., nij
# where µ = overall mean (with constraints ?? ??i = 0, ?? ??j = 0, ??i (????)ij = 0, ??j (????)ij = 0)
# ??i = main effect of Factor A (level of Severity)
# ??j = main effect of Factor B (therapy)
# (????)ij = interaction of level of severity and therapy
# ??ij = random error ~ independent N(0, ??2)
#5. Formally test if there is an interaction between therapy type and severity of phobia.
phobia_model<-lm(fear~D*S,data=myDat)
anova(phobia_model)
## Analysis of Variance Table
##
## Response: fear
## Df Sum Sq Mean Sq F value Pr(>F)
## D 2 17.167 8.583 1.1495 0.326572
## S 1 275.521 275.521 36.8972 3.105e-07 ***
## D:S 2 92.667 46.333 6.2049 0.004355 **
## Residuals 42 313.625 7.467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#(Fear)ij = µ + alphai + Betaj + (alpha*Beta)ij + ??ijk
# Ho: aplpha * Beta = 0
# Ha: aplpha * Beta =/ 0
# Fvalue:6.205
# Df:2, 42
# P value= 0.00436 (less than 0.05 level of signifaicance)
# Conclusion:
#6. Should you test main effects in this analysis. Explain why or why not.
#Since there is an interaction between phobia and therapy,there won't be a need to test for the main effects of the factors. We can however do a pairwise-t-test for phobia as it shows significant effect on fear (p- value<0.05)
pairwise.t.test(myDat$fear, myDat$phobia, p.adj = "none")
##
## Pairwise comparisons using t tests with pooled SD
##
## data: myDat$fear and myDat$phobia
##
## Mild
## Severe 1.8e-06
##
## P value adjustment method: none
#Also, since our phobia variable contains only two levels, there is no need to conduct follow-up comparisons,
#7.Are the model assumptions (equal variance, normality, independence, absence of outliers) met for these data?
plot(phobia_model)


## hat values (leverages) are all = 0.125
## and there are no factor predictors; no plot no. 5

hist(myDat$fear) #Check for outliers

#Checking the distribution of residuals.
residual<-resid(phobia_model)
myDat <- data.frame(fear=c(D_M,D_S,Im_M,Im_S,I_M,I_S), level=S, therapy=D,residual=residual)
hist(myDat$residual)#This again checks for outliers

#8. Reformulate this analysis as a one-way anova model. Fit this model and use Tukey's HSD procedure to determine which means differ. What is your overall conclusion from this.
#First,obtain the response from the replicated treatments
D_M <- c(16, 13, 12, 15, 11, 12, 14, 13)
D_S <- c(16, 10, 11, 12, 6, 8, 14, 12)
Im_M <- c(14, 16, 17, 15, 13, 17, 13, 16)
Im_S <- c(13, 7, 3, 10, 4, 2, 4, 9)
I_M <- c(15, 15, 12, 14, 13, 11, 11, 12)
I_S <- c(15, 10, 11, 7, 5, 12, 6, 8)
#Assign an object to the treatment
Treatment<-rep(c("D_M","D_S","Im_M","Im_S","I_M","I_S"),c(8,8,8,8,8,8))
#Combine the treatments and measured response in dataframe.
myDat2 <- data.frame(fear2=c(D_M,D_S,Im_M,Im_S,I_M,I_S),Treatment=Treatment)
myDat2
## fear2 Treatment
## 1 16 D_M
## 2 13 D_M
## 3 12 D_M
## 4 15 D_M
## 5 11 D_M
## 6 12 D_M
## 7 14 D_M
## 8 13 D_M
## 9 16 D_S
## 10 10 D_S
## 11 11 D_S
## 12 12 D_S
## 13 6 D_S
## 14 8 D_S
## 15 14 D_S
## 16 12 D_S
## 17 14 Im_M
## 18 16 Im_M
## 19 17 Im_M
## 20 15 Im_M
## 21 13 Im_M
## 22 17 Im_M
## 23 13 Im_M
## 24 16 Im_M
## 25 13 Im_S
## 26 7 Im_S
## 27 3 Im_S
## 28 10 Im_S
## 29 4 Im_S
## 30 2 Im_S
## 31 4 Im_S
## 32 9 Im_S
## 33 15 I_M
## 34 15 I_M
## 35 12 I_M
## 36 14 I_M
## 37 13 I_M
## 38 11 I_M
## 39 11 I_M
## 40 12 I_M
## 41 15 I_S
## 42 10 I_S
## 43 11 I_S
## 44 7 I_S
## 45 5 I_S
## 46 12 I_S
## 47 6 I_S
## 48 8 I_S
#One-way analysis of variance with Treatment as a single factor
phobia_one_way<-aov(fear2~Treatment,data=myDat2)
phobia_one_way
## Call:
## aov(formula = fear2 ~ Treatment, data = myDat2)
##
## Terms:
## Treatment Residuals
## Sum of Squares 385.3542 313.6250
## Deg. of Freedom 5 42
##
## Residual standard error: 2.732629
## Estimated effects may be unbalanced
#Post-hoc test to compare the means using TukeyHSD
fear_posthoc <- TukeyHSD(phobia_one_way,conf.level=0.95)
fear_posthoc
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = fear2 ~ Treatment, data = myDat2)
##
## $Treatment
## diff lwr upr p adj
## D_S-D_M -2.125 -6.20378894 1.95378894 0.6314778
## I_M-D_M -0.375 -4.45378894 3.70378894 0.9997730
## I_S-D_M -4.000 -8.07878894 0.07878894 0.0573900
## Im_M-D_M 1.875 -2.20378894 5.95378894 0.7429090
## Im_S-D_M -6.750 -10.82878894 -2.67121106 0.0001785
## I_M-D_S 1.750 -2.32878894 5.82878894 0.7936017
## I_S-D_S -1.875 -5.95378894 2.20378894 0.7429090
## Im_M-D_S 4.000 -0.07878894 8.07878894 0.0573900
## Im_S-D_S -4.625 -8.70378894 -0.54621106 0.0181182
## I_S-I_M -3.625 -7.70378894 0.45378894 0.1069124
## Im_M-I_M 2.250 -1.82878894 6.32878894 0.5733210
## Im_S-I_M -6.375 -10.45378894 -2.29621106 0.0004249
## Im_M-I_S 5.875 1.79621106 9.95378894 0.0013151
## Im_S-I_S -2.750 -6.82878894 1.32878894 0.3525751
## Im_S-Im_M -8.625 -12.70378894 -4.54621106 0.0000020
plot(fear_posthoc)
