STAD29 Assignment 6

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

library(tidyverse)
library(car)

Facts and opinions

5,035 Americans were each given ten news statements, and were asked to classify each news statement as “fact” or “opinion”. The total number of news statements correctly classified was recorded, along with the age group of each respondent. The purpose of the survey was to detect whether there was a difference between older and younger Americans in their ability to distinguish facts from opinions. There are five age groups: 18–27, 28–38, 39–49, 50–65, and 65+. Please ignore the fact that the second-last age group should really be 50–64. The data are in http://datafiles.ritsokiguess.site/fact-opinion.csv.

(1 point) Read in and display (a little of) the data.

I hope you can do this from memory by now:

my_url <- "http://datafiles.ritsokiguess.site/fact-opinion.csv"
survey <- read_csv(my_url)

Rows: 5035 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): age
dbl (1): correct

ℹ 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.

survey

The columns are called age and correct.

Extra: this is based on the fact_opinion data set in the openintro package, which is itself based on a 2018 survey by Pew Research.

(1 point) Create and save a dataframe in which the age column is a factor (not just text).

There is no problem with overwriting the old age column and the dataframe (whatever you called it):

survey %>% mutate(age = factor(age)) -> survey
survey

although there is also no problem with creating a new column and/or a new dataframe should you prefer to. It is smart to display the new dataframe to check that the new column you created is indeed a fctr.

(4 points) The people who did the survey were interested in comparing the following age groups:

18–27 vs. 28–38
50–65 vs. 65+
the average of (18–27 and 28–38) vs 39–49
the average of the three youngest age groups vs. the average of the two oldest age groups.

Set up contrasts that will enable you to test those comparisons (the tests are coming later). You may assume that the age groups are in numerical order, 18–27 first and 65+ last.

I listed the contrasts in ascending order of difficulty to construct them. You have some alternatives (discussed below), and use names that will remind you later of what you are actually comparing:

youngest <- c(1, -1, 0, 0, 0)
oldest <- c(0, 0, 0, 1, -1)
young_vs_middle <- c(1/2, 1/2, -1, 0, 0)
young_vs_old <- c(1/3, 1/3, 1/3, -1/2, -1/2)

Your first quick check should be that the five numbers within each contrast add up to zero (this is actually the definition of a contrast). The second thing to worry about is that the number on the bottom of a fraction (in the last two contrasts) represents the number of means in a group you are comparing. For example, for the last one, we are comparing the three youngest groups (so the weight for each of these should be 1/3) against the two oldest groups (so the weight for each of those should be 1/2 with the opposite sign).

The order of the numbers in each contrast is the same as the order of the age groups that you found with count above: that is, the first number corresponds to the 18–27 age group, and so if it zero, it means that the 18–27 age group does not figure in that contrast.

Alternatives: you can multiply a contrast through by any constant, which could be $-1$ (which will flip all the signs), or by something that will make the weights whole numbers (eg., 6 for the last one, so that the weights become 2 three times followed by $-3$ twice). Doing this won’t stop any contrast from being a contrast (the five weights will still add up to zero) and won’t change the orthogonal-ness of any contrasts (coming up).

(1 point) Pick any two of your contrasts and verify that they are orthogonal.

Multiply the weights in the contrasts together elementwise. I’ll do my last two contrasts:

young_vs_middle * young_vs_old

[1]  0.1666667  0.1666667 -0.3333333  0.0000000  0.0000000

and verify that the results add up to zero, or do that in one step (equally good):

sum(young_vs_middle * young_vs_old)

[1] 0

Extra: in principle, you check every pair of contrasts for orthogonality, but with four contrasts, that is a lot of work (there are six pairs of contrasts).

(3 points) Test all four of your contrasts for significance. (Interpretation is coming up in a moment.)

Three steps:

glue the contrasts together into a matrix
set that as contrasts for your factor column
run the ANOVA as if it is a regression:

m <- cbind(youngest, oldest, young_vs_middle, young_vs_old)
contrasts(survey$age) <- m
survey.1 <- lm(correct ~ age, data = survey)

You can display the summary of survey.1 now, or in the next question. Either is good.

Extra: you can check your matrix of contrasts, if you wish:

     youngest oldest young_vs_middle young_vs_old
[1,]        1      0             0.5    0.3333333
[2,]       -1      0             0.5    0.3333333
[3,]        0      0            -1.0    0.3333333
[4,]        0      1             0.0   -0.5000000
[5,]        0     -1             0.0   -0.5000000

