www.cmm.bristol.ac.uk/forum

Posted: **Sat Aug 15, 2015 4:06 pm**

Hi
I would like to have some feedback on some alternative parameterizations: one using an intercept and one omitting it. My question is a mixture of both Mlwin and R2mlwin.

My model: I want to estimate the level of economic orientation in journals during three eras. I have therefore constructed three dummies: Era_1890to1920 + Era_1921to1985 + Era_1986to2014, which cover all my observations.

The model I want to estimate (preferably in R2mlwin):

# removing the intercept
(Model1 <- runMLwiN(OrgSocThy~ Era_1890to1920 + Era_1921to1985 + Era_1986to2014 + (Era_1890to1920 + Era_1921to1985 + Era_1986to2014|journal) + (1 | Article_ID), estoptions = list(EstM = 0, debugmode =F, resi.store=T, x64=F, optimat=T), data = df.eo.sorted))

# Result: gives This gives warning about random params has gone –ve definite. But it gives sensible numbers.

[see attached file for output]

But I have these two alternative parameterizations:

# Keeping the intercept as is with two dummies
(Model2 <- runMLwiN(OrgSocThy~ 1+ Era_1921to1985 + Era_1986to2014 + (1+ Era_1921to1985 + Era_1986to2014|journal) + (1 | Article_ID), estoptions = list(EstM = 0, debugmode =F, resi.store=T, x64=F, optimat=T), data = df.eo.sorted))

# Result: This gives warning about random params has gone –ve definite

[see attached file for output]

# keeping the 1 in, 1+ Era_1890to1920 + Era_1921to1985 + Era_1986to2014|journal)
(Model3 <- runMLwiN(OrgSocThy~ Era_1 + Era_1921to1985 + Era_1986to2014 + (1+ Era_1890to1920 + Era_1921to1985 + Era_1986to2014|journal) + (1 | Article_ID), estoptions = list(EstM = 0, debugmode =F, resi.store=T, x64=F, optimat=T), data = df.eo.sorted))

# Result: This gives warning about random params has gone –ve definite

[see attached file for output]

My question
1. I know that by removing the intercept one can directly estimate the values of the three eras. (dummies), but why is model 1 giving sensible numbers where as models 2 and 3 are not. They are just alternative parameterizations of the same model?
2. Moroever, in model 3, what does it mean when I keep the constant in (1 + Era_1986to2014|journal)? I am not entirely sure, but I assume that when I remove the intercept I should also remove the “1” in the randome part. Please, correct me if I am wrong?

question with output picture.docx: (195.75 KiB) Downloaded 579 times

Thanks in advance
Adel

Posted: **Mon Aug 17, 2015 4:23 pm**

Hi - I've tried to roughly replicate your scenario with one of the sample datasets (tutorial):

Code: Select all

library(R2MLwiN)
data("tutorial")
my_tutorial <- tutorial

# create new factor variable standlrt_cat from standlrt, consisting of low/mid/high
my_tutorial$standlrt_cat[my_tutorial$standlrt <  -1] <- "low"
my_tutorial$standlrt_cat[my_tutorial$standlrt >=  -1 & my_tutorial$standlrt < 1] <- "mid"
my_tutorial$standlrt_cat[my_tutorial$standlrt >=  1] <- "high"
my_tutorial$standlrt_cat <- as.factor(my_tutorial$standlrt_cat)

# use R2MLwiN's Untoggle function to create dummy for each category of standlrt_cat
my_tutorial = cbind(my_tutorial, Untoggle(my_tutorial$standlrt_cat, 'standlrt_cat'))

# fit model without constant as intercept
(Model1 <- runMLwiN(normexam ~ standlrt_cat_low + standlrt_cat_mid + standlrt_cat_high + (standlrt_cat_low + standlrt_cat_mid + standlrt_cat_high|school) + (1 | student), estoptions = list(EstM = 0, debugmode =F, resi.store=T, x64=F, optimat=T), data = my_tutorial))
#--------------------------------------------------------------------------------------------------- 
#The model formula:
#normexam ~ standlrt_cat_low + standlrt_cat_mid + standlrt_cat_high + 
#    (standlrt_cat_low + standlrt_cat_mid + standlrt_cat_high | 
#        school) + (1 | student)
#Level 2: school     Level 1: student      
#--------------------------------------------------------------------------------------------------- 
#The fixed part estimates:  
#                       Coef.   Std. Err.        z    Pr(>|z|)         [95% Conf.   Interval] 
#standlrt_cat_low    -0.83171     0.04933   -16.86   8.723e-64   ***     -0.92839    -0.73503 
#standlrt_cat_mid    -0.01905     0.04418    -0.43      0.6663           -0.10564     0.06754 
#standlrt_cat_high    0.82042     0.06390    12.84   9.788e-38   ***      0.69519     0.94565 
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
#--------------------------------------------------------------------------------------------------- 

# fit model with constant as intercept
(Model2 <- runMLwiN(normexam ~ 1 + standlrt_cat_mid + standlrt_cat_high + (1 + standlrt_cat_mid + standlrt_cat_high|school) + (1 | student), estoptions = list(EstM = 0, debugmode =F, resi.store=T, x64=F, optimat=T), data = my_tutorial))
#--------------------------------------------------------------------------------------------------- 
#The model formula:
#normexam ~ 1 + standlrt_cat_mid + standlrt_cat_high + (1 + standlrt_cat_mid + 
#    standlrt_cat_high | school) + (1 | student)
#Level 2: school     Level 1: student      
#--------------------------------------------------------------------------------------------------- 
#The fixed part estimates:  
#                       Coef.   Std. Err.        z     Pr(>|z|)         [95% Conf.   Interval] 
#Intercept           -0.83173     0.04934   -16.86    9.197e-64   ***     -0.92843    -0.73503 
#standlrt_cat_mid     0.81269     0.04597    17.68    6.013e-70   ***      0.72259     0.90278 
#standlrt_cat_high    1.65216     0.06566    25.16   1.059e-139   ***      1.52347     1.78086 
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
#---------------------------------------------------------------------------------------------------

The estimates from these models indicate they are reparameterisations of each other:

Model 1
low -0.83171
mid -0.01905
high 0.82042

Model 2
low -0.83173
mid -0.83173 + 0.81269 = -0.01904
high -0.83173 + 1.65216 = 0.82043

I noticed from your output that not all of your models iterated to convergence, however (as indicated towards the top of the model output), and I wonder if that might be causing some confusion re: the estimates returned? The maximum number of iterations before IGLS estimation halts can be altered via maxiter in estoptions (see ?runMLwiN).

Hope that helps.

Posted: **Tue Aug 18, 2015 2:52 pm**

Hi Richard,

Thanks for coming back to me. This makes sense, and part of the confusion arouse from me not noticing that the models did not converge.

As a follow up and as a response to my second question:

When one fits models of the type of model1, one needs to always remove the “1” in the “(1|school)”. I see that, this model do not converge, right?

Thanks
Adel

Posted: **Tue Aug 18, 2015 4:01 pm**

Hi - yes, that was my understanding: I don't think it makes sense to try to estimate an overall intercept variance (and its covariances) at level 2 when the coefficients of all the dummy variables of a categorical predictor (i.e. all three dummies from a categorical variable consisting of three categories) are also allowed to vary at that level. In terms of the categorical predictor, the random part of the model at level 2 would already be saturated (there'd be nothing left to estimate for it) when you include all its dummies.

www.cmm.bristol.ac.uk/forum

Removing the (random) intercept

Removing the (random) intercept

Re: Removing the (random) intercept

Re: Removing the (random) intercept

Re: Removing the (random) intercept