TABLE OF CONTENTS
estimation/fupdatefrailties4 [ Functions ]
NAME
fupdatefrailties4 --- update frailty estimates
FUNCTION
Compute estimates of the cluster- and subject-level frailties. This function acts as a wrapper for the FORTRAN routine fmkfrail. See the fmkfrail documentation for details, or see fupdatefrailties3 for an R implementation.
SYNOPSIS
671 fupdatefrailties4 <- function(m, Ji, datamat, datamatd, alphars, alpharsd, 672 as, asd, betahat, betadhat, dispparams)
INPUTS
m number of clusters
Ji cluster sizes
datamat data matrix generated by makedatamat for event 1
datamatd data matrix generated by makedatamat for event 2
alphars matrix of baseline hazard parameters for event 1
alpharsd matrix of baseline hazard parameters for event 2
as matrix of discretization breakpoints for event 1
asd matrix of discretization breakpoints for event 2
betahat regression coefficient estimates for event 1
betadhat regression coefficient estimates for event 2
dispparams list of dispersion parameter estimates, see updatedisppears
OUTPUTS
Uihat vector of cluster-level frailties for event 1
Vihat vector of cluster-level frailties for event 2
Uijhat vector of subject-level frailties for event 1
Vijhat vector of subject-level frailties for event 2
pqrs list of intermediate values which are useful for computing
dispersion parameter estimates
Uijframe list containing matrix versions of Uijhat and Vijhat
SOURCE
675 { 676 # Extract dispersion parameters 677 sigma2hat <- dispparams$sigma2hat 678 sigma2dhat <- dispparams$sigma2dhat 679 nu2hat <- dispparams$nu2hat 680 nu2dhat <- dispparams$nu2dhat 681 thetahat <- dispparams$thetahat 682 683 ## Initialize frailty vectors 684 Uihat <- rep(0, m);Vihat <- rep(0, m); 685 Uijhat <- rep(0, sum(Ji));Vijhat <- rep(0, sum(Ji)); 686 687 # Initialize matrices for storage of pi, qi, pij, qij, etc. 688 jimax <- max(Ji) 689 pi <- rep(0, m);qi <- rep(0, m);ri <- rep(0, m);si <- rep(0, m) 690 piprime <- rep(0, m);qiprime <- rep(0, m); 691 riprime <- rep(0, m);siprime <- rep(0, m); 692 wi <- rep(0, m) 693 pij <- matrix(0, m, jimax);pijprime <- matrix(0, m, jimax); 694 qij <- matrix(0, m, jimax);qijprime <- matrix(0, m, jimax); 695 rij <- matrix(0, m, jimax);rijprime <- matrix(0, m, jimax); 696 sij <- matrix(0, m, jimax);sijprime <- matrix(0, m, jimax); 697 wij <- matrix(0, m, jimax);zij <- matrix(0, m, jimax); 698 699 # Prepare data for submission to Fortran function fmkfrail 700 Z <- matrix(datamat[, -(1:8)], dim(datamat)[1], dim(datamat)[2] - 8) 701 ncovs1 <- dim(Z)[2] 702 d1 <- dim(Z)[1] 703 index <- datamat[, c("i", "j", "k", "r", "smin", "smax")] 704 delta <- datamat[, "delta"] 705 Zd <- matrix(datamatd[, -(1:8)], dim(datamatd)[1], dim(datamatd)[2] - 8) 706 ncovs2 <- dim(Zd)[2] 707 d2 <- dim(Zd)[1] 708 indexd <- datamatd[, c("i", "j", "k", "r", "smin", "smax")] 709 deltad <- datamatd[, "delta"] 710 nr = dim(alphars)[1] 711 ns = dim(alphars)[2] 712 nsd = dim(alpharsd)[2] 713 714 # Submit to Fortran 715 out <- try(.Fortran("fmkfrail", 716 index = as.integer(index), # indices i,j,k,r,smin,smax 717 indexd = as.integer(indexd), 718 delta = as.double(delta), # event indicators 719 deltad = as.double(deltad), 720 Z = as.double(Z), # covariates 721 Zd = as.double(Zd), 722 alphars = as.double(alphars), # hazards 723 alpharsd = as.double(alpharsd), 724 as = as.double(as), # breakpoints 725 asd = as.double(asd), 726 betahat = as.double(betahat), # regression coefficients 727 betadhat = as.double(betadhat), 728 m = as.integer(m), # clusters and cluster sizes 729 Ji = as.integer(Ji), 730 jimax = as.integer(max(Ji)), 731 Jicum = as.integer(c(0, cumsum(Ji))), 732 jisum = as.integer(sum(Ji)), 733 d1 = as.integer(d1), # dimensions of covariate vectors 734 d2 = as.integer(d2), 735 ncovs1 = as.integer(ncovs1), 736 ncovs2 = as.integer(ncovs2), 737 nr = as.integer(nr), # number of strata 738 ns = as.integer(ns), # number of breakpoints 739 nsd = as.integer(nsd), 740 sigma2 = as.double(sigma2hat), # current dispersion parameters 741 sigma2d = as.double(sigma2dhat), 742 nu2 = as.double(nu2hat), 743 nu2d = as.double(nu2dhat), 744 theta = as.double(thetahat), 745 Uihat = as.double(Uihat), # emtpy matrices to store frailty estimates 746 Vihat = as.double(Vihat), 747 Uijhat = as.double(Uijhat), 748 Vijhat = as.double(Vijhat), 749 pi = as.double(pi), # matrices to store intermediate values 750 qi = as.double(qi), 751 ri = as.double(ri), 752 si = as.double(si), 753 piprime = as.double(piprime), 754 qiprime = as.double(qiprime), 755 riprime = as.double(riprime), 756 siprime = as.double(siprime), 757 pij = as.double(pij), 758 qij = as.double(qij), 759 rij = as.double(rij), 760 sij = as.double(sij), 761 pijprime = as.double(pijprime), 762 qijprime = as.double(qijprime), 763 rijprime = as.double(rijprime), 764 sijprime = as.double(sijprime), 765 wi = as.double(wi), 766 wij = as.double(wij), 767 zij = as.double(zij))) 768 769 # Too small frailty estimates are not allowed 770 out$Uijhat[out$Uijhat<.01] <- 0.01; out$Vijhat[out$Vijhat < 0.01] <- 0.01 771 772 # Extract output 773 pij = matrix(out$pij, m, jimax) 774 qij = matrix(out$qij, m, jimax) 775 rij = matrix(out$rij, m, jimax) 776 sij = matrix(out$sij, m, jimax) 777 778 pijprime = matrix(out$pijprime, m, jimax) 779 qijprime = matrix(out$qijprime, m, jimax) 780 rijprime = matrix(out$rijprime, m, jimax) 781 sijprime = matrix(out$sijprime, m, jimax) 782 783 zij = matrix(out$zij, m, jimax) 784 wij = matrix(out$wij, m, jimax) 785 786 787 # Format frailtyoutput list 788 return(list(Uihat = out$Uihat, 789 Vihat = out$Vihat, 790 Uijhat = out$Uijhat, 791 Vijhat = out$Vijhat, 792 pqrs = list(pij = pij, qij = qij, rij = rij, sij = sij, 793 pijprime = pijprime, qijprime = qijprime, 794 rijprime = rijprime, sijprime = sijprime, 795 pi = out$pi, qi = out$qi, ri = out$ri, si = out$si, 796 piprime = out$piprime, qiprime = out$qiprime, 797 riprime = out$riprime, siprime = out$siprime, 798 wi = out$wi, zij = zij, wij = wij), 799 Uijframe = makeUijframe(m, Ji, out$Uijhat, out$Vijhat))) 800 }