var.components
Posted: Mon Jul 23, 2012 1:06 pmI have two questions regarding the var.components function in Rmark.
1. Why is there an upper limit on the range searched to find the root for equation 2 in Burnham and White (2002)? I find that if I simulate data with normal variability in survival on the logit scale, the solution to this equation is above max(vcv), the upper limit in var.components function. If I increase this upper limit, I get results that agree (almost exactly, see Q2) with Program Mark. I realize that Burnham and White assumed normal variability on the real scale, but I should still be able to fit their model to my simulated data.
2. I've extended var.components to calculate the confidence interval on the process variance, the predicted random effects, and their confidence intervals. It would also be easy to compute the AIC for the random effects model. Everything appears to match up between the extended function and Program Mark except the confidence interval on the process variance in my function is consistently very slightly more narrow. The ratio between the interval in Program Mark to R is ~1:0.995. Any ideas why this could be? I find it strange that the point estimate of the process variance matches up between the two programs, but the confidence interval does not exactly. The only difference between these two tasks is that the confidence intervals involve the chi squared quantile function. Here is the code I'm working with (my additions are the added upper bound for the process variance and the section in the middle separated by comments where I started and stopped editing.
1. Why is there an upper limit on the range searched to find the root for equation 2 in Burnham and White (2002)? I find that if I simulate data with normal variability in survival on the logit scale, the solution to this equation is above max(vcv), the upper limit in var.components function. If I increase this upper limit, I get results that agree (almost exactly, see Q2) with Program Mark. I realize that Burnham and White assumed normal variability on the real scale, but I should still be able to fit their model to my simulated data.
2. I've extended var.components to calculate the confidence interval on the process variance, the predicted random effects, and their confidence intervals. It would also be easy to compute the AIC for the random effects model. Everything appears to match up between the extended function and Program Mark except the confidence interval on the process variance in my function is consistently very slightly more narrow. The ratio between the interval in Program Mark to R is ~1:0.995. Any ideas why this could be? I find it strange that the point estimate of the process variance matches up between the two programs, but the confidence interval does not exactly. The only difference between these two tasks is that the confidence intervals involve the chi squared quantile function. Here is the code I'm working with (my additions are the added upper bound for the process variance and the section in the middle separated by comments where I started and stopped editing.
- Code: Select all
var.components2=function (theta, design, vcv, alpha, upper, LAPACK = TRUE)
{
if (nrow(design) != length(theta))
stop("Number of rows of design matrix must match length of theta vector")
if (nrow(vcv) != length(theta))
stop("Number of rows of vcv matrix must match length of theta vector")
if (ncol(vcv) != length(theta))
stop("Number of columns of vcv matrix must match length of theta vector")
if (length(theta) <= ncol(design))
stop("Length theta must exceed number of columns of design")
rn = (1:nrow(design))[apply(design, 1, function(x) any(as.numeric(x) !=
0))]
theta = theta[rn]
design = design[rn, , drop = FALSE]
vcv = vcv[rn, rn]
sigma = 0
beta.hat = function(Dinv, X, theta) return(solve(qr(t(X) %*%
Dinv %*% X, LAPACK = LAPACK)) %*% (t(X) %*% Dinv %*%
theta))
compute.Dinv = function(sigma, vcv) {
K = nrow(vcv)
sigmamat = matrix(0, nrow = K, ncol = K)
diag(sigmamat) = sigma
return(solve(qr(vcv + sigmamat, LAPACK = LAPACK)))
}
mom.sig = function(sigma, vcv, X, theta) {
Dinv = compute.Dinv(sigma, vcv)
beta = beta.hat(Dinv, X, theta)
return(t(theta - X %*% beta) %*% Dinv %*% (theta - X %*%
beta) - (nrow(X) - length(beta)))
}
soln = try(uniroot(mom.sig, lower = -(min(Re(eigen(vcv)$values)) +
1e-12), upper = upper, vcv = vcv, theta = theta, X = design,
tol = 1e-15))
if (class(soln) == "try-error")
sigma = 0
else if (soln$root < 0)
sigma = 0
else sigma = soln$root
Dinv = compute.Dinv(sigma, vcv)
beta = beta.hat(Dinv, design, theta)
#Begin Augustine additions
eig=eigen(Dinv)
Dinvsqrt=eig$vectors %*% diag(sqrt(eig$values)) %*% qr.solve(eig$vectors)
H=sqrt(sigma)*Dinvsqrt
X=rep(1,length(theta))
A=X%*%qr.solve(t(X)%*%Dinv%*%X)%*%t(X)
I=diag(1,length(theta))
G=H+(I-H)%*%A%*%Dinv
betarand=G%*%theta
VCbetarand=G%*%vcv%*%t(G)
RMSE=sqrt(diag(VCbetarand)+(betarand-theta)^2)
betarandL=betarand-qnorm(1-alpha/2,0,1)*RMSE
betarandU=betarand+qnorm(1-alpha/2,0,1)*RMSE
mom.sigCI = function(sigma, vcv, X, theta,p) {
Dinv = compute.Dinv(sigma, vcv)
beta = beta.hat(Dinv, X, theta)
return(t(theta - X %*% beta) %*% Dinv %*% (theta - X %*%
beta) - qchisq(p,nrow(X) - length(beta)))
}
solnL = try(uniroot(mom.sigCI, lower = -(min(Re(eigen(vcv)$values)) +
1e-12), upper = upper, vcv = vcv, theta = theta, X = design, p=1-alpha/2,
tol = 1e-15))
solnU = try(uniroot(mom.sigCI, lower = -(min(Re(eigen(vcv)$values)) +
1e-12), upper = upper, vcv = vcv, theta = theta, X = design, p=alpha/2,
tol = 1e-15))
betarand=cbind(betarand,RMSE,betarandL,betarandU)
colnames(betarand)=c("est.","RMSE","lower","upper")
sigma = c(sigma,solnL$root,solnU$root)
names(sigma)=c("est.","lower","upper")
#end Augustine additions
rownames(beta) = colnames(design)
vcv.beta = solve(qr(t(design) %*% Dinv %*% design, LAPACK = LAPACK))
rownames(vcv.beta) = colnames(design)
colnames(vcv.beta) = colnames(design)
beta = as.data.frame(beta)
beta$se = sqrt(diag(vcv.beta))
names(beta) = c("Estimate", "SE")
return(list(sigma = sigma, beta = beta,betarand=betarand,vcv.beta = vcv.beta))
}