Walkthrough of an "empirical logit" analysis in R

For this walkthrough, I am using R version 2.12.0 and lme4 package version 0.999375.36. If you are using different versions, your results might differ (hopefully only slightly, unless you're using a really old version, in which case I recommend that you upgrade). For this same reason, they will also slightly differ from the results I presented in the MLR paper. Here's how to find out what versions you are using.

R.Version()$version.string
packageVersion("lme4")

[1] "R version 2.12.0 (2010-10-15)"
[1] ‘0.999375.36’

First, we need to load the lme4 package and then load in the data from the experiment. We have two files (bySubj.RData and byItem.RData) for the two by-subject and by-item analyses that we will perform.

library(lme4)

load("bySubj.RData")
load("byItem.RData")

ls()

[1] "byItem"      "bySubj"

The data frames bySubj and byItem contain the data aggregated by Subject and by Item respectively. The portion of the data we are looking at contains twelve trials per subject, three in each of the four cells of the design. The data was prepared by aggregating over 50 ms time slices (3 frames of eye data) and over all of the trials in a given condition within each subject (or item, for the analysis by item). Thus, each 50 ms bin consists of 12 frames of eye data.

For example, take a look at some of the data for subject 30.

subset(bySubj, SubjID==30)

   SubjID AggID SpeakerT CogLoadT Bin   t1 y  N
1      30   241     Same     None   0 0.00 0 12
2      30   241     Same     None   1 0.05 0 12
3      30   241     Same     None   2 0.10 0 12
4      30   241     Same     None   3 0.15 0 12
5      30   243     Diff     None   0 0.00 3 12
6      30   243     Diff     None   1 0.05 3 12
7      30   243     Diff     None   2 0.10 5 12
8      30   243     Diff     None   3 0.15 4 12
9      30   245     Same     Load   0 0.00 3 12
10     30   245     Same     Load   1 0.05 3 12
11     30   245     Same     Load   2 0.10 7 12
12     30   245     Same     Load   3 0.15 6 12
13     30   247     Diff     Load   0 0.00 3 12
14     30   247     Diff     Load   1 0.05 3 12
15     30   247     Diff     Load   2 0.10 3 12
16     30   247     Diff     Load   3 0.15 3 12

Note that we have aggregated together all of the trials within a given condition for each Subject, and broken up the time variable into a series of 50 ms bins. The bins for a given aggregated condition (henceforth "aggregate") were assigned an ID number (to represent that they are all from the same condition for the same subject). The precise number that is given to a particular "aggregate" does not matter; what matters is that the number given to all the bins in a given aggregate is unique. This is a way of telling the model fitting function lmer that these bins belong together.

The variable y is the number of frames in the 50 ms window for which the gaze fell within the bounds of the target object; The variable N is the total number of frames within each bin.

Variable t1 is the time variable, coded in seconds, from the onset of the analysis window (which was 300 ms after the start of the word).

We can compute the empirical logit for each bin using the formula \( elog = log\left(\frac{Y + .5}{N - Y + .5}\right) \). This is an approximation to log odds. Because we are using a linear approximation, we also will need to calculate weights \( \frac{1}{Y+.5} + \frac{1}{N-Y+.5} \) because the variance of the logit depends on the mean.

bySubj$elog <- log( (bySubj$y + .5) / (bySubj$N - bySubj$y + .5) )
bySubj$wts <- 1/(bySubj$y + .5) + 1/(bySubj$N - bySubj$y + .5)

We will now create condition codes for our categorical predictor variables. We will use ANOVA-style "deviation" coding.

bySubj$S <- ifelse(bySubj$SpeakerT=="Same",-.5,.5)
bySubj$L <- ifelse(bySubj$CogLoadT=="None",-.5,.5)
head(bySubj, 12)

   SubjID AggID SpeakerT CogLoadT Bin   t1 y  N       elog       wts    S    L
1      30   241     Same     None   0 0.00 0 12 -3.2188758 2.0800000 -0.5 -0.5
2      30   241     Same     None   1 0.05 0 12 -3.2188758 2.0800000 -0.5 -0.5
3      30   241     Same     None   2 0.10 0 12 -3.2188758 2.0800000 -0.5 -0.5
4      30   241     Same     None   3 0.15 0 12 -3.2188758 2.0800000 -0.5 -0.5
5      30   243     Diff     None   0 0.00 3 12 -0.9985288 0.3909774  0.5 -0.5
6      30   243     Diff     None   1 0.05 3 12 -0.9985288 0.3909774  0.5 -0.5
7      30   243     Diff     None   2 0.10 5 12 -0.3101549 0.3151515  0.5 -0.5
8      30   243     Diff     None   3 0.15 4 12 -0.6359888 0.3398693  0.5 -0.5
9      30   245     Same     Load   0 0.00 3 12 -0.9985288 0.3909774 -0.5  0.5
10     30   245     Same     Load   1 0.05 3 12 -0.9985288 0.3909774 -0.5  0.5
11     30   245     Same     Load   2 0.10 7 12  0.3101549 0.3151515 -0.5  0.5
12     30   245     Same     Load   3 0.15 6 12  0.0000000 0.3076923 -0.5  0.5

The model that we are fitting is a simple line; the regression analysis will estimate the slope and intercept for the line, as well as how the various factors (and their interaction) modulate the slope and intercept terms (see the MLR paper for details). So the next step is to fit the full model, including all random effects for the Subject analysis.

bySubj.lmer <- lmer(elog ~ t1*S*L +  
                    (1 + t1 | AggID:SubjID) + (1 + t1 | SubjID), 
                    data=bySubj)

Note that we have broken what is actually a single command across three lines. This makes it easier to read (and explain).

First, bySubj.lmer is the name of the variable that will be created by the model fitting procedure lmer. The first argument to the lmer function is a formula, specifically elog ~ t1*S*L + (1 + t1 | SubjID:AggID) + (1 + t1 | SubjID). Note the tilde (~) where an equal sign would normally be; this is the conventional way to write a formula in R. The terms on the first line represent the DV (elog) and the predictor variables, which determine what fixed effects will be estimated. We used the syntax t1*S*L as a shorthand to get all main effects and interactions between time (t1), speaker (S), and load (L), rather than laboriously writing out t1+S+L+S:L+t1:S+t1:L+t1:S:L.

The second line (1 + t1 | SubjID:AggID) + (1 + t1 | SubjID) represents the random effects in the model. The pipe symbol | divides what we want to allow to vary (slope and intercept) from the units over which we will allow them to vary (aggregate nested within subject, SubjID:AggID, and subject, SubjID).

The third line just tells lmer the name of the data frame where our variables are located.

Now let's have a look at the results.

 1:  Linear mixed model fit by REML 
 2:  Formula: elog ~ t1 * S * L + (1 + t1 | AggID:SubjID) + (1 + t1 | SubjID) 
 3:     Data: bySubj 
 4:    AIC  BIC logLik deviance REMLdev
 5:   2088 2160  -1029     2058    2058
 6:  Random effects:
 7:   Groups       Name        Variance  Std.Dev. Corr
 8:   AggID:SubjID (Intercept)  1.601033 1.26532         
 9:                t1          55.756646 7.46704  -0.451 
10:   SubjID       (Intercept)  0.021771 0.14755         
11:                t1           0.399988 0.63245  1.000  
12:   Residual                  0.160052 0.40007         
13:  Number of obs: 896, groups: AggID:SubjID, 224; SubjID, 56
14:  
15:  Fixed effects:
16:              Estimate Std. Error t value
17:  (Intercept) -1.37923    0.08965 -15.385
18:  t1           2.63059    0.55965   4.700
19:  S            0.31955    0.17490   1.827
20:  L           -0.02207    0.17490  -0.126
21:  t1:S        -2.21968    1.10648  -2.006
22:  t1:L        -0.30371    1.10648  -0.274
23:  S:L         -0.20553    0.34980  -0.588
24:  t1:S:L       1.31463    2.21296   0.594

There is a lot of complicated output here, so let's take it step by step.

Line 5 tells us about the quality of the model fit, including the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), the log likelihood at convergence (logLik), the deviance (which is just -2 times the log likelihood), and something called the REMLdev which is the REML deviance from the restricted maximum-likelihood technique that lme4 has used. If we wanted to get values for the likelihood that we could use for likelihood tests, we could fit the model using the option REML=FALSE. However, we will not be doing this in the current example.

