TABLE OF CONTENTS
estimation/getdispmarg3 [ Functions ]
NAME
getdispmarg3 --- workhorse function for computing marginal estimators
FUNCTION
This function is used by updatedispmarg to compute the sums of cross products of event indicators and means. Inputs are data for any two processes, which may be event processes 1 and 2, or two copies of event 1, or two copies of event 2. The sample variances and covariances of the two processes are computed under the auxiliary Poisson model.
SYNOPSIS
1209 getdispmarg3 <- function(m, Ji, datamat1, datamat2, alphars1, alphars2, 1210 as1, as2, betahat1, betahat2, same)
INPUTS
datamat1 data matrix generated by makedatamat for process 1 datamat2 data matrix generated by makedatamat for process 2 alphars1 matrix of baseline hazard parameters for process 1 alphars2 matrix of baseline hazard parameters for process 2 as1 matrix of discretization breakpoints for process 1 as2 matrix of discretization breakpoints for process 2 betahat1 regression coefficient estimates for process 1 betadha2 regression coefficient estimates for process 2 same boolean indicator of whether process 1 and 2 are the same
OUTPUTS
sig2nu2 marginal estimate of sigma^2 + nu^2 (or theta, if the processes
are different)
sig2 marginal estimate of sigma^2
SOURCE
1213 { 1214 1215 #Initialze vectors and matrices 1216 jimax <- max(Ji) 1217 Smu1 <- matrix(0, m, jimax);Sdelta1 <- matrix(0, m, jimax); 1218 Smu2 <- matrix(0, m, jimax);Sdelta2 <- matrix(0, m, jimax); 1219 1220 # Prepare for FORTRAN 1221 covs1 <- matrix(datamat1[, -(1:8)], dim(datamat1)[1], dim(datamat1)[2] - 8) 1222 ncovs1 <- dim(covs1)[2] 1223 d1 <- dim(covs1)[1] 1224 index1 <- datamat1[, c("i", "j", "k", "r", "smin", "smax")] 1225 delta1 <- datamat1[, "delta"] 1226 1227 # First, compute the sums for the first data set. 1228 Sm1 <- try(.Fortran("fsmuij", 1229 betahat = as.double(betahat1), 1230 index = as.integer(index1), 1231 delta = as.double(delta1), 1232 Z = as.double(covs1), 1233 alphars = as.double(alphars1), 1234 as = as.double(as1), 1235 d = as.integer(d1), 1236 ncovs = as.integer(ncovs1), 1237 nr = as.integer(dim(alphars1)[1]), 1238 ns = as.integer(dim(alphars1)[2]), 1239 m = as.integer(m), 1240 maxj = as.integer(max(Ji)), 1241 Smu = as.double(Smu1), 1242 Sdelta = as.double(Sdelta1))) 1243 1244 Smu1 <- matrix(Sm1$Smu, m, jimax) 1245 Sdelta1 <- matrix(Sm1$Sdelta, m, jimax) 1246 1247 # Compute sums for the second data set 1248 covs2 <- matrix(datamat2[, -(1:8)], dim(datamat2)[1], dim(datamat2)[2] - 8) 1249 ncovs2 <- dim(covs2)[2] 1250 d2 <- dim(covs2)[1] 1251 index2 <- datamat2[, c("i", "j", "k", "r", "smin", "smax")] 1252 delta2 <- datamat2[, "delta"] 1253 1254 Sm2 <- try(.Fortran("fsmuij", 1255 betahat = as.double(betahat2), 1256 index = as.integer(index2), 1257 delta = as.double(delta2), 1258 Z = as.double(covs2), 1259 alphars = as.double(alphars2), 1260 as = as.double(as2), 1261 d = as.integer(d2), 1262 ncovs = as.integer(ncovs2), 1263 nr = as.integer(dim(alphars2)[1]), 1264 ns = as.integer(dim(alphars2)[2]), 1265 m = as.integer(m), 1266 maxj = as.integer(max(Ji)), 1267 Smu = as.double(Smu2), 1268 Sdelta = as.double(Sdelta2))) 1269 1270 Smu2 <- matrix(Sm2$Smu, m, jimax) 1271 Sdelta2 <- matrix(Sm2$Sdelta, m, jimax) 1272 1273 # This computation relies on the fact that for vectors x_1, x_2, 1274 # a^Tx_1x_2^Te = (a^Tx_1)(a^Tx_2). Therefore, 1275 # the numerator and denominator of sig2nu2= can be computed as 1276 sig2nu2num <- sum((Sdelta1 - Smu1) * (Sdelta2 - Smu2)) 1277 sig2nu2den <- sum(Smu1 * Smu2) 1278 #! If they come from the same data set, we must subtract the diagonal elements 1279 if(same){ 1280 Smud2 <- try(.Fortran("fsmud2", 1281 betahat = as.double(betahat1), 1282 index = as.integer(index1), 1283 delta = as.double(delta1), 1284 Z = as.double(covs1), 1285 alphars = as.double(alphars1), 1286 as = as.double(as1), 1287 d = as.integer(d1), 1288 ncovs = as.integer(ncovs1), 1289 nr = as.integer(dim(alphars1)[1]), 1290 ns = as.integer(dim(alphars1)[2]), 1291 Smud = as.double(0), 1292 Smu2 = as.double(0))) 1293 Smud <- Smud2$Smud 1294 Smusq <- Smud2$Smu2 1295 # Correct the estimates by removing the diagonal terms 1296 sig2nu2num <- sig2nu2num - Smud 1297 sig2nu2den <- sig2nu2den - Smusq 1298 } 1299 sig2nu2 <- sig2nu2num / sig2nu2den 1300 1301 # Computing the value of sig2 is analogous. 1302 # First, sum over all the subjects in each cluster, then compute the 1303 # numerator and denominator separately, by adding all the cross terms 1304 # and subtracting those for which b = j. 1305 Sdeltai1 <- apply(Sdelta1, 1, sum) 1306 Sdeltai2 <- apply(Sdelta2, 1, sum) 1307 Smui1 <- apply(Smu1, 1, sum) 1308 Smui2 <- apply(Smu2, 1, sum) 1309 sig2 <- (sum((Sdeltai1 - Smui1) * (Sdeltai2 - Smui2)) - sum((Sdelta1 - Smu1) * 1310 (Sdelta2 - Smu2))) / (sum(Smui1 * Smui2) - sum(Smu1 * Smu2)) 1311 1312 # Format for output 1313 return(list(sig2nu2 = sig2nu2, sig2 = sig2)) 1314 }