extracting variance coefficient

Welcome to the forum for R2MLwiN users. Feel free to post your question about R2MLwiN here. The Centre for Multilevel Modelling take no responsibility for the accuracy of these posts, we are unable to monitor them closely. Do go ahead and post your question and thank you in advance if you find the time to post any answers!

Go to R2MLwiN: Running MLwiN from within R >> http://www.bris.ac.uk/cmm/software/r2mlwin/
Post Reply
med151991
Posts: 13
Joined: Tue Dec 11, 2018 11:52 am

extracting variance coefficient

Post by med151991 »

Hello,

I am trying to calculate the median incidence rate ratio for my poisson three level model. I have a calculation and code to do this, however its based on either the glmer or brm package in R. I was wondering if someone might help me apply this to the functions/terms in MLwiN?

this is the code for brm

Code: Select all

tau.info <- summary(data)$random$region
tau <- tau.info[1,c("Estimate","l-95% CI","u-95% CI")]
tau2CI <- tau^2
MRR <- exp(sqrt(2*tau2CI)*qnorm(0.75))
then after estimating the model

Code: Select all

estMRR(mymodel) 
returns the result. I assume random can be used in conjuction with brm, but it does not work when this function and estMRR after running my R2MLwiN model. I tried the following, which did not produce any error messages but it says null when I return the results.

Code: Select all

tau.info <- summary(mymodel)["RP2_var_Intercept"]
tau <- tau.info[1,c("Estimate","l-95% CI","u-95% CI")]
tau2CI <- tau^2
MRR <- exp(sqrt(2*tau2CI)*qnorm(0.75))
this version of the code is from glmmer from a paper by Austin (2018)

Code: Select all

tau.info <‐ summary(mlm)$varcor$hospital
but this would not work for MLwiN because varcor is not a function or slot used in that package. If anyone has any insight into the solution I would greatly appreciate it! :)
ChrisCharlton
Posts: 1285
Joined: Mon Oct 19, 2009 10:34 am

Re: extracting variance coefficient

Post by ChrisCharlton »

R2MLwiN does not currently provide an implementation of the VarCorr() function, however it does return these values within the output of coef(). The following shows the same model being fitted with glmer, brms, MLwiN(IGLS) and MLwiN(MCMC):

Code: Select all

data(mmmec, package = "R2MLwiN")

library(lme4)
glmermod <- glmer(obs ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region), data = mmmec, family = poisson(link = "log"))

library(brms)
brmmod <- brm(obs ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region), data = mmmec, family = poisson)

library(R2MLwiN)
iglsmod <- runMLwiN(log(obs) ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region), D = "Poisson", data = mmmec)
mcmcmod <- runMLwiN(log(obs) ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region), D = "Poisson", estoptions = list(EstM = 1), data = mmmec)
In glmer and brms you could then extract the random part parameter estimates using the VarCorr() function, i.e.

Code: Select all

> VarCorr(glmermod)
 Groups Name        Std.Dev.
 region (Intercept) 0.23970 
 nation (Intercept) 0.40623 

> VarCorr(brmmod)
$nation
$nation$sd
           Estimate Est.Error      Q2.5     Q97.5
Intercept 0.5200079 0.1700911 0.2905706 0.9578975


$region
$region$sd
           Estimate  Est.Error      Q2.5     Q97.5
Intercept 0.2454637 0.02715162 0.1987185 0.3031363
However to get the same results from a R2MLwiN model you have to take a subset of the results from the coef function:

Code: Select all

> sqrt(coef(iglsmod)[c("RP3_var_Intercept", "RP2_var_Intercept")])
RP3_var_Intercept RP2_var_Intercept 
        0.4350172         0.2129688 

> sqrt(coef(mcmcmod)[c("RP3_var_Intercept", "RP2_var_Intercept")])
RP3_var_Intercept RP2_var_Intercept 
        0.4829251         0.2452319
Or from the MCMC version calculate this from the chains directly:

Code: Select all

