gmane.comp.lang.r.linguistics

Re: Simpler model with random slopes

Levy, Roger — 2012-03-31T00:38:23

Hi Sverre,

This is a belated follow-up to your follow-up questions.

On Mar 20, 2012, at 3:49 AM PDT, Sverre Stausland wrote:


Critical question: by "assess the contribution" do you mean assess whether there is a reliable fixed effect of Indep2 (i.e. whether Indep2 has a reliable effect in one direction or another that generalizes across subjects)?  Or assessing whether including Indep2 at all in your model helps its ability to predict Dep, even if Indep2 may matter only in a subset of subjects?


--

Roger Levy                      Email: rlevy-XkckGZ689+c< at >public.gmane.org
Assistant Professor             Phone: 858-534-7219
Department of Linguistics       Fax:   858-534-4789
UC San Diego                    Web:   http://idiom.ucsd.edu/~rlevy

Re: Conflicting p-values from pvals.fnc

Levy, Roger — 2012-03-30T19:34:32

On Mar 28, 2012, at 7:13 PM PDT, Levy, Roger wrote:


On Mar 27, 2012, at 3:19 PM PDT, Tom Gijssels wrote:

Hey Roger,

Thanks for the response!

Let me try to address your questions as good as possible (and in reverse order to make things a bit more concrete).

The reason for including weights is that there's a reasonable amount of variance in the relative importance of the data-points I'm using. To be a bit more specific: the experiment looks at whether people accommodate their speech rate to a conversational partner's speech rate. Speech rate was calculated by taking the words per second for each of the subject's utterances during Baseline and during Conversation. Since a subject's utterances vary in duration (e.g. some are 2s, others 10s), I want to make sure that the speech rate measure of a long utterance gets contributes more strongly to the model than a short utterance. Therefore, I wanted to include the duration for each utterance as a weight for the corresponding speech rate measure. Does this sou

Re: Conflicting p-values from pvals.fnc

Nathaniel Smith — 2012-03-29T13:45:01

I'm not sure what it means that two tests of unknown quality disagree
with a third test of unknown quality... but figuring out which tests
might be accurate is actually very easy. You already wrote code to
simulate data from your model; just run it 1000 times with no real
effects put in, and look at the distribution of p values you get. An
accurate test will produce p values that are uniformly distributed,
i.e., 0.05 * 1000 = ~50 of them should be less than 0.05.

(Also I notice your simulation code uses different weights for the
different conditions, but doesn't actually add more noise to one
condition than the other; you might want to fix that first.. You just
need something like err <- err * 1/sqrt(weights).)

Re: Conflicting p-values from pvals.fnc

Levy, Roger — 2012-03-29T02:13:14

On Mar 27, 2012, at 3:19 PM PDT, Tom Gijssels wrote:

Hey Roger,

Thanks for the response!

Let me try to address your questions as good as possible (and in reverse order to make things a bit more concrete).

The reason for including weights is that there's a reasonable amount of variance in the relative importance of the data-points I'm using. To be a bit more specific: the experiment looks at whether people accommodate their speech rate to a conversational partner's speech rate. Speech rate was calculated by taking the words per second for each of the subject's utterances during Baseline and during Conversation. Since a subject's utterances vary in duration (e.g. some are 2s, others 10s), I want to make sure that the speech rate measure of a long utterance gets contributes more strongly to the model than a short utterance. Therefore, I wanted to include the duration for each utterance as a weight for the corresponding speech rate measure. Does this sound sensible?

As far as I can see, it seems that this

Re: Conflicting p-values from pvals.fnc

Tom Gijssels — 2012-03-27T22:19:42

Hey Roger,

Thanks for the response!

Let me try to address your questions as good as possible (and in reverse
order to make things a bit more concrete).

The reason for including weights is that there's a reasonable amount of
variance in the relative importance of the data-points I'm using. To be a
bit more specific: the experiment looks at whether people accommodate their
speech rate to a conversational partner's speech rate. Speech rate was
calculated by taking the words per second for each of the subject's
utterances during Baseline and during Conversation. Since a subject's
utterances vary in duration (e.g. some are 2s, others 10s), I want to make
sure that the speech rate measure of a long utterance gets contributes more
strongly to the model than a short utterance. Therefore, I wanted to
include the duration for each utterance as a weight for the corresponding
speech rate measure. Does this sound sensible?

As far as I can see, it seems that this is at least conceptually what the
weights argument does

Re: Conflicting p-values from pvals.fnc

Levy, Roger — 2012-03-26T21:56:07

Hi Tom,

