TABLE OF CONTENTS
estimation/updatedispohls [ Functions ]
NAME
updatedispohls --- Ohlsson-type dispersion parameter estimators
FUNCTION
Compute dispersion parameter estimators using the method proposed by Ohlsson et al. See the paper for the reference.
SYNOPSIS
956 updatedispohls <- function(m, Ji, datamat, datamatd, alphars, alpharsd, 957 as, asd, betahat, betadhat, fixzero)
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
fixzero list of estimators that should be fixed at zero (testing only)
OUTPUTS
sigma2hat cluster-level dispersion for process 1
sigma2dhat cluster-level dispersion for process 2
nu2hat subject-level dispersion for process 1
nu2dhat subject-level dispersion for process 2
thetahat subject-level covariance
SOURCE
960 { 961 jimax <- max(Ji) 962 Smu <- matrix(0, m, jimax);Sdelta <- matrix(0, m, jimax); 963 Seta <- matrix(0, m, jimax);SDelta <- matrix(0, m, jimax); 964 965 # Prepare data for Fortran function fsmuij to compute sum for event 1 966 covs <- matrix(datamat[, -(1:8)], dim(datamat)[1], dim(datamat)[2] - 8) 967 ncovs <- dim(covs)[2] 968 d <- dim(covs)[1] 969 index <- datamat[, c("i", "j", "k", "r", "smin", "smax")] 970 delta <- datamat[, "delta"] 971 972 # The FORTRAN function fsmuij computes a matrix of size m x max(Ji), whose 973 # (i, j) entry is sum(mu_rijks) over all r,k,s 974 Sm <- try(.Fortran("fsmuij", 975 betahat = as.double(betahat), 976 index = as.integer(index), 977 delta = as.double(delta), 978 Z = as.double(covs), 979 alphars = as.double(alphars), 980 as = as.double(as), 981 d = as.integer(d), 982 ncovs = as.integer(ncovs), 983 nr = as.integer(dim(alphars)[1]), 984 ns = as.integer(dim(alphars)[2]), 985 m = as.integer(m), maxj = as.integer(max(Ji)), 986 Smu = as.double(Smu), 987 Sdelta = as.double(Sdelta))) 988 989 # extract sums from FORTRAN output 990 Smu <- matrix(Sm$Smu, m, jimax) 991 Sdelta <- matrix(Sm$Sdelta, m, jimax) 992 993 # Repeat FORTRAN call for event 2 994 covs <- matrix(datamatd[, -(1:8)], dim(datamatd)[1], dim(datamatd)[2] - 8) 995 ncovs <- dim(covs)[2] 996 d <- dim(covs)[1] 997 index <- datamatd[, c("i", "j", "k", "r", "smin", "smax")] 998 Delta <- datamatd[, "delta"] 999 1000 # Compute sum of eta_ij analogously 1001 Se <- try(.Fortran("fsmuij", 1002 betahat = as.double(betadhat), 1003 index = as.integer(index), 1004 delta = as.double(Delta), 1005 Z = as.double(covs), 1006 alphars = as.double(alpharsd), 1007 as = as.double(asd), 1008 d = as.integer(d), 1009 ncovs = as.integer(ncovs), 1010 nr = as.integer(dim(alpharsd)[1]), 1011 ns = as.integer(dim(alpharsd)[2]), 1012 m = as.integer(m), 1013 maxj = as.integer(max(Ji)), 1014 Smu = as.double(Seta), 1015 Sdelta = as.double(SDelta))) 1016 1017 Seta <- matrix(Se$Smu, m, jimax) 1018 SDelta <- matrix(Se$Sdelta, m, jimax) 1019 1020 # Compute sum(mu_i*) for both processes 1021 Smui <- apply(Smu, 1, sum) 1022 Setai <- apply(Seta, 1, sum) 1023 1024 # Begin computing the Ohlsson-type estimates 1025 Xtild <- Sdelta / Smu 1026 Ztild <- SDelta / Seta 1027 Xtildi <- apply(Sdelta, 1, sum) / Smui 1028 Ztildi <- apply(SDelta, 1, sum) / Setai 1029 Ztildii <- rep(Ztildi, Ji) 1030 1031 # This is a simplistic correction to attempt to resolve some occasional 1032 # computational problems caused by too large or too small estimated means. 1033 # This line replaces subject estimates which are in the most extreme 5\% by 1034 # the cluster-level estimate for the corresponding cluster. 1035 Ztild[Ztild > quantile(Ztild, .975) | Ztild < quantile(Ztild, .025)] <- 1036 Ztildii[Ztild > quantile(Ztild, .975) | Ztild < quantile(Ztild, .025)] 1037 1038 # Subject-level Ohlsson estimates 1039 nu2hat <- 0;nu2dhat <- 0;thetahat <- 0;thetadenom <- 0;sigma2hat <- 0;sigma2dhat <- 0; 1040 for(i in 1:m){ 1041 for(j in 1:Ji[i]){ 1042 nu2hat <- nu2hat + Smu[i, j] * (Xtild[i, j] - Xtildi[i])^2 1043 nu2dhat <- nu2dhat + Seta[i, j] * (Ztild[i, j] - Ztildi[i])^2 1044 thetahat <- thetahat + (Xtild[i, j] - Xtildi[i]) * 1045 (Ztild[i, j] - Ztildi[i]) 1046 thetadenom <- thetadenom + (1 + 1 / Smui[i] / Setai[i] * 1047 sum(Smu[i, ] * Seta[i, ])) - (Smu[i, j] / Smui[i] + 1048 Seta[i, j] / Setai[i]) 1049 } 1050 } 1051 # Complete the estimates of the dispersion parameters. 1052 nu2hat <- (nu2hat - sum(Ji - 1)) / 1053 (sum(Smui) - sum(apply(Smu^2, 1, sum) / Smui)) 1054 nu2dhat <- (nu2dhat - sum(Ji - 1)) / 1055 (sum(Setai) - sum(apply(Seta^2, 1, sum) / Setai)) 1056 thetahat <- thetahat / thetadenom 1057 1058 # Fix parameters at 0 if desired 1059 if(!is.null(fixzero)){ 1060 if("nu2hat"%in%fixzero) nu2hat <- 0 1061 if("nu2dhat"%in%fixzero) nu2dhat <- 0 1062 if("thetahat"%in%fixzero) thetahat <- 0 1063 } 1064 1065 # The Ohlsson - type estimators do not guarantee that the esimated covariance 1066 # matrix is valid. This code truncates the estimate of thetahat if needed 1067 if(nu2hat<.001) nu2hat <- .001 1068 if(nu2dhat<.001) nu2dhat <- .001 1069 badtheta <- try(thetahat / sqrt(nu2hat * nu2dhat) > 1 | 1070 thetahat / sqrt(nu2hat * nu2dhat)< -1) 1071 if(inherits(badtheta, "try - error") | is.na(badtheta)) {browser(); return(0)} 1072 if(badtheta) thetahat <- try(sign(thetahat) * sqrt(nu2hat * nu2dhat)) 1073 if(inherits(thetahat, "try - error")) return(0) 1074 1075 #! Estimate the cluster - level dispersion parameters 1076 zij <- Smu / (Smu + 1 / nu2hat) 1077 zi <- apply(zij, 1, sum) 1078 zijd <- Seta / (Seta + 1 / nu2dhat) 1079 zid <- apply(zijd, 1, sum) 1080 1081 Xztildi <- apply(zij * Xtild, 1, sum) / zi 1082 Zztildi <- apply(zijd * Ztild, 1, sum) / zid 1083 1084 Xztild <- sum(zi * Xztildi) / sum(zi) 1085 Zztild <- sum(zid * Zztildi) / sum(zid) 1086 1087 sigma2hat <- (sum(zi * (Xztildi - Xztild)^2) - nu2hat * (m - 1)) / 1088 (sum(zi) - sum(zi^2) / sum(zi)) 1089 sigma2dhat <- (sum(zid * (Zztildi - Zztild)^2) - nu2dhat * (m - 1)) / 1090 (sum(zid) - sum(zid^2) / sum(zid)) 1091 1092 # Fix parameters at 0 if desired 1093 if(!is.null(fixzero)){ 1094 if("sigma2hat"%in%fixzero) sigma2hat <- 0 1095 if("sigma2dhat"%in%fixzero) sigma2dhat <- 0 1096 } 1097 1098 # Format for output 1099 return(list(sigma2hat = sigma2hat, 1100 sigma2dhat = sigma2dhat, 1101 nu2hat = nu2hat, 1102 nu2dhat = nu2dhat, 1103 thetahat = thetahat)) 1104 }