[R-sig-ME] Random effects in multinomial regression in R?
René
b|mono@om @end|ng |rom gm@||@com
Sun Mar 24 14:39:26 CET 2019
Ok. I try a short one :)
Also going back to your first mail:
This model is what you want:
score~Time*Group+(1+Time|subject)+(1+Time*Group | words)
No skew in here ;) (but maybe model identification issues depending on the
number of observations you got.) In other words, (random) by-word
variations and (random) by-participant variations are taken into account,
and the "average across stimuli" problem articulated here (Judd, Westfall,
& Kenny, 2012) is appropriately tackled. As you do not explicitly reference
to this paper, I just assume that you mean this kind of "skew" (alpha
inflation due to stimulus aggregation), because I am currently not aware of
other theoretically or statistically relevant 'skews' there could be.
Best, René
Ps: Why not using brms... :) You get parameter likelihood estimates and
Bayesian credible intervals for your effect estimates? The model code and
output-format are basically the same. And you don't need to report p-values
:))
Pps: For making the above formula work, you need to recode your data matrix
to:
Subject word Group Time Score
1 1 control Pre 1
1 1 control Post 1
1 2 control Pre 1
1 2 control Post 1
1 3 control Pre 0
1 3 control Post 0
1 4 control Pre 0
1 4 control Post 0
Am So., 24. März 2019 um 09:36 Uhr schrieb Souheyla GHEBGHOUB <
souheyla.ghebghoub using gmail.com>:
> Dear René,
>
> When I do Time (PrePost), score (response) will return to be dichotomous
> (0, 1) and I will be using logistic regression (glmer) using lme4 package
> instead of brms.
>
> So my previous question was about glmer model and whether its results
> reflect a prepost comparison based on word to word analysis of each
> participant instead of their sums; the latter will skew results.
> I appreciate your constant help,
> Thank you
> Souheyla
>
> On Sat, 23 Mar 2019, 22:01 René, <bimonosom using gmail.com> wrote:
>
>> Souheyla, it's a solvable problem :)
>> But does the model analyse at the level of every word or at the level of
>> sum of words?
>> it's both.
>> Yes, the model tries to predict general tendencies (means/proportions,
>> conditioned on factors). Observation sums are always involved.
>> The second answer is, if there is no factor conditioning, then that's it.
>> With factors, its measuring by item-differences if you want: And it seems
>> you want: Let 'word type' interact with 'treatment' to predict gain
>> proportions (0's and 1's for each participant). The relation between word
>> type and treatment is your domain, but your statistical issues imply the
>> following brms model: It measures differences between single
>> questions/words (pre-post on average/proportion on log-scale). If not, the
>> answer is still 'yes': Hierarchical Bayesian logistic regression can do
>> what you want :)). I assume, in terms of mixed-effect models (brm is not
>> mixed-effects least-mean squares, it's hierarchical Bayes) I would go with:
>>
>> mod5<-brm(score~Group*Time*Stims+(1|subject)+(1+Time*Stims| words),
>> family=bernoulli(),...)
>>
>> What you will get are Bayesian likelihood estimates for each
>> mean/proportion (and difference) across participants, for each Group by
>> every word. After estimation, you can apply traditional statistical
>> thinking like comparing mean estimates by word, or in general. These (both)
>> can be accessed via 'emmeans', or 'emmip', functions and you also can
>> compute Bayes Factors for a change vs no-change hypotheses (or others) for
>> each question/word (using the hypothesis() function, or prior_samples()
>> and posterior_samples() ), if you follow your own rules ;)
>> Best, René
>>
>> Am Sa., 23. März 2019 um 18:58 Uhr schrieb Souheyla GHEBGHOUB <
>> souheyla.ghebghoub using gmail.com>:
>>
>>> Dear René and any other member interested in this discussion,
>>>
>>> I appreciate the long feedback I received from you. But I can tell I
>>> could not well convey my concern.
>>>
>>> The aim of my analysis is to predict the odds of correct/incorrect
>>> responses based on some predictors: group (Treatment1/Treatment2/Control);
>>> Time (pretest/posttest), verbal frequency, concreteness vs.
>>> abstractness...etc. I have hypotheses such as, Time*Group interaction will
>>> show a significant effect as I expect students in Treatment Group 1 will
>>> outperform at the level of Posttest.
>>> I have a problem of predicting Response with Time as an Independent
>>> variable which could be summarised in this example of one participant:
>>>
>>> Subject word Group Pretest Posttest
>>> 1 1 control 0 1
>>> 1 2 control 0 1
>>> 1 3 control 1 1
>>> 1 4 control 1 1
>>> 1 5 control 1 0
>>> 1 6 control 0 0
>>> 1 7 control 0 0
>>> 1 8 control 0 0
>>> =3 =4
>>>
>>> Would the model compare the mean score of pretest against mean score of
>>> posttest? If yes, it will not truly reflect the data, because its not about
>>> the sum,
>>> its about every word on its own. So it should not be about 4 in posttest
>>> against 3 in pretest means 1 gain.
>>> In here, there are 2 gains in word 1,2 and one decline in word 5. But
>>> does the model analyse at the level of every word or at the level of sum of
>>> words?
>>> That’s the main problem I have not solved since 3 months now.
>>> Best,
>>> SOUHEYLA
>>>
>>> On Sat, 23 Mar 2019 at 12:46, René <bimonosom using gmail.com> wrote:
>>>
>>>> Hey, things clear up. Thanks for the picture. I want to challenge your
>>>> concerns:
>>>> First, you see that observing a main effect of PrePost (or "Time") in
>>>> mod3 can only mean two things: 1) performance goes up (positive effect); or
>>>> 2) performance goes down (negative effect), in general. Saying something
>>>> like this, is what the model is defined for. The conclusion will be this
>>>> "general increase" or "general decrease" or "no evidence here". If you have
>>>> a different question from that (i.e. on an item level), then you should
>>>> specify it in more detail (see below).
>>>> There are three issues in-between the lines of your questions:
>>>> 1. is it a statistical concern you have?
>>>> 2. is it an actual theoretical question you have?
>>>> 3. is it a matter of making a non-result to a result?
>>>>
>>>> 1. The statistical view:
>>>> Counter-question: who would ever assume that a 1 score for a word in a
>>>> pre-test will remain constant until eternity. And why should it? The answer
>>>> is: Nobody, and this is the reason statistics (probability theory) exists
>>>> :). So the first simple answer to your question is you do not need to test
>>>> whether observed gains from pre-to-post are 'genuine' (only from 0 to 1,
>>>> without decline cases) because 'nature' guarantees that there will be
>>>> decline somewhere. That's "randomness" :) But the question is, what's
>>>> stronger, gain or decline? See... and there is no problem in it: a general
>>>> main effect (e.g. an overall gain) still is an overall gain, even if some
>>>> cases decline.
>>>>
>>>> 2. Theoretical view:
>>>> If, however, such special item dynamics are theoretized in advance,
>>>> then simply test it :) For instance, the assumption whether the
>>>> treatment leads to gain on abstract words, but to decline on concrete
>>>> words, then should find into the model by coding the factor AbsCon
>>>> (abstract vs. concrete words) as fixed effect (assuming a within
>>>> participants manipulation).
>>>> mod4<-brm(score~PrePost*Group*AbsCon+(1+PrePost*AbsCon| subject) +
>>>> (1+group|words),...)
>>>> (Note: the 1+XX|subject just means random intercept for subject (1+)
>>>> plus slope for XX on subject; and writing (1+group|words) is the same
>>>> as writing (group|words), but you can estimate the slope without the
>>>> (word) intercept by writing (0+group|words))
>>>>
>>>> Without having such a theoretical account testing for can also be done
>>>> via:
>>>> mod4<-brm(score~PrePost*Group*words+(1+PrePost*words| subject),...)
>>>> But you will hardly be able to interpret the interactions in this model
>>>> because words alone has 28 levels.
>>>>
>>>> 3. But, your example seems to suggest a special case... (i.e. an actual
>>>> Null-Effect):
>>>> "If a participant has got 5 correct words out of 28 in both pretest
>>>> and posttest"
>>>> then there would be no general effect of PrePost in the model above
>>>> (generalizing to all participants now). And searching for "deeper" model
>>>> checks looks like rescuing all effects you can get (post hoc). But of
>>>> course, making an argument like there still is an effect, namely for 5
>>>> specific words, which is not observable because there is also a detriment
>>>> for other 5 words is possible, but requires a solid theory which explicitly
>>>> predicts this interaction, and an experiment which was explicitly designed
>>>> to test this interaction. (point 2)
>>>>
>>>> So... if it is point 2 you got... Then go ahead :) test it in a
>>>> meaningful way. Otherwise, simply treat this "effect by words" interaction
>>>> as random slope (1+group|words), or btw. (1+PrePost*Group|words) is also
>>>> possible..., which is basically 'statistically' integrating the idea that
>>>> the treatment*time effects vary (randomly) between stimuli. And doing this
>>>> in the random effects has the notion of "generalizing" estimation error in
>>>> the population, and should be preferred to implementing those in the fixed
>>>> effects, if the 28 words can be seen as a random (non-special) stimulus
>>>> sample. If this is not the case, then consider coding the "special" thing
>>>> about the words as fixed effects (e.g. if you want to use the same design
>>>> again, for testing something, while controlling for stimulus specifics).
>>>>
>>>> Best, René
>>>>
>>>> Am Sa., 23. März 2019 um 11:59 Uhr schrieb Souheyla GHEBGHOUB <
>>>> souheyla.ghebghoub using gmail.com>:
>>>>
>>>>> Dear René,
>>>>>
>>>>> Thank you for the feedback. Actually, my model was originally like you
>>>>> suggested now (except for slopes I had PrePost without 1 in both words and
>>>>> subjects. I called PrePost as "Time". I will read more about the 1+prepost
>>>>> form you mentioned.
>>>>>
>>>>> The reason why I gave up this model and looked for something else is
>>>>> the fact that each subject has 28 words tested twice (pre&post) and I was
>>>>> not too sure whether such model will take into consideration differences
>>>>> between pre & post at the level of every word of each participant (which is
>>>>> what I want) instead of merely comparing every participant's pre mean sores
>>>>> of 28 words against his post mean score (which is what i should avoid), here
>>>>> is a short example as to why:
>>>>>
>>>>> If a participant has got 5 correct words out of 28 in both pretest
>>>>> and posttest, there will be multiple interpretations: e.g. They could
>>>>> refer to the same words (i.e. 0 gain), or they could be totally new words
>>>>> (i.e. 5 gains) ...etc , hence I am not sure whether such model of pretest
>>>>> vs posttest will compare each subject score of each word from pretest to
>>>>> posttest then base its analysis on these score changes instead of comparing
>>>>> the sum scores between pre and post and which likely skew results.
>>>>>
>>>>> I posted about this in stackexchange 3 months ago and was told that it
>>>>> does compare word to word for every participant, but I am still not
>>>>> confident enough to use it because all the accurateness of the results and
>>>>> discussion chapters of my PhD thesis will be based on this decision.
>>>>>
>>>>> I look forward to receive feedback from you and any member reading
>>>>> this,
>>>>> Souheyla
>>>>> University of York
>>>>>
>>>>> On Sat, 23 Mar 2019, 10:01 René, <bimonosom using gmail.com> wrote:
>>>>>
>>>>>> Hi Souheyla,
>>>>>>
>>>>>> Well, I guess in your case it is simply more elegant to leave the
>>>>>> measured predictor out of the fixed effects, because there is also another
>>>>>> implied question (i.e. about the strength of change between pre and post).
>>>>>>
>>>>>> So, another possibility to re-define your model (as logistic
>>>>>> regression) allowing for better interpretations:
>>>>>> mod3<-brm(score~PrePost*Group+(1+PrePost | subject)+(1+group |
>>>>>> words),...)
>>>>>>
>>>>>> score = 0 or 1 for a given testitem
>>>>>> PrePost = Pre vs. Post (basically just an indicator of the
>>>>>> measurement time point)
>>>>>> Thus, the PrePost main effect will tell, whether there is a change
>>>>>> from pre to post (e.g. a gaint), and you can also tell how strong it is (in
>>>>>> odds ratios).
>>>>>> And if PrePost interacts with Group, then the change (e.g. a gain) is
>>>>>> moderated by the treatment, which seems to be your main question.
>>>>>>
>>>>>> Now in this model, you can also have by-subject random slopes for
>>>>>> PrePost of course (because the fixed effect of PrePost is present for every
>>>>>> subject).
>>>>>>
>>>>>> Best, René
>>>>>>
>>>>>>
>>>>>> Am Sa., 23. März 2019 um 10:12 Uhr schrieb Souheyla GHEBGHOUB <
>>>>>> souheyla.ghebghoub using gmail.com>:
>>>>>>
>>>>>>> I read that in multinomial regression, all independent variables
>>>>>>> should be
>>>>>>> variables that we manipulate. Can I still have pretest as IV without
>>>>>>> skewing my results?
>>>>>>>
>>>>>>> Best,
>>>>>>> Souheyla
>>>>>>>
>>>>>>> On Fri, 22 Mar 2019, 23:31 Souheyla GHEBGHOUB, <
>>>>>>> souheyla.ghebghoub using gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>> > Thank you both. I will look into this and see :)
>>>>>>> >
>>>>>>> > Best,
>>>>>>> > Souheyla
>>>>>>> >
>>>>>>> > On Fri, 22 Mar 2019, 22:02 Uanhoro, James, <
>>>>>>> uanhoro.1 using buckeyemail.osu.edu>
>>>>>>> > wrote:
>>>>>>> >
>>>>>>> >> In standard regression models, the assumption is predictor
>>>>>>> variables are
>>>>>>> >> measured without error. Test scores will have measurement error,
>>>>>>> hence
>>>>>>> >> Doran's comment when test scores are used as covariates. See:
>>>>>>> Hausman, J.
>>>>>>> >> (2001). Mismeasured Variables in Econometric Analysis: Problems
>>>>>>> from the
>>>>>>> >> Right and Problems from the Left. *Journal of Economic
>>>>>>> Perspectives*,
>>>>>>> >> *15*(4), 57–67. https://doi.org/10.1257/jep.15.4.57
>>>>>>> >> I will note that many practitioners ignore this issue, and it is
>>>>>>> quite
>>>>>>> >> common to use predictors measured with error. Consider the number
>>>>>>> of times
>>>>>>> >> people use polychotomized income measures, or SES measures as
>>>>>>> predictors,
>>>>>>> >> or some other "construct".
>>>>>>> >> On Mar 22 2019, at 5:39 pm, Souheyla GHEBGHOUB <
>>>>>>> >> souheyla.ghebghoub using gmail.com> wrote:
>>>>>>> >>
>>>>>>> >> Dear Doran,
>>>>>>> >>
>>>>>>> >> Could you explain more this point to me, please?
>>>>>>> >>
>>>>>>> >> Thank you,
>>>>>>> >> Souheyla
>>>>>>> >>
>>>>>>> >> On Fri, 22 Mar 2019, 21:19 Doran, Harold, <HDoran using air.org> wrote:
>>>>>>> >>
>>>>>>> >> Yes, but conditioning on the pre-test means you are using a
>>>>>>> variable
>>>>>>> >> measured with error and the estimates you obtain and now
>>>>>>> inconsistent, and
>>>>>>> >> that¹s a pretty big sin.
>>>>>>> >>
>>>>>>> >> On 3/22/19, 3:49 PM, "Souheyla GHEBGHOUB" <
>>>>>>> souheyla.ghebghoub using gmail.com>
>>>>>>> >> wrote:
>>>>>>> >>
>>>>>>> >> Dear René,
>>>>>>> >>
>>>>>>> >> Thank you for your feedback to me. You are right, dropping the
>>>>>>> pretest
>>>>>>> >> from
>>>>>>> >> covariate if I predict change definitely makes sense to me! But
>>>>>>> the fact
>>>>>>> >> that i need to control for the starting levels of participants
>>>>>>> makes it
>>>>>>> >> obligatory for me to chose the second way, which is predicting
>>>>>>> posttest
>>>>>>> >> instead of change to have pretest scores controlled for.
>>>>>>> >>
>>>>>>> >> You also chose (1+group | word) , which is new to me. Does it
>>>>>>> intend to
>>>>>>> >> assume the effect of group to vary across words, which is
>>>>>>> something
>>>>>>> >> applicable to my data, right?
>>>>>>> >> I will discuss all this with my supervisor, and may reply here
>>>>>>> again in
>>>>>>> >> few
>>>>>>> >> days if you do not mind.
>>>>>>> >> Thank you very much
>>>>>>> >> Souheyla
>>>>>>> >> University of York
>>>>>>> >>
>>>>>>> >>
>>>>>>> >> On Fri, 22 Mar 2019 at 13:42, René <bimonosom using gmail.com> wrote:
>>>>>>> >>
>>>>>>> >> Hi Souheyla,
>>>>>>> >>
>>>>>>> >> it seems to me that you will run into problems with your coding of
>>>>>>> >> change
>>>>>>> >> (gain, no gain and decline) because the 'change' is by
>>>>>>> >> definition/calculation depending on the predictor pretest.
>>>>>>> >> See, according to your coding scheme:
>>>>>>> >> Change = decline can only occur if pretest=1 (not by pretest=0).
>>>>>>> >> Change = gain can only occur if pretest = 0 (not by pretest=1)
>>>>>>> >> Change = No Gain can occur if pretest= 1 or 0
>>>>>>> >> In other words:
>>>>>>> >> If pretest = 1 then the possible outcomes can be decline or no
>>>>>>> gain
>>>>>>> >> If pretest = 0 then the possible outcomes can be gain or no gain
>>>>>>> >>
>>>>>>> >> And if the model result shows you then that the pre-test is
>>>>>>> >> significantly
>>>>>>> >> related to p(change-outcome), I guess there is no surprise in it,
>>>>>>> is it?
>>>>>>> >>
>>>>>>> >> So the first solution to this would be simply kicking the pre-test
>>>>>>> >> predictor out of the model completely, and predict:
>>>>>>> >> mod1 <- brm(Change ~ Group + (1|Subject) + (1+Group|Word),...)
>>>>>>> >> (Btw.: actually the first Hierarchical Bayes Model question I see
>>>>>>> on the
>>>>>>> >> mixed-effects mailing list :))
>>>>>>> >>
>>>>>>> >> Attempt for a further clarification on which random slopes would
>>>>>>> reflect
>>>>>>> >> the model's design:
>>>>>>> >> If you have a within-subjects design, by-subject random slopes are
>>>>>>> >> possible for the within-subject variable (e.g. if there are two
>>>>>>> sets of
>>>>>>> >> words/lists [e.g. abstract vs. concrete words] for each
>>>>>>> participant, and
>>>>>>> >> you test whether there is a performance-difference between these
>>>>>>> >> word-lists, then you can implement by-subject random slopes for
>>>>>>> words,
>>>>>>> >> because each participant has seen both sets.) If each participant
>>>>>>> has
>>>>>>> >> seen
>>>>>>> >> only one list (i.e. between subjects design) by subject random
>>>>>>> slopes
>>>>>>> >> for
>>>>>>> >> words are not appropriate, because there is no 'slope' by
>>>>>>> participant
>>>>>>> >> (i.e.
>>>>>>> >> by definition, having a slope requires at least two
>>>>>>> observations...).
>>>>>>> >> This
>>>>>>> >> is always a good rule of thumb without thinking about it too
>>>>>>> heavily :)
>>>>>>> >> Ans as you see: you can define a random slope for words:
>>>>>>> >> (1+Group|Word),
>>>>>>> >> because each word has been presented in each group (i.e. there
>>>>>>> can be a
>>>>>>> >> slope for each word). And intuitively speaking the
>>>>>>> Treatment-effect can
>>>>>>> >> vary depending on the stimuli you use, and the slope makes sense.
>>>>>>> (You
>>>>>>> >> also
>>>>>>> >> see in this example that the treatment effect can also vary by
>>>>>>> subjects,
>>>>>>> >> but in fact, this subject effect variation IS EQUAL to the effect
>>>>>>> you
>>>>>>> >> want
>>>>>>> >> to test, and having by subject group random slopes would
>>>>>>> eliminate the
>>>>>>> >> fixed effect...)
>>>>>>> >>
>>>>>>> >> Anyway, there is a second possibility to define your model,
>>>>>>> depending on
>>>>>>> >> how you want to interpret it. In the previous model you can say
>>>>>>> >> something
>>>>>>> >> about the type-of-change likelihoods depending on the treatment
>>>>>>> group.
>>>>>>> >> But
>>>>>>> >> you could implement the model as binomial as well (i.e. logistic
>>>>>>> >> regression)
>>>>>>> >>
>>>>>>> >> mod2 <- brm(posttest ~ pretest*Group + (1|Subject) +
>>>>>>> (1+Group|Word),...)
>>>>>>> >>
>>>>>>> >> And what you would expect here would be an interaction between
>>>>>>> pre-test
>>>>>>> >> and Group. For instance; if pretest=0 & treatment 1 then posttest
>>>>>>> larger
>>>>>>> >> than with pretest=0 & treatment 2; but not when pretest=1
>>>>>>> (because this
>>>>>>> >> is
>>>>>>> >> a plausible no gain situation). And so on...
>>>>>>> >> (And in this model there are no also no further random slopes
>>>>>>> hidden in
>>>>>>> >> your design :))
>>>>>>> >> Hope this helps.
>>>>>>> >>
>>>>>>> >> Best, René
>>>>>>> >>
>>>>>>> >>
>>>>>>> >> Am Do., 21. März 2019 um 14:01 Uhr schrieb Souheyla GHEBGHOUB <
>>>>>>> >> souheyla.ghebghoub using gmail.com>:
>>>>>>> >>
>>>>>>> >> Dear Philip,
>>>>>>> >>
>>>>>>> >> I understand , here is the structure of my data in case it could
>>>>>>> help.
>>>>>>> >>
>>>>>>> >> I have 3 groups of participants (control, treatment1,
>>>>>>> treatment2). Each
>>>>>>> >> group was tested twice, once before treatment (pretest) and once
>>>>>>> after
>>>>>>> >> treatment (posttest).
>>>>>>> >> In each test, they were tested on knowledge of 28 words, scores
>>>>>>> are
>>>>>>> >> dichotomous (0 = unknown , 1 = known). Tests are the same.
>>>>>>> >>
>>>>>>> >> I calculated change from pretest to posttest :
>>>>>>> >> if pretest 0 and posttest 0 = no gain
>>>>>>> >> if pretest 1 and posttest 1 = no gain
>>>>>>> >> if pretest 0 and posttest 1 = gain
>>>>>>> >> if pretest 1 and posttest 0 = decline
>>>>>>> >> So I ended up with a dependent variable called Change with 3
>>>>>>> levels
>>>>>>> >> (no_gain, gain, decline) and I tried to predict it using Group and
>>>>>>> >> Pretest
>>>>>>> >> as covariates using multinomial logit model. mod0 <- brm(Change ~
>>>>>>> >> Pretest
>>>>>>> >> +
>>>>>>> >> Group) I would like to add random effects for subjects but don't
>>>>>>> know
>>>>>>> >> what's the best form when Time factor is absent.
>>>>>>> >>
>>>>>>> >> I hope other statisticians who read this could help
>>>>>>> >> Thank you
>>>>>>> >> Souheyla
>>>>>>> >>
>>>>>>> >> [[alternative HTML version deleted]]
>>>>>>> >>
>>>>>>> >> _______________________________________________
>>>>>>> >> R-sig-mixed-models using r-project.org mailing list
>>>>>>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>> >>
>>>>>>> >>
>>>>>>> >>
>>>>>>> >> [[alternative HTML version deleted]]
>>>>>>> >>
>>>>>>> >> _______________________________________________
>>>>>>> >> R-sig-mixed-models using r-project.org mailing list
>>>>>>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>> >>
>>>>>>> >>
>>>>>>> >>
>>>>>>> >>
>>>>>>> >> [[alternative HTML version deleted]]
>>>>>>> >>
>>>>>>> >> _______________________________________________
>>>>>>> >> R-sig-mixed-models using r-project.org mailing list
>>>>>>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>> >>
>>>>>>> >>
>>>>>>>
>>>>>>> [[alternative HTML version deleted]]
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-sig-mixed-models using r-project.org mailing list
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>>>
>>>>>>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list