Variance function plot - get CIs in R?

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
andrewjdbell
Posts: 17
Joined: Mon Jun 03, 2013 3:19 pm

Variance function plot - get CIs in R?

Post by andrewjdbell »

Hi there,

I want to produce variance function plots in R, of the sort that you can get using the "Model - variance function" tool in MLwiN. Crucially, in MLwiN you can get confidence intervals for those variance functions. Whilst I can work out the variance function in R (var(Int) + 2*cov*X + var(slope)*Xsquared), any idea the best way of calculating confidence intervals around that line?

Thanks!
Andy

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

Re: Variance function plot - get CIs in R?

Post by ChrisCharlton »

The calculations that MLwiN does to calculate the variance function and s.e. are something like the following (note that there could be some errors so please compare the results with those in MLwiN):

Code: Select all

# Run a model
mymodel <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), data=tutorial)

# Extract the level-2 random part parameters and co-variance matrix
pos <- grep("RP2", names(coef(mymodel)))
RP2 <- coef(mymodel)[pos]
RP2cov <- vcov(mymodel)[pos, pos]

# Reshape level-2 random parameters into matrix form
RPmat <- matrix(0, 2, 2)
RPmat[lower.tri(RPmat, diag=TRUE)] <- RP2
RPmat[upper.tri(RPmat)] <- t(RPmat)[upper.tri(RPmat)]

# Calculate variance for values in the first row
x <- as.matrix(tutorial[1, c("cons", "standlrt")])
x %*% RPmat %*% t(x)

# Calculate variance for values in the third row
x <- as.matrix(tutorial[3, c("cons", "standlrt")])
x %*% RPmat %*% t(x)

# Calculate variance s.e. for values in the first row
x <- as.matrix(tutorial[1, c("cons", "standlrt")])
# Create terms corresponding to the random part parameters
xrp <- matrix(0, 1, 3)
xrp[1, 1] <- x[1,1] * x[1,1]
xrp[1, 2] <- x[1,1] * x[1,2] * 2
xrp[1, 3] <- x[1,2] * x[1,2]
sqrt(xrp %*% RP2cov %*% t(xrp))

Post Reply