Each column has the name of its contrast, and the contrast weights can be read down the columns. The purpose of the contrasts line is to set lm up so that the “regression slopes” become the contrasts we want to test.

A more mathematical way to check for orthogonality is to work out m-transpose matrix-multiplied by m:¹

t(m) %*% m

                youngest oldest young_vs_middle young_vs_old
youngest               2      0             0.0    0.0000000
oldest                 0      2             0.0    0.0000000
young_vs_middle        0      0             1.5    0.0000000
young_vs_old           0      0             0.0    0.8333333

This is a shorthand way of checking all the contrasts for orthogonality with each other. Ignore the numbers down the diagonal (top-left to bottom-right), since a contrast is not going to be orthogonal with itself; as long as all the off-diagonal numbers are zero, your contrasts are orthogonal.

If you have done linear algebra, you might recognize this sort of thing from there. The reason that orthogonality is important is that each contrast is then tested independently from all the others.

(2 points) What do you conclude from the analysis, in the context of the data?

Here is where you need to look at the summary of the lm, if you haven’t done so already:

summary(survey.1)


Call:
lm(formula = correct ~ age, data = survey)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.5186 -1.6635  0.3104  1.5219  3.3365 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         7.171768   0.032520 220.537   <2e-16 ***
ageyoungest         0.004817   0.057243   0.084    0.933    
ageoldest          -0.013085   0.039719  -0.329    0.742    
ageyoung_vs_middle  0.023785   0.067651   0.352    0.725    
ageyoung_vs_old     0.990442   0.074109  13.365   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.157 on 5030 degrees of freedom
Multiple R-squared:  0.03439,   Adjusted R-squared:  0.03362 
F-statistic: 44.78 on 4 and 5030 DF,  p-value: < 2.2e-16

To take our four contrasts in order, and translating back to what we were actually testing:

the two youngest age groups do not differ significantly (in ability to detect fact from opinion in news statements)
the two oldest age groups do not differ significantly (ditto)
the two youngest age groups (taken together) do not differ significantly from the 39–49 age group
the three youngest age groups do differ significantly (in detecting fact from opinion in news stories) from the two oldest age groups.

Extra: That’s as far as I need you to go, but you might be curious about how the oldest and youngest age groups differ. A boxplot would be a good way to see that:

ggplot(survey, aes(x = age, y = correct)) + geom_boxplot()

There seems to be a sharp dividing line at age 50, with respondents older than that doing a worse job on average at distinguishing fact from opinion than younger respondents. The boxplot, in fact, tells exactly the same story as the contrasts, but the contrasts come with P-values attached. The sample size is large (over 600 observations in each group), so that even though the above-50 vs. below-50 difference is small, it is very definitely significant.

It is just coincidence that in the examples of contrasts we have seen, only one of them has come out significant each time. Any number of the tests could, in principle, come out significant.

Spanish onions

In an experiment involving the production of Spanish onions in two locations in South Australia, onions were grown in 84 plots. For each plot, the yield of onions grown (in grams per plant) was recorded, as well as the location of the plot. It was also suspected that the yield of onions grown would depend on how close together the plants were planted. This is recorded in column dens. The data are in http://datafiles.ritsokiguess.site/onions.csv.

(1 point) Read in and display (some of) the data.

As always:

my_url <- "http://datafiles.ritsokiguess.site/onions.csv"
onions <- read_csv(my_url)

Rows: 84 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): location
dbl (2): dens, yield

ℹ 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.

onions

These data are the onions dataset from package SemiPar. The locations were coded as 0 and 1, so for you I changed these to the actual names of the locations.

(2 points) By drawing a suitable graph, demonstrate that the relationship between yield and planting density is non-linear for each location.

This is no more than a scatterplot with the locations distinguished by colour. However, since we are trying to demonstrate a non-linear relationship, we do not want to draw regression lines here; the default geom_smooth is the right thing to add to the graph here:

ggplot(onions, aes(x = dens, y = yield, 
                   colour = location)) +
  geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

These are downward but non-linear trends: they descend quickly at first and then level off.

You might say “well, why not draw the regression lines and show that the points don’t fit them very well?” There is some value in this but, as I mention below, adding regression lines to the plot will make you think that lines fit better than they actually do. It is better to “let the data decide” what trend there is, and then it will be clear what the shape of the trends is. If that trend happens to be close to linear, geom_smooth() will show that, but if not, geom_smooth() will clearly show that as well.