> colMeans(sqrt(mcmcmod@chains[, c("RP3_var_Intercept", "RP2_var_Intercept")]))
RP3_var_Intercept RP2_var_Intercept 
        0.4618277         0.2436698
Which would also allow you to easily calculate the credible intervals too.

Note that MLwiN returns the variances directly, so for your purposes you can remove the square root to cancel out the later squaring in your calculation.
med151991
Posts: 13
Joined: Tue Dec 11, 2018 11:52 am

Re: extracting variance coefficient

Post by med151991 »

Hi Chris,

Thank you so much thats so helpful!

I have a couple other questions too. Firstly, how do you build up a model from previous values in R? And also if you run the model initially in IGLS do you have to specificy something when running then in MCMC about using those starting values or is it automatic?

And finally, is there a way to let the random variance at level 2 go negative? When I run the empty model with just the variance terms i get a variance of zero for the second level, but in the full model which I have run for many iterations I get a small variance. I wondering if it could also be because some of those cluster sizes are small numbers as its based on household?

Many thanks,

Megan
ChrisCharlton
Posts: 1285
Joined: Mon Oct 19, 2009 10:34 am

Re: extracting variance coefficient

Post by ChrisCharlton »

You can find an example of setting starting values for parameters in section 5.5 of the MCMC manual replication materials. If you do not specify these then the model will be run with (R)IGLS (up to the specified maximum number of iterations) to generate starting values before the MCMC run begins.

For (R)IGLS models you can allow variance estimates at chosen levels to be negative via the reset option. This is specified as a vector with an element for each level. A value of zero causes variances and associated covariances that go negative to be reset to zero, a value of one will set negative variances to zero, but not the associated covariances, a value of two allows any variances or covariances to become negative at that level.

See the R2MLwiN reference files for further details.
med151991
Posts: 13
Joined: Tue Dec 11, 2018 11:52 am

Re: extracting variance coefficient

Post by med151991 »

Thats great thanks so much!
med151991
Posts: 13
Joined: Tue Dec 11, 2018 11:52 am

Re: extracting variance coefficient

Post by med151991 »

Hi Chris sorry to bother you again!

I am getting myself a little confused about the IGLS starting values. From everything I have read, it seems to indicate that when running the model if you have EstM=1, this generates starting values from IGLS, is this right? But I have seen some comments about putting the IGLS values into the MCMC model, which I wondered if there is something that is specifcally written in the code but cannot find an example of that. So far I have been running the IGLS model first then MCMC but is this necessary?

Code: Select all

model<-log(death)~1+gender+education+offset(log(personyears))+(1|region)+(1|household)
modeligls<-runMLwiN((Formula=model, D="Poisson", data=dataset)
modelmcmc<runMLwiN(Formula=model, D="Poisson", data=dataset, estoptions=list(EstM=1, mcmcMeth=list(burnin=1000, nchains=4, iterations=100000)))
Also, can you use the reset option in an MCMC model formula? I have one variance that goes negative during just the IGLS model, but if I am writing out the mcmc code I am now thinking the starting values will be different as i haven't input any of the values from the IGLS model I have run before that.

Thank you!
ChrisCharlton
Posts: 1285
Joined: Mon Oct 19, 2009 10:34 am

Re: extracting variance coefficient

Post by ChrisCharlton »

If you are running the model with MCMC (EstM=1) and do not specify starting values then it runs IGLS first and uses the result from this as starting values. Otherwise it uses the values that you provide. You can see this in lines 1496-1520 of https://github.com/r-forge/r2mlwin/blob ... ite.MCMC.R.

If you wanted to see/check the starting values, or have more control over the model being fitted then you could run IGLS first manually, otherwise this is not necessary.

The MCMC estimation option does not use the reset option and this is ignored for the IGLS model that is used to generate starting values.
med151991
Posts: 13
Joined: Tue Dec 11, 2018 11:52 am

Re: extracting variance coefficient

Post by med151991 »

Thanks so much Chris thats much clearer to me thank you
Post Reply