[R-sig-ME] question about random coefficients model (RCM)
Malcolm Fairbrother
m.fairbrother at bristol.ac.uk
Thu May 5 18:26:20 CEST 2011
Dear Jose,
Did you ever figure out how to replicate Western's (1998) model using lme4?
I've figured out how to replicate precisely the "Interactions Model" in his Table 2, but not the "Hierarchical Model"--see code below. I'm not sure if the issue is the estimation technique (REML rather than a Gibbs sampler), or if I've just not understood the model. I've also tried using MCMCglmm but am not having much more luck with that--other than the intercept, all fixed effects span 0, unlike in the results that Western reports.
I'd be interested to know if you found a way to replicate Western's results more precisely.
Cheers,
Malcolm
library(pcse) # includes the relevant dataset
data(agl)
agl2 <- agl[agl$year!=1970,]
agl2$lgrowth <- as.vector(unlist(by(agl, agl$country, function(x) x$growth[-15]))) # create lagged y variable
for (i in c(5:8, 11)) { agl2[,i] <- as.vector(unlist(by(agl2, agl2$country, function(x) scale(x[,i], scale=F)))) } # to center the time-series predictors within countries
w98.eq11 <- lm(growth ~ lgrowth*central + opengdp*central + openex*central + openimp*central + leftc*central, data=agl2) # matches the "Interactions Model" in Table 2
library(lme4)
w98.eq13 <- lmer(growth ~ lgrowth*central + opengdp*central + openex*central + openimp*central + leftc*central + (lgrowth + opengdp + openex + openimp + leftc | country), data=agl2)
data.frame(round(fixef(w98.eq13), 3), round(coef(w98.eq11), 3)); ranef(w98.eq13)
library(MCMCglmm)
prior <- list(R = list(V = 1, nu = 0), G = list(G1 = list(V = diag(6), nu = 1)))
w98.eq13.MCMC <- MCMCglmm(growth ~ lgrowth*central + opengdp*central + openex*central + openimp*central + leftc*central, random=~us(1 + lgrowth + opengdp + openex + openimp + leftc):country, data=agl2, prior=prior)
> Message: 5
> Date: 30-Jul-2010 16:48:33 EDT
> From: JOSE A ALEMAN <aleman at fordham.edu>
> To: Andrew Dolman <andydolman at gmail.com>
> Cc: dmbates at gmail.com, r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] question about random coefficients model (RCM)
> Message-ID: <OF4E9738AC.B8F4A4AE-ON85257770.00710104 at fordham.edu>
> Content-Type: text/plain; charset=ISO-8859-1
>
> Hi Andy,
>
> You're right, ranef() does actually work. Actually, data.neocorporatism is
> pretty much a time-invariant institutional variable. What I want is to
> estimate a random intercepts random slopes model where the effects of other
> time varying variables enter the model through the institutional variable
> "data.neocorporatism". The effects for neocorporatism, in other words,
> should be different depending on the country, that is, there shouldn't be a
> single coefficient for this variable but 18 different coefficients for the
> 18 different countries. If that's the case, then "data.neocorporatism"
> should not be in the fixed effects equation.
>
> I indeed tried the model having removed this variable, but the result after
> typing coef() was the same. Actually, ranef() gave me what I believe are
> the random intercepts for each country, but I still need to figure out a
> way to get the coefficients.
>
> I'm probably doing something wrong, but I am basically trying to fit a
> model similar to Western's (1998) described in my original message below.
>
> Thank you,
>
> Jose
>
>
>
> Andrew Dolman
> <andydolman at gmail
> .com> To
> JOSE A ALEMAN <aleman at fordham.edu>
> 07/30/2010 03:45 cc
> PM Douglas Bates
> <bates at stat.wisc.edu>,
> dmbates at gmail.com,
> r-sig-mixed-models at r-project.org
> Subject
> Re: [R-sig-ME] question about
> random coefficients model (RCM)
>
>
>
>
>
>
>
>
>
>
> Hi Jose
>
> data.nation is not present in the the fixed effects specification of
> your model, so coef() won't work, but ranef() should. Your model
> assumes that there is no overall relationship between data.gini_net
> and data.nation and yet allows for a relationship within the
> sub-groups defined by data.neocorporatism_std, with a slope that
> varies randomly between groups but is zero overall. That's probably
> not what you want.
>
> I also notice that data.neocorporatism seems to have been standardised
> somehow, which would be odd for a categorical variable used to group
> things.
>
>
>
> andydolman at gmail.com
>
>
>
> On 30 July 2010 18:52, JOSE A ALEMAN <aleman at fordham.edu> wrote:
>> The model I'm using is of the form:
>>
>> RCM <- lmer (data.gini_net ~ data.neocorporatism_std +
>> data.firm_level_coop_std + data.cog + data.gov_employ_std +
>> ? ? ?data.netden_std + data.trade_openness_std +
>> data.capital_liberalization_std + data.lagged_unemploy_std +
>> ? ? ?data.deindustrialization_std + data.pop_over_65_percent_std +
>> (data.nation|data.neocorporatism_std), data=new.data)
>>
>> but when I try to use the extractor coef(), I get the following error
>> message:
>>
>>> coef(RCM)
>> Error in coef(RCM) : unable to align random and fixed effects
>>
>>
>>
>>
>> ? ? ? ? ? ? Douglas Bates
>> ? ? ? ? ? ? <bates at stat.wisc.
>> ? ? ? ? ? ? edu> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? To
>> ? ? ? ? ? ? Sent by: ? ? ? ? ? ? ? ? ?JOSE A ALEMAN <aleman at fordham.edu>
>> ? ? ? ? ? ? dmbates at gmail.com ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?cc
>> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? r-sig-mixed-models at r-project.org
>> ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Subject
>> ? ? ? ? ? ? 07/17/2010 11:30 ? ? ? ? ?Re: [R-sig-ME] question about
>> ? ? ? ? ? ? AM ? ? ? ? ? ? ? ? ? ? ? ?random coefficients model (RCM)
>>
>>
>> On Sat, Jul 17, 2010 at 9:52 AM, JOSE A ALEMAN <aleman at fordham.edu>
> wrote:
>>>
>>> Dear list users,
>>>
>>> I have a question about how to estimate the RCM described by Western
>> (1998)
>>> in his article "Causal Heterogeneity in Comparative Research: A Bayesian
>>> Hierarchical Modeling Approach". While I know how to estimate this model
>> in
>>> R via restricted maximum likelihood, I am at a loss as to how to
>> interpret
>>> the results. The random coefficients are supposed to be different for
>> each
>>> country, yet R provides ?a standard regression-like output with a single
>>> coefficient for the variable with the random slope. How would I go about
>>> extracting the country specific information?
>>
>> Can you be more specific about how you fit the model? ?Just saying "in
>> R" is too general.
>>
>> If you fit the model with lmer from the lme4 package then you may want
>> to look at the output from the ranef() and coef() extractors.
More information about the R-sig-mixed-models
mailing list