library(tidyverse)
library(marginaleffects)STAD29 Assignment 5
You are expected to complete this assignment on your own: that is, you may discuss general ideas with others, but the writeup of the work must be entirely your own. If your assignment is unreasonably similar to that of another student, you can expect to be asked to explain yourself.
If you run into problems on this assignment, it is up to you to figure out what to do. The only exception is if it is impossible for you to complete this assignment, for example a data file cannot be read. (There is a difference between you not knowing how to do something, which you have to figure out, and something being impossible, which you are allowed to contact me about.)
You must hand in a rendered document that shows your code, the output that the code produces, and your answers to the questions. This should be a file with .html on the end of its name. There is no credit for handing in your unrendered document (ending in .qmd), because the grader cannot then see whether the code in it runs properly. After you have handed in your file, you should be able to see (in Attempts) what file you handed in, and you should make a habit of checking that you did indeed hand in what you intended to, and that it displays as you expect.
Hint: render your document frequently, and solve any problems as they come up, rather than trying to do so at the end (when you may be close to the due date). If your document will not successfully render, it is because of an error in your code that you will have to find and fix. The error message will tell you where the problem is, but it is up to you to sort out what the problem is.
Packages
Teak
Teak is a tropical hardwood tree that grows in south and southeast Asia. The wood has been used for shipbuilding because of its durability and resistance to weather. Teak trees are grown from seeds; the seeds are planted in small pots in a greenhouse, and when they have germinated (become small plants), the small plants are planted in the ground outside. In this experiment, the small plants are planted in the ground using one of two methods (in column Plant):
- using a
crowbarto make a hole for the small plant to go into - digging a larger
pitand putting the small plant at the bottom.
In addition, each small plant was planted outside when its roots reached a randomly-chosen length: 4, 6, or 8 inches, denoted R4, R6, or R8 in column Root.
After one year, the Height of each plant was measured, in inches. (Teak trees grow very slowly.) The experimenters were interested in whether the planting method or the root length, or the combination of both, had any effect on Height. The data are in http://datafiles.ritsokiguess.site/teak.csv. There is an additional column Block that you can ignore.
- (1 point) Read in and display (some of) the data.
As is the way:
my_url <- "http://datafiles.ritsokiguess.site/teak.csv"
teak <- read_csv(my_url)Rows: 24 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): Block, Plant, Root
dbl (1): Height
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.teakExtra: There are three root lengths and two planting methods, so there are \(2 \times 3 = 6\) treatment combinations, and there are 24 observations, so there are (apparently) \(24 / 6 = 4\) replicates at each treatment combination:
teak %>% count(Root, Plant)The small plants were planted in four different fields (this is the column Block). You might imagine that growing conditions could vary from one field to another, and so we ought to incorporate Block into the analysis. I’m not asking you to do this, but I will explore later.
- (2 points) Make a suitable graph of the three variables of interest.
These are Plant and Root (categorical) and Height (quantitative), so a grouped boxplot:
ggplot(teak, aes(x = Root, y = Height, fill = Plant)) + geom_boxplot()
Both categorical variables are explanatory, but there are three levels of Root and only two of Plant, so I would put Root on the \(x\)-axis and have only two colours for Plant.
Decide for yourself if you like it better or worse the other way around:
ggplot(teak, aes(x = Plant, y = Height, fill = Root)) + geom_boxplot()
I think having Plant as x wastes a lot of space on the \(x\)-axis, myself, and you have to mentally assess groups of three boxplots rather than groups of two.
Extra: you might be worried about some of the boxplots being much taller than others, which implies unequal spreads, but bear in mind that each of those boxplots is based on only four observations, so there is plenty of room for craziness even if normality and equal spreads is actually fine.
- (2 points) What does your plot of the previous question tell you about a likely interaction between planting method and root length?
Using my first boxplot: when the planting method is Pit, the plant has a greater height on average than when the planting method is Crowbar, and the amount by which it is higher is about the same for each root length. Thus, I would expect to see a Plant effect that is the same for each root length: that is, no interaction between planting method and root length.
Extra: if you want, you can draw an interaction plot and see whether it tells the same story as the boxplots (it should, except for the fact that an interaction plot uses means and the boxplots use medians):
teak %>%
group_by(Root, Plant) %>%
summarize(mean_height = mean(Height)) %>%
ggplot(aes(x = Root, y = mean_height, colour = Plant, group = Plant)) +
geom_line()`summarise()` has grouped output by 'Root'. You can override using the
`.groups` argument.
and those two traces do indeed look pretty close to parallel.
The boxplot version of an interaction plot would be to draw lines through the medians of each boxplot: blue lines through the blue boxplots, and red lines through the red ones.
- (2 points) Run a suitable analysis of variance of the variables of interest, and display the results.
Include an interaction first off, even though you probably think it won’t be significant:
teak.1 <- aov(Height ~ Root * Plant, data = teak)
summary(teak.1) Df Sum Sq Mean Sq F value Pr(>F)
Root 2 203.15 101.57 32.131 1.15e-06 ***
Plant 1 187.60 187.60 59.344 4.18e-07 ***
Root:Plant 2 3.78 1.89 0.598 0.56
Residuals 18 56.90 3.16
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- (2 points) What do you conclude from your analysis of variance?
Test the interaction first. With a P-value of 0.56, this is nowhere near significant. This is where you stop.
You might be thinking that the main effects will be (strongly) significant, but it’s too soon to conclude that now.
- (4 points) Complete your analysis (or explain why it is already complete). What do you conclude, in the context of the data?
Take out the non-significant interaction and try again:
teak.2 <- aov(Height ~ Root + Plant, data = teak)
summary(teak.2) Df Sum Sq Mean Sq F value Pr(>F)
Root 2 203.15 101.57 33.48 4.14e-07 ***
Plant 1 187.60 187.60 61.83 1.52e-07 ***
Residuals 20 60.68 3.03
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now you can conclude that there are significant main effects, of root length (P-value \(4.1 \times 10^{-7}\)) and of planting method (P-value \(1.5 \times 10^{-7}\)).
To find out what kind of effects they are, run Tukey (which you need to do because there are three different root lengths):
TukeyHSD(teak.2) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Height ~ Root + Plant, data = teak)
$Root
diff lwr upr p adj
R6-R4 5.6875 3.48403 7.89097 0.0000067
R8-R4 6.5625 4.35903 8.76597 0.0000008
R8-R6 0.8750 -1.32847 3.07847 0.5825209
$Plant
diff lwr upr p adj
pit-crowbar 5.591667 4.108292 7.075041 2e-07The Plant analysis here doesn’t tell us anything we haven’t already seen: when the Pit planting method is used, the height is on average greater than for the Crowbar method. The Root analysis, however, does tell us something new: the height for plants that had root length 4 inches is significantly less than for plants with either of the other two root lengths (which are not significantly different).
There is no value in doing simple effects here, because the interaction is not significant. The story is the same as on the worksheet problem, where the interaction was also not significant, I had you do some simple effects anyway, and you found there was nothing to see.
Extra: as on the worksheet problem, we have a blocking variable, which is actually called Block here. The conclusions from our analysis above were pretty clear, but it might be instructive to add it to the analysis and see what happens. Remember we add a block main effect (only), so this. I am using update because I am making a small change to a model I already fitted:
teak.3 <- update(teak.1, . ~ . + Block)
summary(teak.3) Df Sum Sq Mean Sq F value Pr(>F)
Root 2 203.15 101.57 75.696 1.45e-08 ***
Plant 1 187.60 187.60 139.806 5.29e-09 ***
Block 3 36.77 12.26 9.135 0.00111 **
Root:Plant 2 3.78 1.89 1.409 0.27499
Residuals 15 20.13 1.34
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The P-value for the interaction is smaller than it was, but still nowhere near significant, so we yank the interaction again. (Note that the block effect is significant, which is not really of interest in itself, but says that it was worth including the blocks in the analysis because they do differ somehow.)
Thus, next:
teak.4 <- update(teak.3, . ~ . - Root:Plant)
summary(teak.4) Df Sum Sq Mean Sq F value Pr(>F)
Root 2 203.15 101.57 72.222 4.9e-09 ***
Plant 1 187.60 187.60 133.391 1.8e-09 ***
Block 3 36.77 12.26 8.716 0.00101 **
Residuals 17 23.91 1.41
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The two main effects are still strongly significant, with P-values even smaller than they were before. We can run Tukey again:
TukeyHSD(teak.4) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Height ~ Root + Plant + Block, data = teak)
$Root
diff lwr upr p adj
R6-R4 5.6875 4.1663509 7.208649 0.0000001
R8-R4 6.5625 5.0413509 8.083649 0.0000000
R8-R6 0.8750 -0.6461491 2.396149 0.3268084
$Plant
diff lwr upr p adj
pit-crowbar 5.591667 4.570203 6.61313 0
$Block
diff lwr upr p adj
II-I -0.4833333 -2.42960184 1.46293517 0.8933553
III-I 1.5500000 -0.39626851 3.49626851 0.1461745
IV-I -1.9166667 -3.86293517 0.02960184 0.0543915
III-II 2.0333333 0.08706483 3.97960184 0.0389207
IV-II -1.4333333 -3.37960184 0.51293517 0.1948952
IV-III -3.4666667 -5.41293517 -1.52039816 0.0005065The conclusions about Plant and Root are as before, albeit with smaller P-values. We now get a comparison of the blocks, in case that is of interest; about all that is worth saying is that some of them differ significantly as well.
This wasn’t really a very good sales job about the value of including blocks in the analysis. The time when it makes a difference is when the significance is less clear; then it can be the difference between a non-significant effect (without blocks) and a significant one (with). All you saw here was that the P-values for the main effects became even smaller than they already were.
Rhizobium
Rhizobium is a bacteria, growing on the roots of clover and alfalfa, that fixes nitrogen from the atmosphere into a chemical form the plants can use. There are two forms of rhizobium, one that grows on the roots of clover and one that grows on the roots of alfalfa.
A study was made on clover plants. The treatments were combinations of bacterial cultures in which the plants were grown. There are two factors:
strain: one of six clover rhizobium culturescomb: ifclover, the strains were as above; ifclover+alfalfa, each strain was mixed with alfalfa rhizobium cultures.
The response variable was the milligrams of nitrogen per gram of dry plant weight, Npg in the dataset. The data are in http://datafiles.ritsokiguess.site/rhiz-clover.csv.
The research question was whether the strains mixed with alfalfa rhizobium cultures had lower nitrogen content than the strains with pure clover rhizobium cultures. If so, this would indicate that the alfalfa and clover rhizobium cultures were competing against each other, a phenomenon known as “antibiosis”. This effect might differ across strains.
Five observations were taken at each of the \(6 \times 2 = 12\) combinations of strain and comb.
- (1 point) Read in and display (some of) the data.
As is the way:
my_url <- "http://datafiles.ritsokiguess.site/rhiz-clover.csv"
clover <- read_csv(my_url)Rows: 60 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): comb, strain
dbl (3): weight, nitro, Npg
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.cloverThere are correctly \(12 \times 5 = 60\) observations, and the variables of interest are as described. (The value of Npg is weight divided by nitro, which is why it has a lot of decimals.)
Extra: these are the rhiz.clover data from the HH package. This package is associated with a book “Statistical Analysis and Data Display” by Heiberger and Holland (2015 edition). A detailed description of the experiment starts on page 400; there were actually two experiments, one with clover plants and one with alfalfa plants, that were conducted in parallel (the book, as befits its title, includes some of the analysis we will do and a number of graphs).
I struggled to come up with a coherent description of what is going on (I am not a biologist).
- (3 points) Create an interaction plot, briefly justifying your choice of which variable goes where on your plot.
This is a plot of the mean value of Npg (\(y\)-axis) for each combination of comb and strain. One of the explanatory variables goes on the \(x\)-axis, and the other is colour.
The first thing is to get the means:
clover %>%
group_by(comb, strain) %>%
summarize(mean_npg = mean(Npg)) -> clover_means`summarise()` has grouped output by 'comb'. You can override using the
`.groups` argument.clover_meansAs a check, there should be twelve means (there are twelve treatment combinations). Each of those means is based on five observations.
For the plot, I think it is best to put the explanatory variable with more levels (strain) on the \(x\)-axis, and have the other one, comb, as colours (so there are not too many colours to distinguish; here there are only two). In the code, you need an x, y, and colour, as you would expect, but remember that you also need a group (to make sure the right points get joined by lines) that should be the same variable as the colour:
ggplot(clover_means, aes(x = strain, y = mean_npg,
colour = comb, group = comb)) + geom_line()
You are also free to do the summaries and the plot in one pipeline (without naming the dataframe of means), or even to use stat_summary from the lecture notes if you can figure out how it works. (I find it easier to actually compute the means myself, and then I know what is going on.)
- (2 points) What is the main message from your interaction plot? Explain briefly.
The major reason for drawing an interaction plot is to see whether you are likely to have an interaction, and you assess that by making a call about whether the lines on your plot are parallel. Here, it seems clear that they are not: the strain 3DOk5 seems very different from the others. So we would expect to see an interaction.
Extra 1: that one very different strain is showing the difference in comb that is consistent with antibiosis (red higher than blue), but for the other strains there is not much difference, or the blue is higher than the red.
Extra 2: another graph I might have asked you to draw is a grouped boxplot. Here, that has a rather disturbing look:
ggplot(clover, aes(x = strain, y = Npg, fill = comb)) + geom_boxplot()
The Npg for strain 3DOk5 and comb clover has nitrogen values that are very much more spread out than for the other treatment groups, and there also appear to be some outliers. If you see this graph, you might be asking for a version of Mood’s median test that applies to two-way ANOVA, except that we don’t have one. Bear in mind, though, that each boxplot is only based on five observations, so there can be outliers more or less by chance (three or four observations close together and one or two that happen to be different). It is perhaps better to plot the individual observations rather than the boxes, which can be done with very few changes to the code:
ggplot(clover, aes(x = strain, y = Npg, colour = comb)) + geom_point()
I replaced the geom_boxplot with geom_point, and the fill with colour, because points don’t have insides that would be coloured with fill.
Two of the 3DOk5 - clover observations appear to be high outliers, and probably also one of the 3DOk1 - clover observations, but there don’t appear to be any other outliers, despite what the boxplots say. Also, the big box for 3DOk5 - clover was caused by the second-highest observation being large, so it looked as if the IQR was large (so that the box is tall but there are no outliers).1
If you discount these outliers, normality and equal spread looks pretty good, but of course you should not just remove them. Possibilities to explore (with a biologist) are that there is something about the biological mechanism that occasionally produces very large nitrogen values, or perhaps there is just something else different about these plants.
As you see, practical data analysis can have a lot of issues beyond just running an ANOVA. But the point here was to get you running an ANOVA correctly, so I didn’t ask you to explore these issues.
- (2 points) Fit a suitable two-way analysis of variance, and display the results.
A suitable two-way ANOVA includes an interaction, one that we are expecting to be significant:
clover.1 <- aov(Npg ~ comb * strain, data = clover)
summary(clover.1) Df Sum Sq Mean Sq F value Pr(>F)
comb 1 6.9 6.88 0.531 0.46955
strain 5 642.3 128.45 9.916 1.47e-06 ***
comb:strain 5 228.2 45.65 3.524 0.00857 **
Residuals 48 621.8 12.95
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1You might care to note that, as you expected, the interaction is significant.
Extra: this is, gratifyingly, identical to the ANOVA table given in Heiberger and Holland’s book.
- (2 points) Given the result in your previous question, what technique would help you understand what is going on? Explain briefly.
The interaction is significant (which is our finding from the ANOVA), so a suitable technique would be to use simple effects to understand what that significant interaction actually means.
- (3 points) For strain
3DOk7, analyze the simple effect ofcomb. Hint: thestraincodes have an uppercase letter O in them, not a zero.)
Grab only the observations where strain is 3DOk7, and then do a one-way ANOVA of Npg as it depends on comb:
clover %>%
filter(strain == "3DOk7") -> d1
d1.1 <- aov(Npg ~ comb, data = d1)
summary(d1.1) Df Sum Sq Mean Sq F value Pr(>F)
comb 1 10.322 10.322 10.51 0.0118 *
Residuals 8 7.854 0.982
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant difference in nitrogen content between the two cultures of 3DOk7. To see which way around they are, look back at the interaction plot to see that the nitrogen content of the pure clover culture is lower than for the clover plus alfalfa culture. Or, if you like, run Tukey (which is otherwise pointless) to see the same thing:
TukeyHSD(d1.1) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Npg ~ comb, data = d1)
$comb
diff lwr upr p adj
clover+alfalfa-clover 2.031939 0.586832 3.477045 0.0118345The mean nitrogen content for clover plus alfalfa is higher than for just clover (by between 0.5 and 3.5 units).
Tukey is otherwise pointless here because there are only two levels of comb (two cultures of this strain), and you know they are significantly different. If you prefer, find another way to compare the means.
Extra: if the antibiosis hypothesis is correct, this one is significant but the wrong way around: that hypothesis says that clover should be higher than clover + alfalfa, and here it is not.
- (2 points) Repeat the previous question, but this time for the strain 3DOk1.
Now grab only the observations where strain is 3DOk1, and then do a one-way ANOVA of Npg as it depends on comb for those:
clover %>%
filter(strain == "3DOk1") -> d2
d2.1 <- aov(Npg ~ comb, data = d2)
summary(d2.1) Df Sum Sq Mean Sq F value Pr(>F)
comb 1 1.01 1.013 0.09 0.772
Residuals 8 90.12 11.265 This is not significant, and there is nothing more to do or say here.
- (2 points) Explain briefly why your results in the previous two questions are the kind of thing you would expect to see, given the result of the earlier ANOVA.
The ANOVA earlier found a significant interaction. This means that the effect of comb will be different for the different strains. In particular, we found from the simple effects analyses that there is no effect of comb for strain 3DOk1, but there is an effect for strain 3DOk7: that is, the effect of comb depends on which strain you are looking at.
Extra 1: you might think that strain 3DOk5 would have the biggest effect, based on what you saw on the interaction plot, but:
clover %>%
filter(strain == "3DOk5") -> d3
d3.1 <- aov(Npg ~ comb, data = d3)
summary(d3.1) Df Sum Sq Mean Sq F value Pr(>F)
comb 1 195.6 195.63 3.521 0.0974 .
Residuals 8 444.5 55.56
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1this big difference is actually not significant. This is because, as we saw in an earlier Extra, two of the clover observations for this strain are much bigger than the others, so it looks as if there is a lot of variability. (Compare the error sum of squares here, 444.5, with the values 7.9 and 90.12 from the simple effects analyses you did before.)
Extra 2: when I drew a grouped boxplot here, I said that there were only five observations per box, which made it hard to judge anything much (like outliers or equal spreads). I might have mentioned in lecture2 that you can get a better sense by treating the ANOVA as if it were a regression, by using lm instead of aov. The analysis gives the same results either way:
clover.2 <- lm(Npg ~ comb * strain, data = clover)
anova(clover.2)(you can check that the three P-values are the same.) You cannot do Tukey this way, but one thing a regression gives you are overall residuals, which you can assess the same way that you do for a regression:
library(broom)
clover.2 %>% augment(clover) -> clover.2a
ggplot(clover.2a, aes(x = .fitted, y = .resid)) + geom_point()
Each treatment group has its own fitted value (so there are only 12 different fitted values), but there seems to be fanning-out, in that the largest fitted values are also the most variable.
Also:
ggplot(clover.2a, aes(sample = .resid)) +
stat_qq() + stat_qq_line()
the distribution of residuals has long tails: there are three observations a long way off the line at the bottom, and three more a long way off the line at the top. So we have a couple of reasons to be suspicious of the results from this ANOVA.
You can also do this to plot residuals against strain:
ggplot(clover.2a, aes(y = .resid, x = strain)) + geom_point()
This partly loses the benefit of running the regression to get a large set of residuals to assess, but what I wanted to find out was where the unusual residuals were coming from. Bear in mind that we are looking for the usual “no pattern” picture that you want from a regression.
This tells you which strain each residual came from. Almost all of the unusual residuals come from strain 3DOk5 (five out of six), and the distribution of that strain is more spread out than the others. So equal spreads definitely fails because of that strain.
Extra 3: I mentioned that we don’t have a version of Mood’s median test that works with two-way ANOVA. Maybe we can make one. (This uses ideas from the log-linear models lecture in week 12.)
Let’s start by working out the grand median:
clover %>%
summarize(med = median(Npg)) %>%
pull(med) -> grand_median
grand_median[1] 25.12672Then I count up how many observations are bigger or smaller than the grand median in each treatment combination. Unfortunately, the obvious thing doesn’t work later:
clover %>%
mutate(bigger = (Npg > grand_median)) %>%
count(comb, strain, bigger)There should be 24 combinations that each have a count: six strains times two values of comb times two values of bigger, but there are only 16 here. This is because count, by default, omits any rows with zero counts. For example, all five observations of clover and 3DOk1 were above the grand median, and there were none below. I need (later) the zero for clover and 3DOk1 and FALSE. Unfortunately, getting it is a bit fiddly:
clover %>%
mutate(bigger = factor(Npg > grand_median)) %>%
count(comb, strain, bigger, .drop = FALSE) -> clover.tab
clover.tabNow I have 24 rows with the zero frequencies as well. What I had to do was:
- figure out where the zero frequencies would be (because I didn’t have both
TRUEandFALSEvalues forbigger) - make sure the offending column
biggerwas afactor(needed for the next thing) - run
countwith.dropset to FALSE. See this answer on Stack Overflow for the process behind these last two steps.
Now there is one row for each of the 24 combinations, and the zero frequencies are included.
When we learned Mood’s median test, we had two categorical variables (the grouping variable, plus an indicator of whether each value was above or below the grand median). We made the test by testing for an association between the two categorical variables using a chi-squared test (which you had either seen before or could learn about without too much trouble).
Now, though, we have three categorical variables, comb, strain, and bigger. We can’t use a simple chi-squared test any more because there are more than two categorical variables. Instead, we use a “log-linear model”, which we learn about properly in week 12. This is a generalization of the chi-squared test to find associations between any number of categorical variables. It is a generalized linear model with Poisson family, and we start by fitting the interaction between all three of our categorical variables. Then we use drop1 to see what if anything we can remove:3
clover.tab.1 <- glm(n ~ comb * strain * bigger,
family = "poisson",
data = clover.tab)
drop1(clover.tab.1, test = "Chisq")The only thing up for grabs is that three-way interaction. It is not significant, so it can come out.
From here, you can do as we did with regression: proceed one step at a time by removing the least significant thing and re-fitting, or you can use step. I’m going to be lazy and use step (you’ll see the other way in week 12):
clover.tab.2 <- step(clover.tab.1)Start: AIC=97.93
n ~ comb * strain * bigger
Df Deviance AIC
- comb:strain:bigger 5 3.9051 91.832
<none> 0.0000 97.927
Step: AIC=91.83
n ~ comb + strain + bigger + comb:strain + comb:bigger + strain:bigger
Df Deviance AIC
- comb:strain 5 5.915 83.842
<none> 3.905 91.832
- comb:bigger 1 6.985 92.912
- strain:bigger 5 56.913 134.841
Step: AIC=83.84
n ~ comb + strain + bigger + comb:bigger + strain:bigger
Df Deviance AIC
- comb:bigger 1 6.985 82.912
<none> 5.915 83.842
- strain:bigger 5 56.913 124.841
Step: AIC=82.91
n ~ comb + strain + bigger + strain:bigger
Df Deviance AIC
- comb 1 6.985 80.912
<none> 6.985 82.912
- strain:bigger 5 57.983 123.911
Step: AIC=80.91
n ~ strain + bigger + strain:bigger
Df Deviance AIC
<none> 6.985 80.912
- strain:bigger 5 57.983 121.911drop1(clover.tab.2, test = "Chisq")What remains is a (very) significant interaction between strain and bigger. That means that there is an association between the two variables, and thus the strains differ in terms of median nitrogen content. (This is the same logic as Mood’s median test.) Note that, according to this model, there is no interaction between strain and comb (if there had been, the model’s three-way interaction comb:strain:bigger would have been significant in the model) and also no comb main effect (otherwise comb:bigger would have been significant in the model).
According to the week 12 notes, the way to investigate a significant association is to draw a graph like this one:
ggplot(clover.tab, aes(x = strain, y = n, fill = bigger)) +
geom_col(position = "fill")
This is a (stacked) bar chart, but with all the bars scaled to the same height, to make it easy to compare how much red or blue there is in each bar. The more blue, the more plants of that strain had nitrogen content above the grand median (so the median was high); the more red, the more plants below the grand median (median was low). So, strains 3DOk1 and 3DOk5 had the highest median nitrogen, and strains 3DOk13 and 3DOk4 have the lowest.
Note that comb has completely disappeared from the model. This means that comb has no effect of any kind on nitrogen fixation. You’ll remember that the research question that went with these data said that if comb was clover+alfalfa, nitrogen would be less than if comb was clover. Our Mood’s-Median-Test-style analysis has found no evidence in favour of this research hypothesis at all. And, because the analysis is based only on counts above and below, those suspected outliers will not cause any problems at all. I think this is a more reliable conclusion than the interaction we found that appeared to be mostly based on those very high values for strain 3DOk5.
Now that we know that the interesting effect is that of strain, and we have the actual nitrogen values, we can draw an ordinary boxplot of Npg and strain, with, we would expect, the same kind of conclusions as above:
ggplot(clover, aes(x = strain, y = Npg)) + geom_boxplot()
which says the same thing as the bar chart did. This is our finding. The researchers may not like it, but this is what we found out.
Points: \(13 + 17 = 30\).
Footnotes
I mentioned “three observations similar and two different ones” above. If the two different ones are both at the same end, as here, you get a big box, but if the two different ones are one higher and one lower, you get a small box with outliers both ends, as for
3DOk7-cloveron the grouped boxplot. This is because, with five observations, Q3 is the second-biggest observation and Q1 is the second-smallest. Such are the hazards of boxplots based on small numbers of observations.↩︎I am writing this before the lecture, and I haven’t looked at my slides for this recently.↩︎
This model requires a dataframe with a frequency for all combinations of our categorical variables’ levels. It is in order to make this work that we went to all that trouble with zero frequencies above.↩︎