Lines 6-12 give us information about the random effects in the model. These values characterize the variance/covariance of two bivariate distributions, one at the subject level and one that describes variation associated with individual aggregates (cells of the design). There is also a residual term.

Finally, lines 15-24 show the fixed effects in the model. Note that we get t values, but no p values. If we want to get p-values we can treat the ts as if they were drawn from a normal distribution, using the pnorm function. The code below does this.

ts <- fixef(bySubj.lmer) / sqrt(diag(vcov(bySubj.lmer)))
print(round(ts,2))
ps <- 2*(1-pnorm(abs(ts)))
print(round(ps,3))

(Intercept)          t1           S           L        t1:S        t1:L 
     -15.39        4.70        1.83       -0.13       -2.01       -0.27 
        S:L      t1:S:L 
      -0.59        0.59
 (Intercept)          t1           S           L        t1:S        t1:L 
      0.000       0.000       0.068       0.900       0.045       0.784 
        S:L      t1:S:L 
      0.557       0.552

We calculate the t values by taking the parameter estimates and dividing them by their standard errors. We have to extract all of this information from the bySubj.lmer object that was created as the result of the call to the lmer function. We extracted the parameter estimates using the fixef function, and then divided through by the standard error, obtained using sqrt(diag(vcov(bySubj.lmer))). The function vcov(bySubj.lmer) gets the variance-covariance matrix for the fixed effects. Their variances are in the diagonal, which we extract using the diag function; and then we take the square root (sqrt) to get the standard errors. Another way to get p-values would be to use likelihood ratio tests, which we don't cover here.