At this point, what I am looking for you to do is think: ask yourself how you can best show that the trends are not linear. (I gave you a hint in the question that this is what you are looking for.)

(2 points) Because of the non-linearity, the researchers decided to predict one over the square root of yield instead of yield itself. Support that decision. The best answer demonstrates how the researchers came up with “one over square root” rather than something else.

The best way to support that decision is to show that you also think that one over the square root of yield is the best transformation to make. As we learned in C32, Box-Cox is a good way to come up with a good transformation:

library(MASS, exclude = "select")


Attaching package: 'MASS'

The following object is masked _by_ '.GlobalEnv':

    survey

boxcox(yield ~ dens * location, data = onions)

Inside the boxcox, put the best model that you know of at the moment (including the interaction, because you don’t know yet whether it can be removed). (That said, the decision is the same if you leave it out on this occasion, but the right thing to do is to put the interaction in “just in case”.)

The best transformation appears to be power $-0.5$, that is, one over square root, which is in the confidence interval, supporting what the experimenters said.

The second-best way to do it is to re-draw your plot of the previous question, but with yield replaced with one over the square root of yield, and demonstrate that the trends are now straight, or close to it. Once again, use a default smooth rather than a regression line (because we want to see whether the trends are straight, not to assume that they are at this point).

ggplot(onions, aes(x = dens, y = 1 / sqrt(yield), 
                   colour = location)) +
  geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

The trends are at least straighter than they were before, but this does not rule out the possibility that a better transformation could be found (whereas the Box-Cox narrows things down a lot, ruling out the log and reciprocal transformations).

(3 points) Fit a model using your transformed response, and demonstrate that the interaction term needs to remain in the model.

You have a couple of options here: either recall that sqrt is square root in R, and work out one over that:

onions.1 <- lm(1 / sqrt(yield) ~ dens * location, data = onions)

or raise yield to a negative power using ^, but then you have to remember that ^ has a special meaning in model formulas, so you need to wrap the power in I() to work around that, as in C32 when we were adding a squared term:²

onions.1 <- lm(I(yield^(-0.5)) ~ dens * location, data = onions)

Then, to see whether you can get rid of the interaction, the best way is this (with an $F$-test, because it is an ordinary regression):

drop1(onions.1, test = "F")

and we see that the interaction, with a P-value of $4.7 \times 10^{-8}$, much less than 0.05, has to stay.

The second-best way is to use summary:

summary(onions.1)


Call:
lm(formula = I(yield^(-0.5)) ~ dens * location, data = onions)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0121384 -0.0031970 -0.0006592  0.0034537  0.0128670 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.673e-02  1.781e-03  31.859  < 2e-16 ***
dens                  4.667e-04  1.992e-05  23.425  < 2e-16 ***
locationVirginia      4.001e-03  2.412e-03   1.659    0.101    
dens:locationVirginia 1.734e-04  2.873e-05   6.037 4.67e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.005379 on 80 degrees of freedom
Multiple R-squared:  0.9519,    Adjusted R-squared:  0.9501 
F-statistic: 527.6 on 3 and 80 DF,  p-value: < 2.2e-16

and it happens that the P-value for the interaction is the same. It worked in this case because there were only two locations, and so comparing the slope for Virginia with that of the baseline location Purnong Landing is actually doing the same test. (If there were three locations, it would be a different story.)

If you use summary, you also need to explain why it will work (that there are two locations and not more).

(3 points) By drawing another graph, explain briefly why the significant interaction is not surprising.

The graph you need is like the first one you drew, but with the transformed yield instead of the yield itself. Also, it is now reasonable to assume straight-line relationships, because the purpose of a transformation is to make the trends linear, or closer to it, or that we have fitted a straight-line model and that is what we are investigating now. Hence:

ggplot(onions, aes(x = dens, y = 1/sqrt(yield), colour = factor(location))) +
  geom_point() + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

Using yield to a negative power also works, but I would put the power in brackets to make sure the negative is done before the raising to a power. There is no need for I() here, because we are no longer inside a modelling function such as lm:

ggplot(onions, aes(x = dens, y = yield^(-0.5), colour = factor(location))) +
  geom_point() + geom_smooth(method = "lm")