I'm not sure how lmer uses the weights argument -- I think I know how it *should* be used (setting a weight of k to a given observation should be equivalent to having seen k replicates of that observation), but what you've got here makes me wonder what's actually going on under the hood.  Could you maybe point us to some documentation as to what lmer claims to be doing with the weights argument, and maybe also explain why you want to use the weights you're using?

Best

Roger

On Mar 26, 2012, at 10:15 AM PDT, Tom Gijssels wrote:


Dear R-langers,

I'm trying to run a mixed effect model using the lmer() function and have
run into some issues in interpreting the p-values generated by
pvals.fnc(). The design is a between-subjects design, with two fixed
effects (condition & block; each with two levels), and one random effect
(subject). Additionally, I have a set of weights that I want to include.

When looking at the pvals.fnc() output,there appears to be a large
discrepancy between the pMCMC values an

Conflicting p-values from pvals.fnc

Tom Gijssels — 2012-03-26T17:15:46

Dear R-langers,

I'm trying to run a mixed effect model using the lmer() function and have
run into some issues in interpreting the p-values generated by
pvals.fnc(). The design is a between-subjects design, with two fixed
effects (condition & block; each with two levels), and one random effect
(subject). Additionally, I have a set of weights that I want to include.

When looking at the pvals.fnc() output,there appears to be a large
discrepancy between the pMCMC values and the t-statistic p-values. Whereas
one of the main effects and the interaction are far from significant
judging by the pMCMC values, they are highly significant when looking at
the t-statistic p-values (e.g. Condition: pMCMC = 0.2294; Pr(>|t|) = 0.0000
& Condition*Block: pMCMC = 0.3296; Pr(>|t|) = 0.0000) . I have read that
the t-statistic based p-values are less conservative, but the difference
between these two values seems really extreme.

Below some code that simulates the model and the data. The original data
set has two precise charac

Re: Simpler model with random slopes

Sverre Stausland — 2012-03-21T12:58:10

An update: I asked below whether one should treat a model with a main
effect Indep2 and random slopes for Indep2 similarly as one should
treat a model with a main effect Indep2 and an interaction between
Indep2 and another main effect. That is, a model with random slopes
for Indep2 should also contain the main effect Indep2 - one should not
remove the main effect Indep2 and leave the random slopes for Indep2
in.

For what it's worth, I just found the same question asked in R-sig-ME
in 2010: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003876.html,
to which Ben Bolker answers in the affirmative (with an explanation).

Sverre

On Tue, Mar 20, 2012 at 11:49 AM, Sverre Stausland
<johnsen-fWAZDB8bsKe+fmr0zi+kZQ< at >public.gmane.org> wrote:

Re: Simpler model with random slopes

Sverre Stausland — 2012-03-20T10:49:41

Thanks for pointing that out, Florian. I noticed just that when I
built my models with and without the covariance, and it makes sense
from Nathaniel's exemplary explanation.

I do have a follow-up question that I'm taking the liberty to ask in
this thread:

As an example, say my model includes random slopes for items with
respect to a variable Indep2, as in (a)

(a) glmer(formula = Dep ~ 1 + (1 | Subject) + (0 + Indep2 | Item) +
Indep1 + Indep2, data = my.data, family = binomial(link = "logit"))

To assess the contribution of variable Indep2 through model
comparison, should I then remove the main effect Indep2 as well as the
random slopes for items, as in (b) below?

(b) glmer(formula = Dep ~ 1 + (1 | Subject) + Indep1, data = my.data,
family = binomial(link = "logit"))

It seems slightly odd to me to allow the model to include random
variation among items with respect to a variable Indep2 if the model
doesn't assume that Indep2 plays a role as a main effect as well,
similar to how higher-order interactions

Re: Simpler model with random slopes

T. Florian Jaeger — 2012-03-20T01:19:31

just to add to Nathaniel's beautifully clear explanation: going from (a) to
(b) is really the opposite of "reducing" the model. going from (b) to (a)
would remove one parameter (the covariance parameter)

flo

On Mon, Mar 19, 2012 at 2:45 PM, Sverre Stausland
<johnsen-fWAZDB8bsKe+fmr0zi+kZQ< at >public.gmane.org>wrote:

Re: Simpler model with random slopes

Sverre Stausland — 2012-03-19T18:45:38

Thanks for the explanation. I think I got it!

On Mon, Mar 19, 2012 at 5:54 PM, Nathaniel Smith <njs-e+AXbWqSrlAAvxtiuMwx3w< at >public.gmane.org> wrote:

Re: Simpler model with random slopes

