##### Helicopter ANOVAs II

# Finally (!) we analyze copter data the correct way. No more dinking around 
# with t-tests and wrong ANOVAs – today we apply the full
# analysis the way we designed the experiment. We'll do that in a couple of
# steps, so you can see how those steps matter.

# Reminder: Three treatments potentially affect helicopter flight times: 
# wing length, body width, and whether wings were folded or not. This was a
# factorial design, where treatments were combined in all possible combinations.
# Replicate copters were organized by Groups (or blocks), and the effect of
# stairs is a covariate.

### Step 1. Import & attach copter data - see weeks 3 or 4 on our class web
# page
# Below I assume you called the imported file "data”.

### Step 2. Make factor variables (e.g., fwl, fbw & ffold – those names are
# assumed below) for WL, BW & Fold. Why? Because we are doing ANOVAs with
# categorical treatments.

### NOTE: Here we defer worries about homogeneity of variance and normality 
# below. We test them, but do so on a final, most-efficient model.

### Step 3. Squint at potential interactions again:

interaction.plot(fwl, trace.factor=fbw, response=Time, fun=mean, type="b")
interaction.plot(fwl, trace.factor=ffold, response=Time, fun=mean, type="b")
interaction.plot(ffold, trace.factor=Group, response=Time, fun=mean, type="b")

### What do you see? What should we then expect in ANOVA results?

### An important note on the “block” effect: Each Group had 10 replicates per
# treatment combination. A simple randomized complete block experiment would
# have only one replicate of each treatment combination in a Group, but then
# could not evaluate a Group x Treatments interaction. So Group here was more
# than a simple block – we can evaluate how Group interacts with a treatment.
# And each Group was free to cut, fold, and drop copters in their own style –
# thus interactions may matter.

### Step 4. Analyze the data according to experimental design:

copters1 <- lm(Time ~ Step + fwl*fbw*ffold + Group + Step)
summary(copters1)
summary.aov(copters1)

# Notice the difference from last week? Including Step here really helped.
# Notice the code shortcut above? You ran a linear regression model (lm) and
# then asked for the ANOVA table that comes from it. We get everything in 3
# lines! And remember that using * above actually means [fwl + fbw + ffold
# + fwl X fbw + fwl X ffold + fbw X ffold + fwl X fbw X ffold].

### What can you conclude? 

### Would you report this complete model in a paper?
# There are (at least) two philosophies about this question:
# A) Standard: Report the most complete model as-is (e.g., copters1 above). That
# is the study 	design, even if multiple terms are not significant. Truth in
# advertising, so to speak. 
# B) Crawley (2013; The R Book) advocates iterative simplification. This means
# eliminating nonsignificant terms to report the most parsimonious model that
# represents significant effects. 

# What should you do in your own research? Explain what you did in analyses,
# the same way you would explain a method to get the data. If you simplify,
# don't just present the last, "winning" model as if that was everything. 

### So let's simplify. We use stepwise AIC - which uses Akaike Information 
# Criterion to ID the most parsimonious model.

library(MASS) # this should already be installed, but if not, install it first.
stepAIC(copters1) # iteratively removes terms and compares by AIC 

# What is shown? It computed AIC for the starting model (copters1)
# It then iteratively tried alterative models with one predictor and listed 
# outcomes in a table, where The Lowest AIC Is Most Plausible. It shows that 
# model (lacking the 3-way interaction term) is most plausible among other 
# choices. Then it shows all the final model coefficients.

### You can make your own table to look like copter1 by using the formula shown
# in output, like this:
copters2 <- lm(Time ~ Step + fwl + fbw + ffold + Group + fwl:fbw + 
    fwl:ffold + fbw:ffold)
summary(copters2)
summary.aov(copters2)

### Now we check the final, most-parsimonious model for assumptions.
# Install the package called performance. This takes a while but is worth it.
install.packages('performance')
library(performance)
check_model(copters2)

### What do you see? Does your final, most-parsimonious model meet assumptions?
# Notice it is the RESIDUALS (i.e., left-over scatter) we worry about here.

### Summary:
# you efficiently combined linear regression and ANOVA models
# you used AICs to obtain the most parsimonious, complete model
# you efficiently evaluated model assumptions
# I call that a good day