`geom_smooth()` using formula = 'y ~ x'

The interaction was significant because those lines are not parallel.

Extra: The lines may look close to parallel, but the scatter of points around the lines is relatively small. Especially on the right, the blue points are all close to the blue line, and the red points close to the red one, and hence we can be pretty sure that the blue line goes up faster than the red one, and thus the P-value for the interaction comes out really small.

The linear trends do look pretty good here, but bear in mind that putting a line on the plot tends to bias you towards thinking that a line fits well. For comparison, look at my second-best answer to Question 9; the smooth trend for Purnong Landing looks more or less linear, but the one for Virginia seems to increase more sharply at the beginning and end and not so fast in the middle.

The last thing I would typically ask you on a question like this is something about which location is better, and whether that depends on planting density. The confusing thing here is that with the transformation we used here, a high value of yield transforms to a low value, with the result that the graph here goes uphill, but the one we drew first goes downhill:

ggplot(onions, aes(x = dens, y = yield, 
                   colour = location)) +
  geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

On the original scale, the better yield at all planting densities is from Purnong Landing. The biggest difference seems to be at a dens of about 50, and the difference elsewhere is smaller.

The other thing to say is that both these trends are downhill (on the original scale), and this is telling us that the more close together you plant the onions (higher dens), the less the yield is. Remember that the yield is measured per plant, so there is presumably an optimal density of plants that maximizes the total yield over a plot of land.

In the original source, dens was measured as plants per square metre, so you can work out the total yield per square metre by multiplying yield by dens. This is the kind of thing a farmer would be interested in: they have a certain number of hectares of land, and they want to know how close together to plant the onions so that they will end up with the most total onions to sell at the end of the growing season. We can investigate that:

onions %>% 
  mutate(yield_per_sqm = yield * dens) -> onions
ggplot(onions, aes(x = dens, y = yield_per_sqm, colour = location)) +
  geom_point() + geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

On this scale, there seems to be more uncertainty. At Virginia, the highest yield per square metre seems to be around a density of 120 plants per square metre, but at Purnong Landing, the yield per square metre seems to keep going up. Another reason we should not be too confident about inferring anything from a graph like this is that most of the observed plant densities were small, say around 50; if we had more observations at a dens value around 150, we would get a clearer picture of whether the yield per square metre was actually coming down again. From a practical point of view, you would imagine that if you had too many onion plants planted too close together, they would be competing for light, rain, fertilizer etc., and most of them would not grow very well.³ But I don’t think that’s what this experiment was about.

Uniformity trial of wheat

In a uniformity trial, something (in this case wheat) is grown over a large number of equally-sized fields, and the idea is that the amount of wheat grown should be about the same (ie., “uniform”) in all the fields. In this particular trial, wheat is grown over 25 fields arranged in a square array, with rows labelled A through E and columns labelled J through M. In this trial, two variables were measured:

the amount of grain, in pounds. This is the edible (by humans) part of wheat. (For example, it can be made into flour for baking.)
the amount of straw, in pounds. This is the dry stalks of wheat plants after the grain has been removed. It is used as bedding or feed for animals such as cattle or horses.

Twenty observations were obtained from each field (by subdividing each of the 25 large fields into 20 smaller ones). The data are in http://datafiles.ritsokiguess.site/mercer-wheat.csv.

(1 point) Read in and display (some of) the data.

No surprises here:

my_url <- "http://datafiles.ritsokiguess.site/mercer-wheat.csv"
wheat <- read_csv(my_url)

Rows: 500 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): row, col
dbl (2): grain, straw

ℹ 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.

wheat

There are indeed $25 \times 20 = 500$ observations.

Extra: These data are based on the mercer.wheat.uniformity data from the agridat package,⁴ and those data came from a trial done at Rothamsted Experimental Station in 1910 (!). This is where the famous Fisher worked for many years, and developed a lot of his ideas about analysis of variance.

(2 points) Why would it be appropriate to run a multivariate analysis of variance here?

There are two (quantitative) response variables (the two variables that were measured), namely grain and straw.

Make sure you say what the two response variables are here, so that it is clear that you know. Only one point for “there are two response variables” without naming them or otherwise making it clear which variables are responses.

(2 points) Run a suitable multivariate analysis of variance, and display the results.

Make a response variable by cbind-ing the two quantitative response variables, either first or in the manova, and then display the summary:

wheat.1  <-  manova(cbind(grain, straw) ~ row * col, data = wheat)
summary(wheat.1)

           Df  Pillai approx F num Df den Df    Pr(>F)    
row         4 0.22035   14.704      8    950 < 2.2e-16 ***
col         4 0.21046   13.966      8    950 < 2.2e-16 ***
row:col    16 0.51195   10.214     32    950 < 2.2e-16 ***
Residuals 475                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

If you want to try using Manova, you can:

wheat.2 <- lm(cbind(grain, straw) ~ row * col, data = wheat)
Manova(wheat.2)


Type II MANOVA Tests: Pillai test statistic
        Df test stat approx F num Df den Df    Pr(>F)    
row      4   0.22035   14.704      8    950 < 2.2e-16 ***
col      4   0.21046   13.966      8    950 < 2.2e-16 ***
row:col 16   0.51195   10.214     32    950 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The strategy here is to fit a “multivariate regression” using lm with your two-column response variable, and then to display the Manova (uppercase M) of that.⁵ The result, including the test statistics and degrees of freedom, is the same, so either way is good.

Extra: I didn’t ask you to run Box’s M test. If you want to check for the multivariate version of equal spreads, you can, borrowing the example from the peanuts data:

library(MVTests)


Attaching package: 'MVTests'

The following object is masked from 'package:datasets':

    iris

wheat %>% mutate(combo = str_c(row, "-", col)) -> wheatx
response <- with(wheat, cbind(grain, straw))
summary(BoxM(response, wheatx$combo))

       Box's M Test 

Chi-Squared Value = 75.79796 , df = 72  and p-value: 0.357

Absolutely no problem here.

(2 points) What do you conclude from your multivariate analysis of variance, in the context of the data?

The interaction is significant, so one or both of the amount of grain or straw depends on the combination of row and col in some way.

This is about as far as you can go so far, and once again, it is an error to try to interpret main effects in the presence of a significant interaction. You’ll also note that the conclusion so far is very vague.

(2 points) By drawing a suitable graph, add some insight to your conclusion of the previous question. (Hint: if you need something from beyond the course materials to answer this, cite your trustworthy source.)

We are in the happy position of having two response variables rather than more, so we can make a non-three-dimensional graph. Let’s summarize our variables:

grain: quantitative, response
straw: quantitative, response
row: categorical, explanatory
col: categorical, explanatory.

If we only had one categorical variable, we could draw a scatterplot with points identified by colour. I think the best way here is to find out how to distinguish the other categorical variable by something else other than colour, such as plotting symbol (that is, plot the points using different symbols rather than just circles).

I wasn’t quite sure what to search for. ggplot plotting symbol gave me all sorts of things, including how to label points (not really the same idea), but one of the search results was this one, which seems like a good source. I think this site is a good source because it has all sorts of tips on statistics and data handling (a good sign). I looked at what they had to say about survival analysis, and it says many of the same things we said in our lecture notes. Being able to look at other information from the same source is a great way to decide whether you trust a source (see, for instance, the Crash Course series of videos on this), but using an LLM is appallingly bad in comparison because you have no idea where the information came from.

Anyway: Scroll down to “multiple groups” in the page I found (that’s what we have) to see that you can use shape inside aes, so we can try this:

ggplot(wheat, aes(x = grain, y = straw, colour = row, shape = col)) + geom_point()

The points at the bottom right seem to be different from the others (off the trend of the others), and they are all light blue plus signs. Look at the legends to the right to see that they are all from row D (light blue) and column M (plus signs). So the field in this particular row and column seems to be different from the others. Because it’s a combination of row and column, this shows up as a significant interaction in the analysis. It’s not true, for example, that row D is always different from the other rows (which would be a row D main effect); it’s only different when combined with column M.

You have a lot of freedom here in making this graph. You can put grain and straw on the other axes, or you can have col as colour and row as shape. However, because it’s the fourth row and fourth column that are different, your points away from the trend should be light blue plus signs.⁶

You also have freedom in the type of graph you choose. One way to go that stays within what we have learned in this course is to consider that we have an “extra” categorical variable and use facets for it, like this (for example):

ggplot(wheat, aes(x = grain, y = straw, colour = row)) +
  geom_point() + facet_wrap(~ col)

These are all regular scatters of points except that for column M, there is a collection of light blue points down and to the right (row D), so the field in row D and column M is once again different from the rest.