Nathaniel Smith — 2012-03-19T16:54:25

On Mon, Mar 19, 2012 at 4:37 PM, Sverre Stausland
<johnsen-fWAZDB8bsKe+fmr0zi+kZQ< at >public.gmane.org> wrote:

(1 | Item) means that items are assumed to have different intercepts

Simpler model with random slopes

Sverre Stausland — 2012-03-19T16:37:15

Hi all,

this question is ultimately based on Florian's lecture1 slides here:
http://hlplab.wordpress.com/2010/05/10/mini-womm/

I'm doing a mixed model logistic regression, with random intercepts
for items and random slopes for items with respect to the fixed effect
Indep2 (cf. slide 85):

(a) glmer(formula = Dep ~ 1 + (1 | Item) + (0 + Indep2 | Item) +
Indep1 + Indep2, data = my.data, family = binomial(link = "logit"))

As per slide 88, I can also reduce the random effects to (1 + Indep2 | Item):

(b) glmer(formula = Dep ~ 1 + (1 + Indep2 | Item) + Indep1 + Indep2,
data = my.data, family = binomial(link = "logit"))

It's not exactly clear to me what (1 + Indep2 | Item) does, since the
output of both (a) and (b) includes random intercepts for items and
random slopes for items by Indep2. At the same time, model (a) and (b)
differ in their exact estimates.

I would appreciate if someone could explain what the difference
between model (a) and (b) is.

Thanks
Sverre

Re: Questions about reporting mixed-effects results

T. Florian Jaeger — 2012-02-23T20:19:35

Hi Ariel,

