Variance function plot - get CIs in R?

[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]

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))