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
Variance function plot  get CIs in R?

 Posts: 1279
 Joined: Mon Oct 19, 2009 10:34 am
Re: Variance function plot  get CIs in R?
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 level2 random part parameters and covariance matrix
pos < grep("RP2", names(coef(mymodel)))
RP2 < coef(mymodel)[pos]
RP2cov < vcov(mymodel)[pos, pos]
# Reshape level2 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))