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 }