#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)