./blupsurv/R/bivrec.r

1. ZZdebug/fupdatefrailties3

ZZdebug/fupdatefrailties3 [ Functions ]

NAME

    fupdatefrailties3 --- update frailty estimates

FUNCTION

Drop-in replacement for fupdatefrailties4, which updates the frailties using only Fortran calls to low-level functionality like fsmuij. See fupdatefrailties4 for details.

SYNOPSIS

480 fupdatefrailties3 <- function(m, Ji, datamat, datamatd, alphars, alpharsd, as, asd, betahat, betadhat, dispparams)

SOURCE

483 {
484   #  cat("\tUpdating frailty estimates\n")
485         
486     ## Extract dispersion parameter estimates
487     sigma2hat <- dispparams$sigma2hat
488     sigma2dhat <- dispparams$sigma2dhat
489     nu2hat <- dispparams$nu2hat
490     nu2dhat <- dispparams$nu2dhat
491     thetahat <- dispparams$thetahat
492     
493     ## Initialize frailty vectors
494     Uihat <- rep(0, m);Vihat <- rep(0, m);
495     Uijhat <- rep(0, sum(Ji));Vijhat <- rep(0, sum(Ji));
496 
497    ## Initialize matrices for storage of pi, qi, pij, qij, etc.
498     jimax <- max(Ji)
499     pi <- rep(0, m);qi <- rep(0, m);ri <- rep(0, m);si <- rep(0, m)
500     piprime <- rep(0, m);qiprime <- rep(0, m);
501     riprime <- rep(0, m);siprime <- rep(0, m);
502     wi <- rep(0, m)
503     pij <- matrix(0, m, jimax);pijprime <- matrix(0, m, jimax);
504     qij <- matrix(0, m, jimax);qijprime <- matrix(0, m, jimax);
505     rij <- matrix(0, m, jimax);rijprime <- matrix(0, m, jimax);
506     sij <- matrix(0, m, jimax);sijprime <- matrix(0, m, jimax);
507     wij <- matrix(0, m, jimax);zij <- matrix(0, m, jimax);
508 
509     Smu <- matrix(0, m, jimax);Sdelta <- matrix(0, m, jimax);
510     Seta <- matrix(0, m, jimax);SDelta <- matrix(0, m, jimax);
511   
512     covs <- matrix(datamat[, -(1:8)], dim(datamat)[1], dim(datamat)[2] - 8)
513     ncovs <- dim(covs)[2]
514     d <- dim(covs)[1]
515     index <- datamat[, c("i", "j", "k", "r", "smin", "smax")]
516     delta <- datamat[, "delta"]
517     
518     Sm <- try(.Fortran("fsmuij", 
519         betahat = as.double(betahat),
520         index = as.integer(index),
521         delta = as.double(delta),
522         Z = as.double(covs),
523         alphars = as.double(alphars), 
524         as = as.double(as),
525         d = as.integer(d),
526         ncovs = as.integer(ncovs),
527         nr = as.integer(dim(alphars)[1]),
528         ns = as.integer(dim(alphars)[2]),
529         m = as.integer(m),
530         maxj = as.integer(max(Ji)),
531         Smu = as.double(Smu),
532         Sdelta = as.double(Sdelta)))
533     Smu <- matrix(Sm$Smu, m, jimax)
534     Sdelta <- matrix(Sm$Sdelta, m, jimax)
535     
536     covs <- matrix(datamatd[, -(1:8)], dim(datamatd)[1], dim(datamatd)[2] - 8)
537     ncovs <- dim(covs)[2]
538     d <- dim(covs)[1]
539     index <- datamatd[, c("i", "j", "k", "r", "smin", "smax")]
540     Delta <- datamatd[, "delta"]
541     
542     Se <- try(.Fortran("fsmuij",
543          betahat = as.double(betadhat),
544          index = as.integer(index),
545          delta = as.double(Delta),
546          Z = as.double(covs),
547          alphars = as.double(alpharsd),
548          as = as.double(asd),
549          d = as.integer(d),
550          ncovs = as.integer(ncovs),
551          nr = as.integer(dim(alpharsd)[1]),
552          ns = as.integer(dim(alpharsd)[2]),
553          m = as.integer(m),
554          maxj = as.integer(max(Ji)),
555          Smu = as.double(Seta),
556          Sdelta = as.double(SDelta)))
557     Seta <- matrix(Se$Smu, m, jimax)
558     SDelta <- matrix(Se$Sdelta, m, jimax)
559 
560     
561     for(i in 1:m){
562         for(j in 1:Ji[i]){
563             Smuij <- Smu[i, j]
564             Setaij <- Seta[i, j]
565             Sdeltaij <- Sdelta[i, j]
566             SDeltaij <- SDelta[i, j]
567             
568             wij[i, j] <- nu2hat - thetahat^2 * Setaij / (1 + nu2dhat * Setaij)
569             zij[i, j] <- nu2dhat - thetahat^2 * Smuij / (1 + nu2hat * Smuij)
570             pij[i, j] <- Smuij / (1 + wij[i, j] * Smuij)
571             qij[i, j] <- -thetahat * Smuij * Setaij / ((1 + nu2hat * Smuij) * 
572                 (1 + nu2dhat * Setaij) - thetahat^2 * Smuij * Setaij)
573             rij[i, j] <- qij[i, j]
574             sij[i, j] <- Setaij / (1 + zij[i, j] * Setaij)
575             pijprime[i, j] <- Sdeltaij / (1 + wij[i, j] * Smuij)
576             qijprime[i, j] <- -thetahat * Smuij * SDeltaij / 
577                 ((1 + nu2hat * Smuij) * (1 + nu2dhat * Setaij) - 
578                 thetahat^2 * Smuij * Setaij)
579             rijprime[i, j] <- -thetahat * Sdeltaij * Setaij / 
580                 ((1 + nu2hat * Smuij) * (1 + nu2dhat * Setaij) - 
581                 thetahat^2 * Smuij * Setaij)
582             sijprime[i, j] <- SDeltaij / (1 + zij[i, j] * Setaij)
583         }
584         ## Compute pi, qi, etc, as well as wi.
585         piprime[i] <- sum(pijprime[i, ])
586         qiprime[i] <- sum(qijprime[i, ])
587         riprime[i] <- sum(rijprime[i, ])
588         siprime[i] <- sum(sijprime[i, ])
589         
590         pi[i] <- sum(pij[i, ])
591         qi[i] <- sum(qij[i, ])
592         ri[i] <- sum(rij[i, ])
593         si[i] <- sum(sij[i, ])
594         
595         wi[i] <- 1 / ((1 + sigma2hat * pi[i]) * 
596             (1 + sigma2dhat * si[i]) - sigma2hat * sigma2dhat * ri[i]^2)
597     
598         ## Compute cluster frailty estimates
599         Uihat[i] <- 1 + sigma2hat * wi[i] * ((1 + sigma2dhat * si[i]) * 
600             (piprime[i] - pi[i] + qiprime[i] - qi[i]) - sigma2dhat * ri[i] * 
601             (riprime[i] - ri[i] + siprime[i] - si[i]))
602         Vihat[i] <- 1 + sigma2dhat * wi[i] * ((1 + sigma2hat * pi[i]) *
603             (riprime[i] - ri[i] + siprime[i] - si[i]) - sigma2hat * ri[i] *
604             (piprime[i] - pi[i] + qiprime[i] - qi[i]))
605         
606         ## Compute individual frailty estimates
607         Jicum <- c(0, cumsum(Ji))
608         for(j in 1:jimax){
609             Uijhat[Jicum[i] + j] <- Uihat[i] - (nu2hat * pij[i, j] + thetahat *
610                 rij[i, j]) * Uihat[i] - (nu2hat * qij[i, j] + thetahat * sij[i, j]) *
611                 Vihat[i] + nu2hat * (pijprime[i, j] + qijprime[i, j]) + thetahat *
612                 (rijprime[i, j] + sijprime[i, j])
613             Vijhat[Jicum[i] + j] <- Vihat[i] - (thetahat * qij[i, j] + nu2dhat *
614                 sij[i, j]) * Vihat[i] - (thetahat * pij[i, j] + nu2dhat * rij[i, j]) *
615                 Uihat[i] + thetahat * (pijprime[i, j] + qijprime[i, j]) + nu2dhat *
616                 (rijprime[i, j] + sijprime[i, j])
617         }      
618     }
619     
620     Uijhat[Uijhat<.01] <- 0.01; Vijhat[Vijhat < 0.01] <- 0.01
621     
622     Uijframe <- makeUijframe(m, Ji, Uijhat, Vijhat)
623     
624     cat("mean(Uijhat): ", mean(Uijhat), "   mean(Vijhat): ", mean(Vijhat), "\n")
625     if(mean(Uijhat)<.5 | mean(Vijhat)<.5) browser()
626 
627     return(list(Uihat = Uihat,
628                 Vihat = Vihat,
629                 Uijhat = Uijhat,
630                 Vijhat = Vijhat,
631                 pqrs = list(pij = pij, qij = qij, rij = rij, sij = sij,
632                     pijprime = pijprime, qijprime = qijprime,
633                     rijprime = rijprime, sijprime = sijprime,
634                     pi = pi, qi = qi, ri = ri, si = si,
635                     piprime = piprime, qiprime = qiprime,
636                     riprime = riprime, siprime = siprime,
637                     wi = wi, zij = zij, wij = wij),
638                 Uijframe = Uijframe))
639 }

TABLE OF CONTENTS

ZZdebug/fupdatefrailties3 [ Functions ]