# What regression model to use for judgment data?

I am new to statistic modeling and I need help with choosing the right model for an analysis.

We run an experiment on word order variation in German. We presented sentences to 27 participants who had to rate the acceptability of the sentences on a 7-point scale. The sentences all contained three arguments: subject (S), direct object (DO) and indirect object (IO) with all possible their 8 permutations. They also varied with respect to whether each argument was a full noun phrase (like “the dog”) or a pronoun (“he”), so that there was a total of 8 x 6 = 48 critical items. The data looks like this:

Item WordOrder PhraseTyp   Rating Participant
1    S-IO-DO  Pro-NP-Pro    3      A
2    S-IO-DO  NP-NP-Pro     2      A
3    S-IO-DO  Pro-Pro-Pro   2      A
4    S-IO-DO  NP-NP-NP      4      A
…      …         …          …      …
15   S-DO-IO  Pro-NP-NP     6      A
16   S-DO-IO  NP-Pro-Pro    7      A
…      …         …          …      …


What I would like to test is if the word order of the argument (“WordOrder”) or the distribution of pronouns/nominal phrases (“PhraseType”) has an effect on the rating.

I guess I can use some kind of regression modeling. I am not sure, what I should use. According to what I read, I seems as if I need some kind of mixed effects model, with random effects for participants (not sure about items). Is that correct? I am also not sure if I need linear regression or ordinal logistic regression? What is the best function to use with R? I tried the following:

glmer(Rating ~ WordOrder + PhraseType + (1|Participant), data = RatingData)


This gives me an warning message:

In glmer(Rating ~ WordOrder + PhraseType + (1 | Participant), data = RatingData): calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly


If I use this, I get the following summary of the the model:

Linear mixed model fit by REML ['lmerMod']
Formula: Rating ~ WordOrder + PhraseType + (1 | Participant)
Data: RatingData

REML criterion at convergence: 4292.9

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.4997 -0.6709 -0.0844  0.5890  3.5763

Random effects:
Groups      Name        Variance Std.Dev.
Participant (Intercept) 0.3809   0.6172
Residual                1.7083   1.3070
Number of obs: 1248, groups:  Participant, 26

Fixed effects:
Estimate Std. Error t value
(Intercept)            1.89503    0.18013  10.521
WordOrderDO-S-IO       0.43269    0.12816   3.376
WordOrderIO-DO-S      -0.06250    0.12816  -0.488
WordOrderIO-S-DO       0.67308    0.12816   5.252
WordOrderS-DO-IO       2.12019    0.12816  16.543
WordOrderS-IO-DO       2.46635    0.12816  19.244
PhraseTypeNP-NP-Pro   -0.55128    0.14799  -3.725
PhraseTypeNP-Pro-NP    0.50000    0.14799   3.379
PhraseTypeNP-Pro-Pro  -0.03846    0.14799  -0.260
PhraseTypePro-NP-NP    0.98077    0.14799   6.627
PhraseTypePro-NP-Pro  -0.30769    0.14799  -2.079
PhraseTypePro-Pro-NP   0.42949    0.14799   2.902
PhraseTypePro-Pro-Pro -0.40385    0.14799  -2.729

Correlation matrix not shown by default, as p = 13 > 12.
Use print(x, correlation=TRUE)  or
vcov(x)     if you need it


Judging just from the intercepts, this looks as I would have expected it; e.g. having the subject at the beginning has a positive effect. (If I interpret the results correctly). What I do not really understand is why I do not get those significance stars (“**” etc.).

Is this the right approach, or should I use ordinal logistic regression? (I tried that as well but didn't get anything working).

Sorry for the longish post.

• If I were you I'd prefer to treat the rating as continuous (not ordinal) and use the lme function from the nlme package. Sep 16, 2016 at 9:51
• re: the lack of "significance stars" - when you use glmer() and don't specify a family (and link function) it effectively assumes you just want linear regression and so uses lmer() instead. lmer() doesn't give p values Sep 16, 2016 at 9:53
• @mark999 For the types of fully crossed designs common in (psycho) linguistics, lmer from lme4 is much better than lme. Most of the published literature in (pyscho)linguistics that uses mixed models uses lme4, including the more tutorial-type publications. Sep 18, 2016 at 1:17
• @Livius Thanks for the comment. I tend to use lme4 only if nlme doesn't have what I want (e.g. generalised linear mixed models). Nice answer. Sep 18, 2016 at 2:37

There are several aspects to the problem, some of which you addressed explicitly in your question, some of which you didn't.

# What type of scale are ratings?

Following Stevens (1946), ratings are ordinal and thus not well suited to linear regression (or its mixed-effects extensions). However, as a rule of thumb, the regression methods that can work on ordinal scales (ordinal logistic, arguably some forms of quantile regression) are more "exotic" and thus less well supported by the various software packages. When you combine "exotic" regression with "mixed-effects", the support is even worse. If you're willing to go the Bayesian route, then the wonderful MCMCglmm package does support mixed-effects ordinal logistic regression and will even output something like to p-values based on its MCMC estimates.

That said, many people in your position would just treat ratings as an interval scale and thus suitable for linear regression. It is up to you to determine which interpretation is appropriate and to make a convincing case for your interpretation to potential reviewers and other consumers of your scientific output. In the following, I will more or less assume that you went the linear regression route, although much of what I say is still applicable for ordinal logistic and other regression techniques.

# lme4: Where are my p-values?

This is an R FAQ. A short summary would be

degrees of freedom aren't as trivial as you were taught when dealing with models where the different parts can partially or fully overlap (i.e. hierarchical/multilevel models with partially or fully crossed levels).

In traditional repeated-measures ANOVA (which is actually a special case of the linear mixed-effects model with a particular symmetry assumption (sphericity)) and other simpler models, it is possible to calculate the degrees of freedom analytically (i.e. there is a formula for them). However, for more general designs, this is not possible and we must approximate the degrees of freedom. There are several different ways to do this (Satterthwaite, Kenward-Roger, etc.) but all are approximations and can at times be quite computationally intensive. Without degrees of freedom, it is not possible to convert t-values to p-values, nor F-values to p-values for the ANOVA case. (After all, the denominator degrees of freedom for the F-test are the square of the degrees of freedom for the t-test.) Given all the issues with p-values and their interpretation (finding the relevant literature is left as an exercise for the reader, but the psychological science literature has taken to self-flagellation trying to come to terms with it, so it's a good place to start), why would you bother with a computationally expensive approximation?

If you absolutely must have p-values, there are several ways to compute them beyond these approximations, including likelihood-ratio tests, although those may be less interesting for your particular model.

You could also recall that the t-distribution approaches the normal distribution as its degrees of freedom approach infinity. This means that for sufficiently large degrees of freedom (which you often have for (psycho)linguistic data), you can treat t-values as z-values. The traditional 5% significance level occurs for (two-sided, which you have/want) at $z = 1.96$. To be a little bit more conservative, you can just say that you reach significance at $t = 2$ (cf. Baayen, Davidson and Bates, 2008, amongst others for a similar train of thought).

But again, I would discourage you from getting too attached to statistical significance. Think about the size of differences and how well we can see that size amongst the noise/error, and not about whether there is a difference!

# How do I interpret my regression coefficients?

There are lots of good explanations of this online and in various textbooks, so I won't go into much detail here, but I think your fundamental problem is in your choice of contrasts. A linguistically motivated discussion of contrast coding (with links to other good discussions) can be found here, but you might want to look at a textbook discussion of contrast coding first (e.g. in one of Andy Field's Discovering Statistics books, which rely on slightly less mathematical intuition than other explanations I've seen). If you're looking at coefficient directly, then you should probably be using sum coding, and I would recommend using contr.Sum from the car package as it gives better names to the coefficients than the built-in contr.sum (note the capitalization).

Finally, it looks like you may actually be more interested in certain pairwise contrasts than anything else. In that case, you may want to consider the lsmeans package, as this is the most convenient way to do post-hoc pairwise comparison. For your case, this has advantages beyond being easier to interpret than regression contrasts/coefficients:

1. It's not sensitive to the choice of contrast coding (i.e. treatment vs. sum coding, obviously if you structure your predictors completely differently that will make a difference)
2. It will give you p-values for the pairwise estimates along with a multiple comparisons correction for the doing so many tests.

# On the structure of you model

Your model is a good start, but based on my understanding of your experimental question, it still needs some improvement:

1. You have no interaction term: Surely you care about how word order and phrase type interact with each other?
2. Depending on what type of contrasts you're interested in, it may make sense to have several categorical predictors for word-order and phrase-type, e.g. three predictors first.arg, second.arg, third.arg (with possibilities S, DO, IO for each of them) or maybe subject, direct, indirect (with possibilities first, second, third) for each of them. Depending on the exact contrasts you set up, some things will be easier to read from the coefficients and some things will be harder. You may also get a warning about dropping columns due to a rank-deficient matrix -- this is not actually a real problem, but just means the model doesn't compute coefficients for combinations which don't occur (e.g. both subject and direct simultaneously have the value first).
3. You need a random-effects term for Item: this isn't like classic ANOVA-by-subjects vs. ANOVA-by-items ($F_1$, $F_2$, cf. Clark 1973), you can and should model both sources of variation (subjects/participants and items) simultaneously, otherwise you may come to the wrong conclusions (cf. Judd, Westfall and Kenny 2012).
4. You may want a random-effects structure that includes slopes and intercepts, but there is a lot of debate on this. See e.g. Barr, Levy and Sheepers 2013; Barr 2013; Bates, Kliegel, Vasishth and Baayen 2015/arXiv as well as the GLMM FAQ linked above.

# One last thing to get started ...

A final bit of advice: it looks like you're completely new to mixed models, maybe even completely new to statistics. That's fine, we all had to start somewhere. But you might want to do some more background reading. :) Although slightly out-of-date on the software front, there is a 2008 special issue of the Journal of Memory and Language on mixed models. Other popular (book) references for this are Gelman & Hill (2006), Fox (2016 -- newest edition has much more on mixed models, but there are also several good older editions) and Pinheiro & Bates (2000). Finally, Bodo Winter has an online tutorial that many of my students have also found quite good.

• Thanks a lot, @Livius! Indeed, I am a newbie; but this looks really helpful! I will work through it. Sep 18, 2016 at 4:57