11
$\begingroup$

The bounty I placed on this question expires in the next 24 hours.

I have a psychological data set which, traditionally, would be analysed using a paired samples t test. The design of the experiment is $39 (subjects) \times 7 (targets) \times 2 (conditions)$, and I'm interested in the difference in a given variable between the conditions.

The traditional approach has been to average across targets so that I have 2 observations per participant, and then compare these averages using a paired t test.

I wanted to use a mixed models approach, as has become increasingly popular in this field (i.e. Baayen, Davidson & Bates, 2008), and so the first model I fit, which I thought should approximate the results of the t test, was one with $condition$ as a fixed effect, and random intercepts for $subjects$ (i.e. $var = \alpha + \beta*condition + Intercept(subject) + \epsilon$. Obviously, the full model would also include random intercepts for $targets$.

However, I'm struggling to understand why I achieve pretty divergent results between the two approaches. Can anyone explain what's going on here? I've also seen (what I understand to be) a similar question asked here, with an answer about correlation structure which I'm not equipped to understand. If this is also what's at issue here, I would appreciate if anyone could suggest some resources to read up on this.

Edit: I've posted the example data, and R script, here.

Edit #2 - Bounty added

Some additional points:

  • I'm only analysing the correct responses (think of it as analogous to reaction time), so there are missing cases - not every participant provides 7 data points per condition.
    • When I analyse all responsees, rather than just the correct ones, the difference between the two results is reduced, but not eliminated. This suggests to me that the missing cases are a factor here.
  • The variable isn't normally distributed. In my final model, I scale it using a Box-Cox transformation, but I omit that here for consistency with the t test.
  • As pointed out by @PeterFlom, the $df$s differ hugely between the two approaches, but I assume this to be because the t test is being applied to the aggregate data (2 observations per participant, 1 per condition), while the mixed model is applied to raw scores ($<14$ observations per participant, $<7$ per condition).
  • @BenBolker notes that the t values also differ pretty considerably.

My analysis code is below.

>library(dplyr)
>subject_means = group_by(data, subject, condition) %>% summarise(var=mean(var))
>t.test(var ~ condition, data=subject_means, paired=T)

    Paired t-test

data:  var by condition
t = -1.3394, df = 37, p-value = 0.1886
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.14596388  0.02978745
sample estimates:
mean of the differences 
            -0.05808822 

>library(lme4)
>lm.0 = lmer(var ~ (1|subject), data=data)
>lm.1 = lmer(var ~ condition + (1|subject), data=data)
>anova(lm.0, lm.1)

Data: data
Models:
object: var ~ (1 | subject)
..1: var ~ condition + (1 | subject)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
object  3 489.09 501.23 -241.55   483.09                           
..1     4 485.81 502.00 -238.90   477.81 5.2859      1     0.0215 *

>library(lmerTest)
>summary(lm.1)$coef

              Estimate Std. Error        df  t value     Pr(>|t|)
(Intercept) 0.11862462 0.02878027  98.60659 4.121734 7.842075e-05
condition   0.09580546 0.04161237 400.27441 2.302331 2.182890e-02

Notice, specifically, the jump in the p value from $p = .188$ in the t test, to $p = .021$ from either lmer method.


I've tried, and failed to provide a reproducible example of this, using the anorexia dataset in the MASS package, so I would assume the problem is something idiosyncratic to my data, but I don't understand what.

# Borrowing from http://ww2.coastal.edu/kingw/statistics/R-tutorials/dependent-t.html
>data(anorexia, package="MASS")
>ft = subset(anorexia, subset=(Treat=="FT"))
>wgt = c(ft$Prewt, ft$Postwt)
>pre.post = rep(c("pre","post"),c(17,17))
>subject = rep(LETTERS[1:17],2)
>t.test(wgt~pre.post, data=ft.new, paired=T)

    Paired t-test

data:  wgt by pre.post
t = 4.1849, df = 16, p-value = 0.0007003
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  3.58470 10.94471
sample estimates:
mean of the differences 
               7.264706 

>m = lmer(wgt ~ pre.post + (1|subject), data=ft.new)
>summary(m)$coef

             Estimate Std. Error       df   t value     Pr(>|t|)
(Intercept) 90.494118   1.689013 26.17129 53.578096 0.0000000000
pre.postpre -7.264706   1.735930 15.99968 -4.184908 0.0007002806
$\endgroup$
5
  • 1
    $\begingroup$ Something strange is going on. Look at the df in the mixed model vs. the t-test. Condition has df = 400. In the t it is 37. That is surely odd. $\endgroup$
    – Peter Flom
    Commented Jul 5, 2014 at 22:49
  • $\begingroup$ ... and it's not just the df. The paired t-test gives you a t-statistic of -1.3394 , while lmer gives you a t-statistic of 2.302331 . Can you post data somewhere? $\endgroup$
    – Ben Bolker
    Commented Jul 6, 2014 at 6:34
  • $\begingroup$ Thanks for the comments. I've made a major edit to the question, addressing both statistical points. To summarise: a) I assume the change in df is down to the use of aggregate data in the t test, and raw data in the mixed model; b) I also presume the crux of the matter is the difference in t values, and suspect it might be to do with the missing cases; and c) I'll post the data in the next few minutes. $\endgroup$
    – Eoin
    Commented Jul 6, 2014 at 17:08
  • $\begingroup$ To continue - d) the discrepancy between the two results is less when I analyse all trials, not just the correct responses, if this casts any light. Finally, I've posted the essentials of my code and R script here. Thanks again to everyone taking an interest here. $\endgroup$
    – Eoin
    Commented Jul 6, 2014 at 17:42
  • $\begingroup$ Just wondering, have you ever figured out what was wrong here? $\endgroup$
    – fsociety
    Commented Oct 6, 2015 at 0:56

1 Answer 1

4
+50
$\begingroup$

I think the problem is the way the paired t test is computed. Try this:

t.test(all_by_sub$var[all_by_sub$condition==1], all_by_sub$var[all_by_sub$condition==0], paired=TRUE)

This gives:

data:  all_by_sub$var[all_by_sub$condition == 1] and all_by_sub$var[all_by_sub$condition == 0]
t = 2.0529, df = 37, p-value = 0.0472
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.0009400428 0.1435297005
sample estimates:
mean of the differences 
             0.07223487 

Obviously, t.test computes x minus y in the paired t test. That's why the sign of the estimate and the t value was reversed. Beyond that, both the estimate (0.072 vs. 0.082) and the t value (2.05 vs. 2.19) of the mixed model are very close to the results of the t test:

            Estimate Std. Error      df t value Pr(>|t|)
(Intercept)    0.118      0.028 106.772   4.159    0.000
condition      0.082      0.037 462.992   2.192    0.029
$\endgroup$
1
  • 1
    $\begingroup$ You're right on all this, but this is analysing all of the data points (which I do in the first part of the script posted, and note that they only diverge a little). It's when I exclude the incorrect responses that the two approaches diverge, as in the second part of the script, the results of which are in the original question. Thanks though! $\endgroup$
    – Eoin
    Commented Jul 8, 2014 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.