##### Code for Interactions in R


### Importing data


datum=read.csv(file.choose()) ### creates a data frame named datum from csv file

head(datum) ### column headers of data file plus first 6 rows of data 


###Plot Data

plot(Mass~Protein,data=datum) ### scatterplot of data used in example


### Run Standard ANCOVA - treat Protein as continuous (linear) variable 

results=lm(Mass~Protein+Supplement,data=datum) ### linear model of example data 
summary(results) ### prints a summary of 'results'


### Run ANCOVA with interactions

results2=lm(Mass~Protein+Supplement+Protein:Supplement,data=datum) #model with interactions
summary(results2) 
### Note that interaction term is significant
### Meaning of interaction term is complicated - difference in slope
### Easiest way to deal with interactions is demonstrate significant
### Then pull apart the data

results2b=lm(Mass~Protein,data=datum,subset=c(Supplement=="Yes"))
summary(results2b)
results2c=lm(Mass~Protein,data=datum,subset=c(Supplement=="No"))
summary(results2c)
###Focus on Effect of Protein under the two supplemental levels



### Run Multi-factor ANOVA with interactions

datum$ProCat=as.factor(datum$ProCat)
results3=lm(Mass~ProtCat+Supplement+ProCat:Supplement,data=datum)
summary(results3) ### Note the multiple interaction terms
anova(results3) ### Test whether interaction is significant as a whole
### Note that because we only want interaction term; don't care that type I SOS

### Interaction term is significant
### Pull apart - Can look at effect of Protein under 2 supplement levels:

results3b=lm(Mass~ProCat,data=datum,subset=c(Supplement=="Yes"))
summary(results3b)
results3c=lm(Mass~ProCat,data=datum,subset=c(Supplement=="No"))
summary(results3c)

### Or Can look at effect of supplement under different Protein levels:

results4b=lm(Mass~Supplement,data=datum,subset=c(ProCat=20)
results4c=lm(Mass~Supplement,data=datum,subset=c(ProCat=25)


### You can also take advantage of the emmeans package
### Won't give you differences. Instead will calculate response
### at different levels (i.e., for graphing)

library(emmeans)
emm=emmeans(results3, specs=~ProCat+Supplement+ProCat:Supplement)
plot(emm)

### Alternative
emm=emmeans(results3,specs="ProCat",by="Supplement")
plot(emm)

### Second Alternative
emmip(results3,Supplement~ProCat)