# Which model to choose for ANOVA with tutors and assignments?

This semester I am the teaching assistant for lecture with graded homework assignments. There are around 100 students, six assignments and nine tutors grading the submissions. Each homework gives the same number of points such that they all have the same weight.

Several years ago each tutor would grade the homework of the students in their tutorials. We realized that this likely is unfair as tutors have different levels of strictness. Therefore we randomly assign the submissions to the tutors for each assignment. With more tutors than assignments the chances should be good that they get a different tutor for each assignment.

Now I would like to see whether this method is sufficiently fair. My concrete questions are:

1. Is there a measurable difference between the grading styles of the different tutors?
2. Do all students on average have a fair result or did students get graded more strict than others?

For the first question I have figured that the ANOVA is the sensible method. I have started with the simplest model: Points ~ Tutor. This is the syntax of R's lm function and means that the dependent variable Points can be explained by the independent variable Tutor.

The output of the ANOVA is the following:

Analysis of Variance Table

Response: Points
Df  Sum Sq Mean Sq F value    Pr(>F)
Tutor       8  145.89  18.236  6.5267 4.646e-08 ***
Residuals 447 1248.94   2.794
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


From the P-value alone, this already is significant and tells me that the tutor has a significant impact on the points given.

Also looking at a boxplot suggests that there are differences:

However, the Sum Sq column has a very large unexplained residual. From this I deduce that I could do better and the model does not really fit the data well. Looking at the Mean Sq however, the residuals are not so bad in comparison. Which of the two numbers do I have to look at? In the infinite statistics limit, I expect only the Mean Sq to be meaningful as the Sum Sq of the residuals will just grow linearly.

It could be possible that the tutors that appear strict just have been assigned the weak students by chance. If one takes into account that students have different levels intrinsically, then perhaps the effect of the tutor is reduced. So I have added the student as an independent variable as well such that the model reads Points ~ Tutor + Student. The output now is this:

           Df Sum Sq Mean Sq F value    Pr(>F)
Tutor       8 145.89 18.2358 11.5730 1.299e-14 ***
Student    96 695.86  7.2485  4.6001 < 2.2e-16 ***
Residuals 351 553.08  1.5757


The sum of squares for the Tutor variable has not changed, but it's F-value has increased, making this more significant now. Also the residuals have reduced a lot, so this seems like a good addition to the model?

Does the lower P-value for the Student dependence compared to the Tutor dependence tell me that it more depends on the student than on the tutor, an indication that the grading procedure is fair?

The assignments did not all have the same difficulty, so perhaps this has also an influence. Doing Points ~ Tutor + Student + Assignment gives me this:

            Df Sum Sq Mean Sq F value    Pr(>F)
Tutor        8 145.89 18.2358 12.6994 4.964e-16 ***
Student     96 695.86  7.2485  5.0479 < 2.2e-16 ***
Assignment   4  54.80 13.7007  9.5412 2.497e-07 ***
Residuals  347 498.28  1.4360


Apparently the assignments were also different enough for this to get a low P-value. There are only four degrees of freedom as only five assignments have been fully graded up to this point.

Just out of curiosity I have added the interaction effect of tutor and assignment. Perhaps each tutor interprets each assignment a bit different. That is the result:

                  Df Sum Sq Mean Sq F value    Pr(>F)
Tutor              8 145.89 18.2358 14.3168 < 2.2e-16 ***
Student           96 695.86  7.2485  5.6907 < 2.2e-16 ***
Assignment         4  54.80 13.7007 10.7563 3.471e-08 ***
Tutor:Assignment  31  95.78  3.0895  2.4256 6.637e-05 ***
Residuals        316 402.50  1.2737


All of these things have extremely low P-values, which starts to sound dubious.

My background is Physics, not Psychology. Therefore I am not so sure which of these models actually is the one to choose and how to interpret is. Which one best matches my actual question at hand? And how would I interpret these P-values?

• Independent of your statistical model have you considered graphing your data e.g. in a boxplot to look at the actual differences? Even if your p value is below 0.01 or whatever other cut-off level you use the actual size of the difference may be pretty small (but statistically significant nonetheless) Commented Jul 14, 2018 at 12:24
• @MaartenPunt: Yes, I had tried to do boxplots of the points with respect to the tutors and assignments. The plot is added in the question now. It seems that there are tutors who have given less points. But the plot does not show whether they just have been assigned weak students by chance. Commented Jul 18, 2018 at 8:01
• Nice. I see the differences, although they are not particularly large if you ask me, except tutor no. 7. I've been chewing on the weak student thing. What I think may be a solution for visualization is the following: For each student tutor combination calculate the average grade of that student on the assignments NOT graded by that tutor. Then classify these average grades into deciles (or another categorization of weak to strong). Then create a bar graph or some other frequency chart how often tutors got assigned students that were weak according to other tutors. Hope that helps Commented Jul 19, 2018 at 12:14
• The low p-values don't seem dubious to me. You are tapping in to the factors affecting student grades. It seems to me that you need to include all of Tutor, Assignment, and Student, for this model in this context. And I would include interactions as make sense and improve the model. You might also look in to estimated marginal means (emmeans package). I'd have to think about how to do this sensibly, but the EM mean for Student would give their average score had the other factors in the model been balanced for that student. Commented Jul 19, 2018 at 14:51

For model selection I suggest running pairwise comparisons using R function anova(), which will tell you if the models differ significantly in the amount of variance they explain. Alternatively I'd compare their AIC.

Did students take more than one test each? If so it makes sense to add them as explanatory variable and that should leave almost no residual unexplained variance.

Then again:

1-you know that different students perform differently and that's not really your focus;

2-multiple test results from the same student are not independent;

so I'd run a mixed effect model with tutor and assignment (and their interaction?) as fixed variables and student ID as random variable. This way you'll obtain an estimate of the effect of the fixed variables on the marks, taking into account that some of the remaining "background noise" is simply due to students performing differently. This should work:

library(lmerTest) # does what lme4 does, plus pseudo-p-values
lmm.interaction <- lmer(test.score~tutor*assignment+(1|student.ID), data=yourdata)
lmm.onlytutor <- lmer(test.score~tutor+(1|student.ID), data=yourdata)
lmm.onlyassignment <- lmer(test.score~assignment+(1|student.ID), data=yourdata)
lmm.null <- lmer(test.score~1+(1|student.ID), data=yourdata)
# pairwise comparisons with anova()
# see which one has the lowest AIC: