1.1 Getting the data

Open the W3 Test on Blackboard and click on the link in Question Q1 to download your ‘personal’ version of the data. This will download your copy of the CSV file blooms.csv with the data for this experiment. (In my analysis, I’ve read my CSV into a data frame called Blooms.)

1.2 Setting up for the analysis

Run summary (or skim or some other summary tool) on the data. You’ll note that the categorical EVs are coded by numbers:

• bed: 1, 2, 3;
• water: 1, 2, 3 for once, twice, three times a week;
• shade: 1, 2, 3, 4 for no shade, 1/4, 1/2 and fully shaded;

So you’ll need to explicitly declare them as factors. As usual, we’ll be using sum contrasts rather than the default treatment contrasts — so don’t forget to set up the contrasts as well.

1.3 Fitting the model

TASK: Let’s start by fitting a model that is appropriate for this experimental design. In a blocked design, the blocking factor (here, bed) should always come first in the model formula (can you think why?²), so:

m.1 <- lm(sqblooms ~ bed + water + shade, data = Blooms)

As you can see, I have called my model m.1 and I suggest you do the same, because we will be referring back to this model later.

2.1 Using just the anova function in base R.

Recall that the anova function of R gives sequential SSQ. For this assessment of orthogonality, you need to fit another version of the model with the same full set of EVs, but change the order of the EVs in the model equation. You could do this like so – the first line is not needed as we have already fitted m.1, this is only to help you see the difference:

m.1  <- lm(sqblooms ~ bed + water + shade, data = Blooms)  # not needed!
m.1a <- lm(sqblooms ~ water + shade + bed, data = Blooms)

# or, simpler and more elegantly, you could `update` the model `m.1`:
m.1a <- update(m.1, . ~ water + shade + bed)

Here, I have switched the order to water + shade + bed – you could switch it to something else, since with three EVs, there are six ways of arranging them.³

TASK: Print the ANOVA tables from the two versions of the model to obtain the sequential SSQs for both. Are they different? Write down in your notes what this means for the orthogonality of the data.

BB TEST Q1 [5 marks]: What is the value of the sequential SSQ for bed in your your model m.1?

2.2 Using Anova from the car package

The Anova function from the car package makes it easy to get adjusted SSQ for all EVs in one single ANOVA table — this is also known as a “Type II” ANOVA. However, Anova provides only the adjusted SSQ, not the sequential SSQ.

For this orthogonality check, you need to compare the sequential and the adjusted SSQ from the same model.

BB TEST Q2 [5 marks]: What is the value of the adjusted SSQ for bed in your your m.1?

Compare this value to that of the sequential SSQ from Q1 – are the two values the same or different? What about the SSQs for water and shade?

BB TEST Q3 [10 marks]: What is your overall assessment of the orthogonality of the Blooms dataset?

3.1 A first answer…

You can get some kind of answer from the \(P\)-value for bed in the ANOVA table for m.1. Instead of applying a simplistic ‘significance cut-off’ of \(\alpha = 0.05\), I want you to think of the \(P\)-value as a measure of the “strength of the evidence” against the Null hypothesis, like this:

• \(P \geq 0.1:\) little or no evidence.
• \(0.1 > P \geq 0.01:\) some evidence.
• \(0.01 > P \geq 0.001:\) strong evidence.
• \(P < 0.001:\) compelling evidence.

In our case, the Null hypothesis is that bed explains none of the observed variation in flower-head counts; if true, the blocked design is unnecessary and it would be safe to do the experiments without blocking in the future. Thus, the \(P\)-value measures the evidence against this idea that blocking is unnecessary – the smaller \(P\), the stronger the evidence against.

But here is the catch: the \(P\)-value cannot give the farmer the answer he would like. A small \(P\)-value is the answer the farmer does not want to hear (“It would be ill-advised to not block your experiments in the future”); but a large \(P-\)-value doesn’t absolve him from blocking in the future. This is because a large \(P\) is not evidence for the Null, only lack of evidence against the Null; you can never use \(P\) to prove the Null. For our farmer, a large \(P\) only means that he is none the wiser: it leaves him with no good evidence that the beds are different, but this doesn’t prove that they are the same, and he should still block his experimental designs.

BB TEST Q4 [10 marks]: Given the \(P\)-value for bed in your data, what advice should you give the flower grower?

3.2 A deeper look…

To better understand whether blocking for bed was worthwhile, the farmer re-analysed his the data without using bed as a blocking factor. The logic here is that we can pretend in the analysis that blocking did not happen, to see how much information we would have lost if we hadn’t used a blocked design (and thus how much we gained by doing a blocked design).

