marginal effect lmer
Dear Andigo, On Mon, 23 Jul 2012 07:42:49 -0700 (PDT)
Andigo <Andreas.Goldberg at unige.ch> wrote:
Dear John, yes, the superscript shall say that the quadratic term of Catholicproportion is included as well (plus the interaction terms with it). In my dataset the quadratic Catholicproportion is already included as a variable of its own, so it should be fine or do I have to use the syntax with "poly..." you mentioned? If so how do I have to rewrite the command for the whole model?
It's hard to know where to begin, particularly since you deleted the earlier discussion. Retrieving that from my reply to you, the model you fit was
m1<- lmer(vote2011~ Catholic + Attendance+ logthreshold + West + Catholicproportion+ (Catholic * Catholicproportion) + (Attendance*Catholicproportion) + Catholicproportion?+ (Catholic *Catholicproportion?)+ (Attendance* Catholicproportion?) + (1 + Attendance+ Catholic | canton), data=dat2011, verbose = T)
I don't really know what these variables are but I'll assume that all are numeric variables and not, e.g., 0/1 dummy regressors that you put in the model instead of factors in the same way that you generated the quadratic regressor for Catholicproportion -- both of these are bad ideas in R because they disguise the structure of the model. Anyway, I suppose that what you want in the fixed-effects part of the model is logthreshold + West + (Catholic + Attendance)*poly(Catholicproportion, 2) This will give you an orthogonal quadradic in Catholicproportion that interacts both with Catholic and Attendance and will allow effect() to figure out the structure of the model; if you want a raw quadratic (which would provide the same fit to the data but be somewhat less numerically stable), you could substitute poly(Catholicproportion, 2, raw=TRUE).
Out of your description of the package I just tried the example you
mentioned for a lmer model, so I tried:
plot(effect("Catholic:Catholicproportion", m1), grid=TRUE)
I'm not sure what you read, but what you got is what you asked for; with the respecified model, you could ask for the "Catholic:poly(Catholicproportion, 2)" effect, or more simply plot(allEffects(m1)) to plot all high-order fixed-effects terms in the model.
The result are ten (dk why 10 and not just 1) little plots, which look rather linear and not the curvilinear pattern it should be.
Catholic and Catholicproportion are presumably numeric variables, and as documented in ?effect, each will by default be set to 10 levels spanning the range of the variable. There are then 10*10 == 100 combinations of values of these variables at which the effect will be estimated. That you expected one plot suggests to me that you don't quite understand what effect() computes (though all the lines could be put on one plot by setting the argument multiline=TRUE -- ?plot.eff). The lines are straight because effect() couldn't read your mind to know that Catholicproportion2 and Catholicproportion are related, which is why using poly() instead is a good idea.
So far I always calculated the marginal effects in Stata, which works fine for rather simple hierarchical models. But when the models become more complicated (more random effects, more interaction terms,..) Stata often fails to calculate the models, so that I cannot use it for displaying the marginal effects. Thats why I try to find a solution in R to calculate the marginal effects exactly the same way as I do in Stata. In Stata I followed the syntax by Brambor, Clark and Golder (https://files.nyu.edu/mrg217/public/interaction.html). Now I just wonder if there is a way to calculate the marginal effects in R exactly the same way as they do in Stata?
I already explained that effect() doesn't compute marginal effects. There may well be an R package that does, but I'm not aware of one. It should, however, be a simple matter to compute marginal effects from the lmer() results if you find marginal effects interesting. Best, John
Best, Andigo -- View this message in context: http://r.789695.n4.nabble.com/marginal-effect-lmer-tp4637421p4637452.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.