Or you might consider facetting on both categorical variables. Since there are now two of them, you need to use facet_grid rather than facet_wrap, for example like this:⁷

ggplot(wheat, aes(x = grain, y = straw)) +
  geom_point() + facet_grid(row ~ col)

Now you have to cast your eye over the 25 subplots and find ones where the trend is somehow different. Most of the trends are more or less linear and arranged towards the top left of their facets, but the trend for row D and column M is towards the bottom right of its facet, so that again this combination of row and column is different from the others.⁸

My take is that you have to work harder with this graph to see the different row and column, so this choice of graph is not quite as good as the previous two.

Extra: here is a different take, one that might actually work for more than two responses (if they are on a similar scale numerically): make boxplots of grain and straw for each combination of row and col (plotted in facets as in the previous graph). To do that, we need the grain and straw values in one column, with a second column labelling which one each number is a measurement of:

wheat %>% 
  pivot_longer(grain:straw, 
               names_to = "yname", 
               values_to = "yval")

This, you will note, is the same kind of idea we used for plotting residuals against all the explanatory variables in regression, but now the two variables are both responses, so I called the new columns yname and yval. Then, make boxplots of yval for each yname:

wheat %>% 
  pivot_longer(grain:straw, 
               names_to = "yname", 
               values_to = "yval") %>%
  ggplot(aes(x = yname, y = yval)) + geom_boxplot()

and then do this for each combination of row and column as we did it above:

wheat %>% 
  pivot_longer(grain:straw, 
               names_to = "yname", 
               values_to = "yval") %>%
  ggplot(aes(x = yname, y = yval)) + geom_boxplot() +
  facet_grid(row ~ col)

Now, what you see is that there is more straw than grain produced in all the fields, and by about the same amount, except for the one field that is in row D and column M, where there is more grain and less straw produced than in the other fields, so that the amounts of grain and straw are almost equal in this field. (Remember that MANOVA is a comparison of means, so that you are looking for a field that is different in terms of, well, medians rather than thinking about outliers or spreads or anything like that.)

This kind of graph will work if all the response variables are the same kind of numbers (as here, about 4 for grain and about 6 for straw) so that any kind of deviation from the overall pattern will show up. It would not be so good if one response was about 4 and the other one was about 800. For a graph, it might be useful in that case to standardize the response variables.

Extra extra: I have to come clean and admit that I faked up these data (a little bit, to give you something to do). In the original dataset (that I mentioned above), the MANOVA has nothing significant, so I took out the interaction and there was still nothing significant. This is actually a success for a uniformity trial: the amount of grain and straw, and their combination, originally showed no significant difference across rows and columns.⁹

I thought that was rather dull for an assignment question, so I went back in and, just for row D and column M, added 1 to the grain values and subtracted 2 from the straw values. This was enough to make the interaction significant, and to make it pretty clear, whichever of the above graphs you made, that it was row D and column M that was different. (That is to say, making this change was literally what produced the significant interaction.)

Thinking about ways in which a uniformity trial could actually fail: you might expect something like a “fertility gradient”, something like (a) both grain and straw being less at the bottom left and more at the top right, in a continuous trend, or (b) the relationship between grain and straw changing in a similar kind of way, so that some fields produce more grain and some of them more straw. Those kinds of things would have been more realistic, but I thought having only one of those fields be different from all the others was something I could expect you to be able to detect.

Points: $$.

Footnotes

t() applied to a matrix works out its transpose, and %*% does matrix multiplication. If you were to use * here, it would multiply the matrices elementwise, which is not what you want.↩︎
It looks as if it works here without the I(), but this is not something to rely on.↩︎
It is also possible that a farmer would do what is called “thinning out”: plant a lot of plants to start with, and then pull out the ones that are not growing very well.↩︎
I fiddled with some things, so it is not quite the original data.↩︎
This is why I loaded the car package earlier.↩︎
I could probably annoy one of our number by calling these “Manchester City plus signs”.↩︎
It would be an error to use scales = "free" here, because that would hide the very thing you are looking for.↩︎
I did the facetting this way around so that the cols are actually columns on the page, so this is like an aerial view of the fields as they were originally labelled.↩︎
Apart from that “absence of evidence is not evidence of absence”. There might have actually been a difference that my analysis didn’t find, a type II error, or the fields at Rothamsted are not the same as your typical farmer’s field.↩︎