Do what the farmer did, and compare the results from two models:

• the full model (m.1), which matches the actual experimental design and therefore analyses the experiment as a randomised block design with bed as the blocking factor; and
• a model that analyses the same experiment as a fully randomised design by omitting the blocking factor from the analysis.

TASK: In your notes, write down how dropping bed has changed the model results.

• How have the residual SSQ changed, and why?
• How have the residual DoF changed, and why?
• What about the adjusted SSQ and the \(F\)-ratios for water and shade?
• What happened to the \(P\)-values for water and shade?

In particular, I want you to think about the changes in the standard errors (SEs) for the model coefficients – what do they tell you? The best way to think about how dropping bed from the model affects the SEs is in terms of “fold-changes”. For example: if the SE for the effect of water is 0.11017 with bed in the model, and 0.13777 after dropping bed, then the fold-change in SE is simply \[SE_\mathrm{after}/SE_\mathrm{before} = 1.25\] so you can say that the SE is \(1.25\times\) larger when bed is dropped from the model.

BB TEST Q5 [10 marks]: What is the fold-change in the standard error of the effect estimates for water due to removing the blocking factor from the model? Use the round function to round your answer to two decimals.

NOTE: You don’t have to write code for extracting the SE values from the model fits (although you could). Your R console output for SE will be accurate to five decimals, so you can simply copy/paste the SE values from the table back into the console, like this (these values are from my version of the data):

round(0.13777/0.11017, digits = 2)

In your notes, write down what the blocked design does for the precision of the estimated effects.

5.1 Know your coefficients

We get the coefficients from the model fit like this:

cf <- coef(m.1)
cf  # print to see what is what

## (Intercept)        bed1        bed2      water1      water2      shade1 
##  4.02902778  0.06197222  0.38047222 -0.41102778  0.37313889  0.09652778 
##      shade2      shade3 
##  0.29341667 -0.11913889

Now we need to remember what the model coefficients actually mean, and how we get from the coefficients to the means that our model has estimated for the different groups.

For plotting how watering levels affect the expected bloom count, we want the estimates to be kind of “independent of the amount of shade” and “regardless of the bed” — but of course shade does affect bloom count, and the beds differ too, so what to do? The most useful prediction of the bloom counts are for an average amount of watering, and in a ‘typical’ average bed.

NOTE: Now it pays off that we have used sum contrasts for all EVs: Remember that with sum contrasts, all estimated effects are relative to the grand mean \(\mu\).

Therefore, instead of working out the answer to “What is the expected bloom count in bed 2 if watered twice a week in full shade?” (which would require coefficients for all EVs), here we can simply drop the other EVs from the calculation:

• mean bloom in ‘watered once’ group = \(\mu + \alpha_1\), i.e., simply not including the coefficients for bed and shade; likewise,
• mean bloom in ‘watered twice’ group = \(\mu + \alpha_2\), and
• mean bloom in ‘watered thrice’ group = \(\mu - \alpha_1 - \alpha_2\)

Can you see why this would be trickier if we had used treatment contrasts?

5.2 The R craft side of things

To get the values of the individual coefficients out of cf you have two options. Either use numerical indexing, or use the associated names. E.g., to get the value of the mean in the ‘watered once’ group:

cf[1] + cf[4]

## (Intercept) 
##       3.618

cf["(Intercept)"] + cf["water1"]

## (Intercept) 
##       3.618

Numerical indexing makes for compact but unreadable code, and it is easy to mis-count the index and thus mis-calculate. I’d use the names.

For subsequent plotting by ggplot, it is useful to collect and organise your calculations into a new data frame, like this:

water.pred <- data.frame(
  water = factor(c("once", "twice", "thrice"),
                 # the next bit is just to ensure order along x-axis when plotting
                 levels = c("once", "twice", "thrice")),
  sqbloom = c(
    cf["(Intercept)"] + cf["water1"],             
    cf["(Intercept)"] + cf["water2"],
    cf["(Intercept)"] - cf["water1"] - cf["water2"]
  )
)

Normally, I would not be seen dead using bar plots, but we’ll make an exception here…

library(ggplot2)

# I'm saving the beginnings of the plot in `p` and add the `geometry` later,
# so that we can reuse the plot using different plot styles:
p <- ggplot(water.pred, aes(x=water, y=sqbloom)) + 
  labs(x="Weekly watering regime", y="bloom count (square root)")

p + geom_col()  # this is ggplot2 function for `bar charts`

p + geom_point(size=4) + geom_line(group=1)

5.3 Over to you!

TASKS:

1. Replicate my steps and create the plot for water.
2. Create the corresponding plot for shade.