TABLE OF CONTENTS


ZZdebug/makeSens2 [ Functions ]

NAME

    makeSens2 --- construct sensitivity matrix

FUNCTION

Construct the sensitivity matrix at the end of the procedure. This is an R implementation of the fmksens2 FORTRAN routine used for debugging.

SYNOPSIS

1976 makeSens2 <- function(m, Ji, datamat, datamatd, alphars, alpharsd, betahat, betadhat,
1977             as, asd, frailtyoutput, dispersionoutput, K, Kd, full = FALSE)

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
    betahat            regression coefficient estimates for event 
    betadhat           regression coefficient estimates for event 2
    as                 matrix of discretization breakpoints for event 1
    asd                matrix of discretization breakpoints for event 2
    frailtyoutput  output from fupdatefrailties4
    dispersionoutput  output from one of the dispersion estimation routines
    K                  number of discretization intervals for event 1
    Kd                 number of discretization intervals for event 2
    full           boolean indicating whether full matrix should be computed.
                   If FALSE, only the sensitivity for betahat is returned. 

OUTPUTS

    S              Sensitivity matrix

SOURCE

1980 {
1981     # Extract intermediate values from output
1982     pi <- frailtyoutput$pqrs$pi; qi <- frailtyoutput$pqrs$qi;
1983     ri <- frailtyoutput$pqrs$ri; si <- frailtyoutput$pqrs$si;
1984     wij <- frailtyoutput$pqrs$wij;zij <- frailtyoutput$pqrs$zij;
1985     wi <- frailtyoutput$pqrs$wi
1986     sigma2hat <- dispersionoutput$sigma2hat;
1987     sigma2dhat <- dispersionoutput$sigma2dhat
1988     nu2hat <- dispersionoutput$nu2hat; nu2dhat <- dispersionoutput$nu2dhat
1989     thetahat <- dispersionoutput$thetahat
1990    
1991     ncovs <- length(betahat)
1992 
1993     Kcum <- cumsum(c(1, K))
1994     Kdcum <- cumsum(c(1, Kd))
1995     if(full){
1996         betaind <- sum(K) + 1:ncovs
1997         betadind <- sum(Kd) + 1:ncovs
1998     }else{
1999         betaind <- 1:ncovs
2000         betadind <- 1:ncovs
2001     }
2002     
2003     # Allocate storage
2004     if(full){
2005     S <- matrix(0, sum(K) + sum(Kd) + 2 * ncovs, sum(K) + sum(Kd) + 2 * ncovs)
2006     }else{
2007     S <- matrix(0, 2 * ncovs, 2 * ncovs)
2008     }
2009     
2010     # The sensitivity matrix is expressed as the sum of matrices corresponding
2011     # to each cluster
2012     for(i in 1:m)
2013     {    
2014         # Allocate storage for intermediate matrices
2015         Smuij <- rep(0, Ji[i]); Setaij <- rep(0, Ji[i])
2016         if(full){
2017             Smuxij <- matrix(0, Ji[i], sum(K) + ncovs);
2018             Setaxijd <- matrix(0, Ji[i], sum(Kd) + ncovs)
2019             S11 <- matrix(0, sum(K) + ncovs, sum(K) + ncovs);
2020             S12 <- matrix(0, sum(K) + ncovs, sum(Kd) + ncovs);
2021             S13 <- matrix(0, sum(Kd) + ncovs, sum(Kd) + ncovs);
2022             S21 <- rep(0, sum(K) + ncovs); S22 <- rep(0, sum(Kd) + ncovs);
2023             S31 <- rep(0, sum(K) + ncovs); S32 <- rep(0, sum(Kd) + ncovs);
2024         }else{
2025             Smuxij <- matrix(0, Ji[i], ncovs);
2026             Setaxijd <- matrix(0, Ji[i], ncovs)
2027             S11 <- matrix(0, ncovs, ncovs);
2028             S12 <- matrix(0, ncovs, ncovs);
2029             S13 <- matrix(0, ncovs, ncovs);
2030             S21 <- rep(0, ncovs); S22 <- rep(0, ncovs);
2031             S31 <- rep(0, ncovs); S32 <- rep(0, ncovs);
2032         }
2033                 
2034         # Find the index of the data matrices where cluster i begins
2035         iind0 <- match(i, datamat[, "i"])
2036         if(i < m) iind1 <- match(i + 1, datamat[, "i"]) else 
2037             iind1 <- dim(datamat)[1] + 1
2038         iind0d <- match(i, datamatd[, "i"])
2039         if(i < m) iind1d <- match(i + 1, datamatd[, "i"]) else 
2040             iind1d <- dim(datamatd)[1] + 1
2041         
2042         # Compute Smuij, Setaij, Smuxij, Setaxijd (first pass)
2043         for(ind in iind0:(iind1 - 1))
2044         {
2045             Z <- matrix(datamat[ind, -(1:8)], ncovs, 1)
2046             bZ <- t(Z)%*%as.matrix(betahat)
2047             r <- datamat[ind, "r"]
2048             smax <- datamat[ind, "smax"]
2049             k <- datamat[ind, "k"]
2050             j <- datamat[ind, "j"]
2051             time <- datamat[ind, "time"]
2052              if(smax > 0){
2053                 for(s in 1:smax){
2054                     mu <- exp(bZ) * alphars[r, s] * A(time, as, r, s)
2055                     Smuij[j] <- Smuij[j] + mu
2056                     Smuxij[j, betaind] <- Smuxij[j, betaind] + mu * as.vector(Z)    
2057                     if(full) Smuxij[j, Kcum[r] + s - 1] <- 
2058                         Smuxij[j, Kcum[r] + s - 1] + mu        
2059                 }
2060             }
2061         }
2062         for(ind in iind0d:(iind1d - 1))
2063         {     
2064             Z <- matrix(datamatd[ind, -(1:8)], ncovs, 1)
2065             bdZ <- t(Z)%*%as.matrix(betadhat)
2066             r <- datamatd[ind, "r"]
2067             smaxd <- datamatd[ind, "smax"]
2068             k <- datamatd[ind, "k"]
2069             j <- datamatd[ind, "j"]
2070             time <- datamatd[ind, "time"]
2071              if(smaxd > 0){
2072                 for(s in 1:smaxd){
2073                     eta <- exp(bdZ) * alpharsd[r, s] * A(time, asd, r, s)
2074                     Setaij[j] <- Setaij[j] + eta
2075                     Setaxijd[j, betadind] <- Setaxijd[j, betadind] + eta * as.vector(Z)
2076                     if(full) Setaxijd[j, Kdcum[r] + s - 1] <- 
2077                         Setaxijd[j, Kdcum[r] + s - 1] + eta        
2078                 }
2079             }
2080         }  
2081         
2082         phi <- rep(0, Ji[i])
2083         for(j in 1:Ji[i]){
2084             phi[j] <- (1 + nu2hat * Smuij[j]) * (1 + nu2dhat * Setaij[j]) - 
2085                 thetahat^2 * Smuij[j] * Setaij[j]
2086         }         
2087         
2088         
2089         ## Compute the various components of the matrix (second pass) 
2090         for(ind in iind0:(iind1 - 1))
2091         {
2092             Z <- matrix(datamat[ind, -(1:8)], ncovs, 1)
2093             bZ <- t(Z)%*%as.matrix(betahat)
2094             Z2 <- Z%*%t(Z)
2095             Z <- as.vector(Z)
2096             r <- datamat[ind, "r"]
2097             smax <- datamat[ind, "smax"]
2098             k <- datamat[ind, "k"]     
2099             j <- datamat[ind, "j"]
2100             time <- datamat[ind, "time"]
2101             if(smax > 0){
2102                 for(s in 1:smax){   
2103                     mu <- exp(bZ) * alphars[r, s] * A(time, as, r, s)
2104                     S11[betaind, betaind] <- S11[betaind, betaind] + 
2105                         as.vector(mu) * Z2
2106                     if(full){
2107                         S11[Kcum[r] + s - 1, Kcum[r] + s - 1] <- 
2108                             S11[Kcum[r] + s - 1, Kcum[r] + s - 1] + mu
2109                         S11[Kcum[r] + s - 1, betaind] <- 
2110                             S11[Kcum[r] + s - 1, betaind] + mu * Z
2111                         S11[betaind, Kcum[r] + s - 1] <- 
2112                             S11[betaind, Kcum[r] + s - 1] + mu * Z
2113                     }
2114                     w <- 1 / (1 + wij[i, j] * Smuij[j])
2115                     S21[betaind] <- S21[betaind] + w * mu * Z
2116                     if(full) S21[Kcum[r] + s - 1] <- S21[Kcum[r] + s - 1] + w * mu
2117                     
2118                     w<- -thetahat * Setaij[j] / phi[j]
2119                     S31[betaind] <- S31[betaind] + w * mu * Z 
2120                     if(full) S31[Kcum[r] + s - 1] <- S31[Kcum[r] + s - 1] + w * mu
2121                 }
2122             }
2123         }
2124         for(ind in iind0d:(iind1d - 1))
2125         {
2126            Z <- matrix(datamatd[ind, -(1:8)], ncovs, 1)
2127            bdZ <- t(Z)%*%as.matrix(betadhat)
2128            Z2 <- Z%*%t(Z)
2129            Z <- as.vector(Z)
2130            r <- datamatd[ind, "r"]
2131            smaxd <- datamatd[ind, "smax"]
2132            k <- datamatd[ind, "k"]     
2133            j <- datamatd[ind, "j"]
2134            time <- datamatd[ind, "time"]
2135            if(smaxd > 0){
2136                 for(s in 1:smaxd){   
2137                     eta <- exp(bdZ) * alpharsd[r, s] * A(time, asd, r, s)
2138                     S13[betadind, betadind] <- S13[betadind, betadind] +
2139                         as.vector(eta) * Z2
2140                     if(full){
2141                         S13[Kdcum[r] + s - 1, Kdcum[r] + s - 1] <- 
2142                             S13[Kdcum[r] + s - 1, Kdcum[r] + s - 1] + eta
2143                         S13[Kdcum[r] + s - 1, betadind] <-
2144                             S13[Kdcum[r] + s - 1, betadind] + eta * Z
2145                         S13[betadind, Kdcum[r] + s - 1] <-
2146                             S13[betadind, Kdcum[r] + s - 1] + eta * Z
2147                     }
2148                     w<- -thetahat * Smuij[j] / phi[j]
2149                     S22[betadind] <- S22[betadind] + w * eta * Z
2150                     if(full) S22[Kdcum[r] + s - 1] <- S22[Kdcum[r] + s - 1] + w * eta
2151                     
2152                     w <- 1 / (1 + zij[i, j] * Setaij[j])
2153                     S32[betadind] <- S32[betadind] + w * eta * Z
2154                     if(full) S32[Kdcum[r] + s - 1] <- S32[Kdcum[r] + s - 1] + w * eta
2155                 }
2156             }
2157         }
2158         
2159         for(j in 1:Ji[i]){
2160             S11 <- S11 - wij[i, j] / (1 + wij[i, j] * Smuij[j]) * 
2161                 Smuxij[j, ]%*%t(Smuxij[j, ])
2162             S12 <- S12 - thetahat / phi[j] * Smuxij[j, ]%*%t(Setaxijd[j, ])
2163             S13 <- S13 - zij[i, j] / (1 + zij[i, j] * Setaij[j]) * 
2164                 Setaxijd[j, ]%*%t(Setaxijd[j, ])
2165         }
2166         
2167         S2 <- c(S21, S22)
2168         S3 <- c(S31, S32)
2169                 
2170         S <- S - rbind(cbind(S11, S12), cbind(t(S12), S13)) +
2171             wi[i] * (sigma2hat * (1 + sigma2dhat * si[i]) * S2%*%t(S2) -
2172             sigma2hat * sigma2dhat * qi[i] * (S3%*%t(S2) + S2%*%t(S3)) + 
2173             sigma2dhat * (1 + sigma2hat * pi[i]) * S3%*%t(S3))
2174     }
2175     return(S)
2176 }