The listing includes main effects and interactions. The effects S, L, and S:L are the effects of speaker, load, and their interaction at the intercept (where t1 equals 0), and thus represent anticipatory baseline effects. Note that in line 19 we see what looks like a trend toward a main effect of speaker. Because we coded the speaker variable as Same=-.5 and Diff=.5, a positive parameter estimate of .31955 implies that the chance of looking at the target at the intercept was higher when the speaker was different than when he/she was the same; we obtained a p-value of .068 for this effect. (See the MLR paper for a possible explanation for this baseline effect). What we also see in line 21 is a speaker:time interaction, with the slope steeper for the same-speaker condition (since "same" was coded as negative) by 2.22 logits, with a p-value of .045.

So, we have an anticipatory baseline effect of Speaker, as well as an effect of Speaker on the slope that goes in the opposite direction (high slopes for same speaker), with no other effect approaching significance. This supports the idea that listeners were facilitated by hearing the precedent spoken in the same voice.

One can replicate the same analysis by Item using the code below. This yields basically the same pattern of results.

byItem$elog <- log( (byItem$y + .5) / (byItem$N - byItem$y + .5) )
byItem$wts <- 1/(byItem$y + .5) + 1/(byItem$N - byItem$y + .5)
byItem$S <- ifelse(byItem$SpeakerT=="Same",-.5,.5)
byItem$L <- ifelse(byItem$CogLoadT=="None",-.5,.5)

byItem.lmer <- lmer(elog ~ t1*S*L +  
                    (1 + t1 | AggID:ItemID) + (1 + t1 | ItemID), 
                    data=byItem)

print(summary(byItem.lmer),corr=F)

ts <- fixef(byItem.lmer) / sqrt(diag(vcov(byItem.lmer)))
print(round(ts,2))
ps <- 2*(1-pnorm(abs(ts)))
print(round(ps,3))

 Linear mixed model fit by REML 
Formula: elog ~ t1 * S * L + (1 + t1 | AggID:ItemID) + (1 + t1 | ItemID) 
   Data: byItem 
  AIC  BIC logLik deviance REMLdev
 1087 1151 -528.5     1061    1057
Random effects:
 Groups       Name        Variance   Std.Dev. Corr   
 AggID:ItemID (Intercept)  1.5442939 1.242696        
              t1          42.4860732 6.518134 -0.505 
 ItemID       (Intercept)  0.0844539 0.290609        
              t1           0.0019619 0.044293 -1.000 
 Residual                  0.1228154 0.350450        
Number of obs: 512, groups: AggID:ItemID, 128; ItemID, 32

Fixed effects:
            Estimate Std. Error t value
(Intercept) -1.41348    0.12400 -11.399
t1           3.05843    0.63933   4.784
S            0.41460    0.22571   1.837
L           -0.04374    0.22571  -0.194
t1:S        -2.72025    1.27856  -2.128
t1:L        -0.25802    1.27856  -0.202
S:L         -0.18244    0.45142  -0.404
t1:S:L       0.66504    2.55713   0.260
 (Intercept)          t1           S           L        t1:S        t1:L 
     -11.40        4.78        1.84       -0.19       -2.13       -0.20 
        S:L      t1:S:L 
      -0.40        0.26
 (Intercept)          t1           S           L        t1:S        t1:L 
      0.000       0.000       0.066       0.846       0.033       0.840 
        S:L      t1:S:L 
      0.686       0.795