I hope that we're moving towards researchers giving full model summaries in
electronic appendices for those of us interested in the details, while
keeping the main text focused on a concise (but interpretable) summary of
the results. The appendices, in my view, should contain detailed
information about how variables were coded and entered into the regression.
Since many researchers in psycholinguistics and its neighboring disciplines
had their statistical training focused on ANOVA, I think using coding that
makes regressions maximally comparable to ANOVA is preferable. This means
that I would recommend to contrast code (or where justified Helmert or
polynomial code) factors rather than to employ the basic treatment coding
(R's default). I think centering of treatment coded factors, while it
doesn't address all concerns, usually doesn't result in reduced
interpretability either although for pretty balanced data sets you might as
well use contrast coding (thereby perhaps removing an issue that makes

Re: Questions about reporting mixed-effects results

2012-02-20T21:27:41

With respect to power, once you have settled on a model and have the  
parameter estimates in hand, you can always run a simulation of the  
model a jillion times using those paraemter estimates and see what  
proportion of the jillion yielded significant results on each of the  
tests you care about.  Those proportions will be the estimated power  
of those tests.

Quoting Jason Kahn <jmkahn-vEoVa6Jh+/J4piUD7e9S/g< at >public.gmane.org>:




Richard S. Bogartz
Professor of Psychology
UMASS, Amherst 01003

Afghanistan:
100 Al Quaeda
100,000 U.S. troops
$139,200,000,000 a year

OUT NOW!!

"When I was 5 years old, my mother always told me? that happiness was  
the key to life. When I went to school, they asked me what I wanted to  
be when I grew up. I wrote down ?happy?. They told me I didn?t  
understand the assignment, and I told them they didn?t understand life."

? John Lennon

Re: Questions about reporting mixed-effects results

Goldberg, Ariel M — 2012-02-20T20:59:22

Dear Jason,

Thanks so much for your advice.  I definitely will check out your website.

I should mention that I made a mistake in my previous message… all of my models contain random intercepts, but a few contain random slopes for the variables of interest.  These random slopes were added when model comparisons indicated that they were warranted.

In the cases where the model contained random slopes, I have reported the t-values as well as the results of a model comparison looking at whether the addition of the one variable of interest significantly increases the fit of the model.

Thanks again!
AG


On Feb 20, 2012, at 3:55 PM, Jason Kahn wrote:

 associated with the parameter estimates, assuming a reasonably large (at the very least 500 observations, preferably 1000) dataset. A slightly more conservative estimate is on my website, and comes courtesy of Roger Levy (if I remember right), who cites Hox (2010).

Re: Questions about reporting mixed-effects results

Jason Kahn — 2012-02-20T20:55:03

Hi Ariel,

Having been through some of these same issues with reviewers recently, I
hope I can say something helpful. Most of this is summarized on my website,
which I would appreciate feedback/commentary on from anyone and everyone,
and which you can find here: http://jmkahn.web.unc.edu/modeling/

You should definitely report centered predictors if their theoretical
interpretation or collinearity-with-interactions depends on that
transformation. Centering will reduce (or, ideally, eliminate) the
collinearity between "first-order" predictors and the interactions they
participate in, unless you have unbalanced conditions. It will also change
the interpretation of your intercept, and the parameter estimates
associated with those variables. If you and/or your readers care about any
of those things, reporting the centered predictors seems necessary.
Otherwise, I don't believe it matters.

Things may have changed since I last used pvals.fnc(), but it used to work
for models with random intercepts, and even models

Questions about reporting mixed-effects results

Goldberg, Ariel M — 2012-02-20T20:15:19

Dear all,

I have a few questions about how to report the results of mixed-effects analyses for publication. I have been perusing the Jaeger & Kuperman presentation but a few questions remain.  

I have been asked by the reviewers to include a full regression table, which I take to comprise coefficient estimates, MCMC-based confidence intervals and MCMC-based p-value estimations.
-Should the model that I use to report these values contain uncentered predictors, centered predictors, or centered and scaled predictors?
-A few of my models involve random intercepts, and I believe that pvals.fnc() is not currently defined for models with random intercepts.  Do you have any suggestions for how I should report these models? 

My models contain many control variables and only one or two variables that I am actually concerned with.  As such, I have not worried about multicollinearity among the control variables.  I suppose I should just state this somewhere to facilitate the interpretation of the regression tables?

Re: negative deviances (again)

T. Florian Jaeger — 2012-01-28T19:14:37

Argh, sorry, I was too trigger happy. I just adapted your example to mixed
models:

library(stats)
library(lme4)
N<-100
z <- rep(1:5, 20)
x <- runif(N)
y1 <- rnorm(N,mean=x,sd=3) + rnorm(2,0,sd=)[z]
y2 <- rnorm(N,mean=x,sd=.3) + rnorm(2,0,sd=)[z]
logLik(lmer(y1 ~ x + (1|z)))
logLik(lmer(y2 ~ x + (1|z)))

and get the expected behavior (positive LLs for low residual SDs and
negative ones for larger SDs). Seems like I simply got unlucky in choosing
the test data sets I ran before (I should have run the sim you suggested
immediately). In retrospect, perhaps this is due that all the test data
contained log-transformed DVs ;). If one unlogs Baayen's RT measure in the
example I sent, the residuals are, of course, large, and then the deviance
is positive. phew.

So, I guess, higher LLs still are better fits since positive deviance are
in principle legit. Sorry, for pressing the panic button.

Florian

On Sat, Jan 28, 2012 at 1:48 PM, piantado <piantado-3s7WtUTddSA< at >public.gmane.org> wrote:

Re: negative deviances (again)

T. Florian Jaeger — 2012-01-28T18:51:47

Hi Steve,

I know that this can happen. I am worried because it seems to occur in ALL
models I've tried on different data sets.

I also meant to say that I am even more confused (troubled) about the
following. As reminder, here are the two models - notice that the negative
Deviance (deviance at ML max) goes down (i.e. becomes larger) as the
additional predictor is added to the model.

Linear mixed model fit by REML
Formula: RT ~ Frequency + (1 | Subject)
   Data: lexdec
    AIC    BIC logLik deviance REMLdev
 -858.4 -836.8  433.2   -880.9  -866.4
[snip]

Linear mixed model fit by REML
Formula: RT ~ Frequency + Trial + (1 | Subject)
   Data: lexdec
    AIC  BIC logLik deviance REMLdev
 -846.1 -819  428.1   -887.2  -856.1
[snip]

When I run an ANOVA over the two models in the previous email (and I've
replicated this for other data), it'll use the log-Likelihood estimate from
the deviance of the ML-fitted (rather than REML-fitted) model, so that the
more complex model has a larger (here: more positive; if the s

Re: negative deviances (again)

piantado — 2012-01-28T18:48:27

(Sorry my email crashed and I'm not sure this went through the first
time)


Hi Florian, 

What makes you sure there's a problem? You can get positive log
likelihoods (and negative deviances) from the residual distribution
having a low enough standard deviation, since the normal distribution
goes above the line y=1. 

It looks like that's what's up here:

l <- lmer(RT ~ Frequency + (1 | Subject), lexdec)
[1] 0.1798177

plot(dnorm(seq(-3,3,0.01), sd=0.18), type="l")


Here's an example where you can twiddle the SD and see the log
likelihood change sign:

library(stats)
N<-100
x <- runif(N)
y <- rnorm(N,mean=x,sd=0.18)
logLik(lm(y ~ x))


++Steve