Page 1 of 1

extracting variance coefficient

Posted: Thu Sep 16, 2021 1:55 pm
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! :)

Re: extracting variance coefficient

Posted: Tue Sep 21, 2021 11:37 am
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.

Re: extracting variance coefficient

Posted: Wed Sep 22, 2021 12:51 pm
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

Re: extracting variance coefficient

Posted: Wed Sep 22, 2021 1:37 pm
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.

Re: extracting variance coefficient

Posted: Wed Sep 22, 2021 1:42 pm
by med151991
Thats great thanks so much!

Re: extracting variance coefficient

Posted: Mon Oct 04, 2021 8:30 am
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!

Re: extracting variance coefficient

Posted: Tue Oct 05, 2021 9:51 am
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.

Re: extracting variance coefficient

Posted: Tue Oct 05, 2021 11:26 am
by med151991
Thanks so much Chris thats much clearer to me thank you