var.components

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

var.components

Postby benaug » Mon Jul 23, 2012 1:06 pm

I 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.

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))
}
benaug
 
Posts: 6
Joined: Fri Jul 20, 2012 1:27 pm

Re: var.components

Postby jlaake » Mon Jul 23, 2012 11:15 pm

There is an upper limit because the root finding function needs a lower and upper limit that contains the zero. Can't remember now what I chose but probably chose something that seemed to be more than what you would see in practice. I leave for the field tomorrow but can look at it when I get back in a week. I will gladly wrap your code mods into RMark if you want. Not sure why the ci's are different. You can ask Gary to see if he is doing anything differently but he doesn't use R so might not get much from your code.

regards --jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: var.components

Postby benaug » Wed Jul 25, 2012 1:08 pm

Thanks. I guess my question should have been why is the upper limit that value, not why is there an upper limit. Sure, I'll add some more code to get the AIC for the random effects model.
benaug
 
Posts: 6
Joined: Fri Jul 20, 2012 1:27 pm

Re: var.components

Postby jgoldberg » Fri Mar 22, 2019 5:08 pm

Hi all,

I think there's a small error in the code for calculating the random effects estimates (the betarand), at least when specifying a model with more than an intercept-only fixed effects design matrix. A small substitution fixes the issue easily.

Swapping
Code: Select all
 X = rep(1, length(theta))

for
Code: Select all
 X = design

works for me.

Cheers, Josh
jgoldberg
 
Posts: 1
Joined: Fri Mar 22, 2019 1:38 pm

Re: var.components

Postby jlaake » Tue Apr 16, 2019 3:08 pm

This error crept in with version 2.1.4 in 2012. I have made the suggested change and Josh tested for me. You can find a windows zip file for v2.2.7 at https://drive.google.com/open?id=1bbMw9DIcQhvo_N-ws_yNxKoLrAjBwXMM. For Mac or Linux users, the tar.gz file is at https://drive.google.com/open?id=12212j3zbWYcYyYxB0PIUQwYm85N2yTaM. Eventually I'll push to CRAN with possibly other changes.

Thanks to Josh for finding this error and the solution.

--jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to RMark

Who is online

Users browsing this forum: No registered users and 12 guests

cron