######################################################################################## ######## R code for "Generalized accelerated recurrence time model ###### ######## for multivariate recurrent event data with missing event type" ###### ######## For academic research only ###### ######################################################################################## library(quantreg) ######################################################################################## ############## First Part: R functions ################################# ######################################################################################## ######################################################################################## #### The kernel function K() used in the paper ################################## ######################################################################################## K <- function(x) { ### the Epanechnikov kernal function # 0.75*(1-x^2)*(abs(x)<1) ### the Gaussian kernel function (2*pi)^(-1/2)*exp(-1/2*x^2) } ######################################################################################## ######################################################################################## #### The function calpi() to calculate pi_i^j in the paper ###################### ######################################################################################## calpi <- function(A, X, T, h, ZETA) { #### N: a scalar denoted number of observations #### A: a N*1 vector denoted A_i^j #### X: a N*p matrix denoted by X_i #### T: a N*1 vector denoted T_i^j #### h: a scalar denoted the bandwidth #### ZETA: a N*1 vector denoted the perturbed variable zeta_i ZETA2 = matrix(rep(ZETA, times=N), ncol=N) z1z = matrix(rep(X[, 2], times=N), ncol=N) z2z = matrix(rep(X[, 3], times=N), ncol=N) T2 = matrix(rep(T, times=N), ncol=N) A2 = matrix(rep(A, times=N), ncol=N) tmp1 = apply( t(ZETA2)*(t(z1z)==z1z)*K((t(T2)-T2)/h)/h*K((t(z2z)-z2z)/h)/h*t(A2), 1, sum) tmp2 = apply( t(ZETA2)*(t(z1z)==z1z)*K((t(T2)-T2)/h)/h*K((t(z2z)-z2z)/h)/h, 1, sum) return(tmp1/tmp2) } ######################################################################################## ######################################################################################## #### The function calpp() to calculate p_i1^j and p_i2^j ###################### ######################################################################################## calpp <- function(A, X, T, delta1, delta2, h, ZETA) { #### N: a scalar denoted number of observations #### A: a N*1 vector denoted A_i^j #### X: a N*p matrix denoted by X_i #### T: a N*1 vector denoted T_i^j #### delta1: a N*1 vector denoted delta_i1^j #### delta2: a N*1 vector denoted delta_i2^j #### h: a scalar denoted the bandwidth #### ZETA: a N*1 vector denoted the perturbed variable zeta_i ZETA2 = matrix(rep(ZETA, times=N), ncol=N) delta1_2 = matrix(rep(delta1, times=N), ncol=N) delta2_2 = matrix(rep(delta2, times=N), ncol=N) z1z = matrix(rep(X[, 2], times=N), ncol=N) z2z = matrix(rep(X[, 3], times=N), ncol=N) T2 = matrix(rep(T, times=N), ncol=N) A2 = matrix(rep(A, times=N), ncol=N) tmp1 = apply( t(ZETA2)*(t(z1z)==z1z)*K((t(T2)-T2)/h)/h*K((t(z2z)-z2z)/h)/h*t(A2), 1, sum) tmp1_1 = apply( t(ZETA2)*(t(z1z)==z1z)*K((t(T2)-T2)/h)/h*K((t(z2z)-z2z)/h)/h*t(A2)*t(delta1_2), 1, sum) tmp1_2 = apply( t(ZETA2)*(t(z1z)==z1z)*K((t(T2)-T2)/h)/h*K((t(z2z)-z2z)/h)/h*t(A2)*t(delta2_2), 1, sum) return( cbind(tmp1_1/tmp1, tmp1_2/tmp1) ) } ######################################################################################## ######################################################################################## #### The function estf() to find estimators through minimization ################ ######################################################################################## estf <- function(wt, X, T, grids, link, U, x, L, R, zeta, ZETA) { #### wt: weight function, can be IPW, EEP, CC, and Full weights #### N: a scalar denoted number of observations #### X: a N*p matrix denoted by X_i #### T: a N*1 vector denoted T_i^j #### grids: number of equally grids in (0, U] #### link: the scale function g(s) #### U: the upper bound of time interval (0, U] #### x: a n by p matrix denoted the covariate #### L: a n by 1 vector denoted the left censoring time #### R: a n by 1 vector denoted the right censoring time #### zeta: a n by 1 vector denoted the perturbed variable zeta #### ZETA: a N*1 vector denoted the perturbed variable zeta_i rqy = matrix(NA, N+2, 1) rqx = matrix(NA, N+2, p) rqw = matrix(NA, N+2, 1) rqy[1:N] = log(T) rqy[N+1] = 10^10 rqy[N+2] = 10^10 rqw[1:N] = wt*ZETA rqw[(N+1):(N+2)] = 1 rqx[1:N, ] = X rqx[N+1, ] = -apply(X*matrix(rep(rqw[1:N, ], times=p), ncol=p), 2, sum) betaest = matrix(NA, grids, p) lc = matrix(0, n, 1) for (k in 1:grids) { if (k == 1) { lc = matrix(0, n, 1) if (link == "1") { lc[which(L==0)] = U/grids } else if (link == "u") { lc[which(L==0)] = (U/grids)^2/2 } else if (link == "exp") { lc[which(L==0)] = exp(U/grids)-1 } } if (k > 1) { tmp = which(L <= exp(x%*%betaest[k - 1, ]) & exp(x%*%betaest[k - 1, ])<= R) if (link == "1") { lc[tmp] = lc[tmp] + U/grids } else if (link == "u") { lc[tmp] = lc[tmp] + (U/grids*k)^2/2 - (U/grids*(k-1))^2/2 } else if (link == "exp") { lc[tmp] = lc[tmp] + exp(U/grids * k) - exp(U/grids*(k - 1)) } } rqx[N + 2, ] = 2*apply(x*matrix(rep(lc*zeta, times=p), ncol=p), 2, sum) est <- rq(rqy ~ rqx - 1, tau = 0.5, weights = rqw)$coefficients if (max(abs(est)) < 10^5) { betaest[k, ] = est } } return(betaest) } ######################################################################################## ######################################################################################## #### The function BJf() to calculate matrixes B and J at a given time point ###### ######################################################################################## BJf <- function(wt, bb, X, T, x, L, R, ID, id) { #### B and J are p*p matrixes #### wt: weight function, can be IPW, EEP, CC, and Full weights #### bb: a p*1 vector, equals to \hat{beta}_k(u) with u ranges from 1 to grids #### N: a scalar denoted number of observations #### X: a N*p matrix denoted by X_i #### T: a N*1 vector denoted T_i^j #### x: a n by p matrix denoted the covariate #### L: a n by 1 vector denoted the left censoring time #### R: a n by 1 vector denoted the right censoring time #### ID: a N*1 vector listed the ID of recurrent events #### id: a n*1 vector listed subject id tmp = wt*(T<=exp(X%*%bb)+10^-10) ln = matrix(NA, n, p) for (i in 1:n) { ln[i, ] = x[i, ]*sum( tmp[ which(ID==id[i]) ] ) } Ln = n^{-1/2}*apply(ln, 2, sum) Omega = t(ln)%*%ln/n e = eigen(Omega) E = (e$vectors)%*%diag(sqrt(e$values))%*%t(e$vectors) D = matrix(NA, nrow=p, ncol=p) bi = matrix(NA, nrow=p, ncol=p) for (j in 1:p) { wstar = rbind(matrix(wt), 1) logTstar = rbind(matrix(log(T)), 10^10) Xstar = rbind(X, -apply(X*matrix(rep(wt, times=p), ncol=p), 2, sum)+ 2*sqrt(n)*(Ln + E[, j])) bi[, j] = matrix(rq(logTstar ~ Xstar-1, tau=0.5, weights=wstar)$coefficients) if ( t(bi[, j])%*%Xstar[length(T)+1, ] >= 10^9 ) { E[, j] = -E[, j] wstar = rbind(matrix(wt), 1) logTstar = rbind(matrix(log(T)), 10^10) Xstar = rbind(X, -apply(X*matrix(rep(wt, times=p), ncol=p), 2, sum)+ 2*sqrt(n)*(Ln + E[, j])) bi[, j] = matrix(rq(logTstar ~ Xstar-1, tau=0.5, weights=wstar)$coefficients) if ( t(bi[, j])%*%Xstar[length(T)+1, ] >= 10^9 ) { D[, j] = 0 } else { D[, j] = bi[, j] - bb } } else { D[, j] = bi[, j] - bb } } if ( abs(det(D)) < 10^-10 ) { for (j in 1:p) { E[, j] = 2*E[, j] wstar = rbind(matrix(wt), 1) logTstar = rbind(matrix(log(T)), 10^10) Xstar = rbind(X, -apply(X*matrix(rep(wt, times=p), ncol=p), 2, sum)+ 2*sqrt(n)*(Ln + E[, j])) bi[, j] = matrix(rq(logTstar ~ Xstar-1, tau=0.5, weights=wstar)$coefficients) if ( t(bi[, j])%*%Xstar[length(T)+1, ] >= 10^9 ) { E[, j] = -E[, j] wstar = rbind(matrix(wt), 1) logTstar = rbind(matrix(log(T)), 10^10) Xstar = rbind(X, -apply(X*matrix(rep(wt, times=p), ncol=p), 2, sum)+ 2*sqrt(n)*(Ln + E[, j])) bi[, j] = matrix(rq(logTstar ~ Xstar-1, tau=0.5, weights=wstar)$coefficients) if ( t(bi[, j])%*%Xstar[length(T)+1, ] >= 10^9 ) { D[, j] = 0 } else { D[, j] = bi[, j] - bb } } else { D[, j] = bi[, j] - bb } if ( abs(det(D)) > 10^-10 ) { break } } } B = matrix(NA, p, p) J = matrix(NA, p, p) if ( abs(det(D)) >= 10^(-10)) { B = n^{-1/2} * E %*% solve(D) tildeE = matrix(NA, nrow=p, ncol=p) for (j in 1:p) { tildeE[, j] = apply(x*matrix(rep((L < exp(x%*%bi[, j]))&(exp(x%*%bi[, j]) <= R), times=p), ncol=p), 2, sum)/sqrt(n) - apply(x*matrix(rep((L < exp(x%*%bb))&(exp(x%*%bb) <= R), times=p), ncol=p), 2, sum)/sqrt(n) } J = n^{-1/2}*tildeE%*%solve(D) } return(list(B = B, J = J)) } ######################################################################################## ######################################################################################## #### The function VARf() to calculate variances of proposed estimators ########### ######################################################################################## VARf <- function(betaest, wt, wtVAR, X, T, grids, U, x, L, R, id, ID) { #### betaest: the estimator of beta #### wt: weight function, can be IPW, EEP, CC, and Full weights #### wtVAR: the AIPW weight function in the variance estimation #### N: a scalar denoted number of observations #### X: a N*p matrix denoted by X_i #### T: a N*1 vector denoted T_i^j #### grids: number of equally grids in (0, U] #### U: the upper bound of time interval (0, U] #### x: a n by p matrix denoted the covariate #### L: a n by 1 vector denoted the left censoring time #### R: a n by 1 vector denoted the right censoring time #### ID: a N*1 vector listed the ID of recurrent events #### id: a n*1 vector listed subject id VAR = matrix(NA, grids, p) Bmat = Jmat = array(NA, dim=c(grids, p, p)) for (j in 1:grids) { if (sum(is.na(betaest[j,]))==0) { Bmat[j,,] = BJf(wt, betaest[j,], X, T, x, L, R, ID, id)$B Jmat[j,,] = BJf(wt, betaest[j,], X, T, x, L, R, ID, id)$J } } for (i in 1:4) { if (is.na(Bmat[i, 1, 1]) == 1) { for (j in (i+1):5) { if (is.na(Jmat[j, 1, 1])==0) { Bmat[i,,] = Bmat[j,,] Jmat[i,,] = Jmat[j,,] break } } } } maxu = 0 for (i in 1:grids) { if (sum(is.na(Bmat[i,,]))==0) { maxu = i } else { break } } ### First calculate I if (maxu > 0) { I = array(NA, dim=c(grids+1, grids+1, p, p)) for (i in 1:maxu) { for (j in (i+1):(maxu+1)) { # if ( (j==(i+1)) & (1-is.na(det(Bmat[j-1,,]))) & (abs(det(Bmat[j-1,,])) > 10^-10 ) ) if (j==(i+1)) { I[i, j,,] = diag(p)+Jmat[j-1,,]%*%solve(Bmat[j-1,,])*U/grids ### suppose g(u)=1 } # if ( (j > (i+1)) & (1-is.na(det(Bmat[j-1,,]))) & (abs(det(Bmat[j-1,,])) > 10^-10 ) ) if ( j > (i+1) ) { I[i, j,,] = (diag(p)+Jmat[j-1,,]%*%solve(Bmat[j-1,,])*U/grids)%*%I[i, j-1,,] } } ## end of for j } ## end of for i ### Then calculate xi xi = array(NA, dim=c(n, grids, p)) for (i in 1:n) { for (j in 1:maxu) { temp1 = wtVAR*(T <= exp(X%*%betaest[j, ])) N1 = sum( temp1[which(ID==id[i])] ) if (j==1) { M1 = I(L[i]==0)*U/grids ### suppose g(u)=1 } if (j > 1) { M1 = M1 + (L[i] <= exp(x[i, ]%*%betaest[j-1, ]))*(exp(x[i, ]%*%betaest[j-1, ]) <= R[i])*U/grids } xi[i, j, ] = x[i, ]*(N1-M1) } ### end of loop j } ### end of loop i ### Calculate Phi Phi = array(0, dim=c(n, grids, p)) for (k in 1:n) { for (i in 1:maxu) { if (i==1) { Phi[k, 1, ] = I[1, 2, ,] %*% xi[k, 1, ] } if (i > 1) { Phi[k, i, ] = I[1, i+1, ,] %*% xi[k, 1, ] for (j in 2:i) { Phi[k, i, ] = Phi[k, i, ] + I[j, i+1, ,] %*% (xi[k, j, ]-xi[k, j-1, ]) } ### end for j } ### end if i>1 } ### end for i } ### end for k ### Calculate eta eta = array(0, dim=c(n, grids, p)) for (i in 1:n) { for (j in 1:maxu) { if (sum(is.na(Bmat[j,,])) == 0) { eta[i, j, ] = solve( Bmat[j,,] ) %*% Phi[i,j,] } } } ### Calculate the variance for (j in 1:maxu) { VAR[j, ] = diag( t(eta[, j, ])%*%eta[, j, ])/n } } ### end of if (maxu>0) return( sqrt(VAR/n) ) } ######################################################################################## ######################################################################################## ######### Second Part: The data Setup in Simulation Studies ##################### ######################################################################################## set.seed(100) n = 200 ##### sample size n NS = 500 ##### number of simulation NS U = 3 ##### the upper bound of time interval (0, U] grids = 150 ##### number of equally grids in (0, U] p = 3 ##### dimension of the covariate X link = "1" ##### the scale function g(s) = 1 rho01 = rho11 = rho21 = 1.5 ##### parameters to control the number of type 1 events rho02 = rho12 = rho22 = 2 ##### parameters to control the number of type 2 events alpha = c(1, 0.15) ##### the parameter to control the missing rate h0 = 0.8 ##### the bandwidth h in the kernel function #### the proposed IPW and EEP estimators for type 1 events estbeta1IPW = estbeta1EEP = array(NA, dim=c(NS, grids, p)) #### the proposed IPW and EEP estimators for type 2 events estbeta2IPW = estbeta2EEP = array(NA, dim=c(NS, grids, p)) #### the CC estimator and Full estimator for type 1 events estbeta1CC = estbeta1FF = array(NA, dim=c(NS, grids, p)) #### the CC estimator and Full estimator for type 2 events estbeta2CC = estbeta2FF = array(NA, dim=c(NS, grids, p)) #### the perturbation-based standard error estimators for type 1 events ESEmat1IPW = ESEmat1EEP = ESEmat1CC = ESEmat1FF = array(NA, dim=c(NS, grids, p)) #### the perturbation-based standard error estimators for type 2 events ESEmat2IPW = ESEmat2EEP = ESEmat2CC = ESEmat2FF = array(NA, dim=c(NS, grids, p)) #### the sample-based standard error estimators for type 1 events ESEsam1IPW = ESEsam1EEP = ESEsam1CC = ESEsam1FF = array(NA, dim=c(NS, grids, p)) #### the sample-based standard error estimators for type 2 events ESEsam2IPW = ESEsam2EEP = ESEsam2CC = ESEsam2FF = array(NA, dim=c(NS, grids, p)) for (k in 1:NS) { Z1 = rbinom(n, 1, 0.5) Z2 = runif(n, -0.5, 0.5) x = cbind(1, Z1, Z2) #### the n*p covariates L = rbinom(n, 1, 0.8)*runif(n, 0, 1) #### the n*1 left censoring time L R = rep(0, times=n) #### the n*1 right censoring time R gam1 = gam2 = rep(1, times=n) #### Case 1: gamma = 1 # gam1 = rgamma(n, shape=2, scale=0.5) #### Case 2: Gamma(2, 1/2), mean 1, variance 0.5 # gam2 = rgamma(n, shape=2, scale=0.5) #### Gamma(2, 1/2), mean 1, variance 0.5 data = NULL for (i in 1:n) { R[i] = runif(1, L[i], 12) Tstar1 = cumsum( rexp(100) ) Tstar2 = cumsum( rexp(100) ) T1 = rho01*Tstar1/gam1[i]*exp( pmin(1, rho11*Tstar1/1.5/gam1[i])*Z1[i] + rho21*Z2[i] ) T2 = rho02*Tstar2/gam2[i]*exp( pmin(1, rho12*Tstar2/1.5/gam2[i])*Z1[i] + rho22*Z2[i] ) T1 = T1[(T1 >= L[i])&(T1 <= R[i])] T2 = T2[(T2 >= L[i])&(T2 <= R[i])] dat = rbind(cbind(T1, rep(1, times=length(T1)), rep(0, times=length(T1))), cbind(T2, rep(0, times=length(T2)), rep(1, times=length(T2)))) sdat = matrix(dat[order(dat[, 1]), ], ncol=3) mi = length(sdat[, 1]) dai = cbind(matrix(rep(cbind(i, L[i], R[i], 1, Z1[i], Z2[i]), times=mi), ncol=6, byrow=TRUE), sdat) covaA = cbind(dai[, 5], dai[, 7]) Ai = rbinom(mi, 1, prob=1/( 1 + exp(covaA %*% alpha) ) ) data = rbind(data, cbind(dai, Ai)) } colnames(data) = c("id", "L", "R", "1", "Z1", "Z2", "T", "delta1", "delta2", "Amis") X = data[, 4:6] #### the N*p covariates T = data[, 7] #### the N*1 recurrent event time T_i^j delta1 = data[, 8] #### the N*1 vector delta_i1^j delta2 = data[, 9] #### the N*1 vector delta_i2^j A = 1 - data[, 10] #### the N*1 vector A_i^j N = length(A) #### the scalar denotes number of observations p = dim(X)[2] #### the dimension of the covariate X id = 1:n #### the n*1 vector listed subject id ID = data[, 1] #### the N*1 vector listed the ID of each observation zeta0 = rep(1, times=n) #### no perturbation in calculating the estimators ZETA0 = rep(1, times=N) #### no perturbation in calculating the estimators ppi = calpi(A, X, T, h=h0, ZETA0) #### the N*1 vector pi_i^j ppi[is.na(ppi)] = 1e-10 ppi[ppi==0] = 1e-10 pp = calpp(A, X, T, delta1, delta2, h=h0, ZETA0) #### the N*2 vector cbind(p_i1^j, p_i2^j) pp[is.na(pp)] = 0 wt1IPW = A/ppi*delta1 #### IPW weight for type 1 events wt1EEP = A*delta1+(1-A)*pp[, 1] #### EEP weight for type 1 events wtVAR1IPW = wtVAR1EEP = A/ppi*delta1+(1-A/ppi)*pp[, 1] #### AIPW weight for type 1 events wt1CC = wtVAR1CC = A*delta1 wt1FF = wtVAR1FF = delta1 wt2IPW = A/ppi*delta2 #### IPW weight for type 2 events wt2EEP = A*delta2+(1-A)*pp[, 2] #### EEP weight for type 2 events wtVAR2IPW = wtVAR2EEP = A/ppi*delta2+(1-A/ppi)*pp[, 2] #### AIPW weight for type 2 events wt2CC = wtVAR2CC = A*delta2 wt2FF = wtVAR2FF = delta2 estbeta1IPW[k,,] = estf(wt=wt1IPW, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta1EEP[k,,] = estf(wt=wt1EEP, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta2IPW[k,,] = estf(wt=wt2IPW, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta2EEP[k,,] = estf(wt=wt2EEP, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta1CC[k,,] = estf(wt=wt1CC, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta1FF[k,,] = estf(wt=wt1FF, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta2CC[k,,] = estf(wt=wt2CC, X, T, grids, link, U, x, L, R, zeta0, ZETA0) estbeta2FF[k,,] = estf(wt=wt2FF, X, T, grids, link, U, x, L, R, zeta0, ZETA0) NSBW = 100 boot1IPW = boot1EEP = boot1CC = boot1FF = array(NA, dim=c(NSBW, grids, p)) boot2IPW = boot2EEP = boot2CC = boot2FF = array(NA, dim=c(NSBW, grids, p)) for (l in 1:NSBW) { zeta = rexp(n, rate=1) ## exponetional variable to the perturbed variance est ZETA = rep(0, times=N) for (i in 1:n) { ZETA[which(data[, 1]==i)]=zeta[i] } ppizeta = calpi(A, X, T, h=h0, ZETA) ppizeta[is.na(ppizeta)] = 1e-10 ppizeta[ppizeta==0] = 1e-10 ppzeta = calpp(A, X, T, delta1, delta2, h=h0, ZETA) ppzeta[is.na(ppzeta)] = 0 wt1IPWzeta = A/ppizeta*delta1 wt1EEPzeta = A*delta1+(1-A)*ppzeta[, 1] wt2IPWzeta = A/ppizeta*delta2 wt2EEPzeta = A*delta2+(1-A)*ppzeta[, 2] boot1IPW[l,,] = (estf(wt=wt1IPWzeta, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta1IPW[k,,])^2 boot1EEP[l,,] = (estf(wt=wt1EEPzeta, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta1EEP[k,,])^2 boot2IPW[l,,] = (estf(wt=wt2IPWzeta, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta2IPW[k,,])^2 boot2EEP[l,,] = (estf(wt=wt2EEPzeta, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta2EEP[k,,])^2 boot1CC[l,,] = (estf(wt=wt1CC, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta1CC[k,,])^2 boot1FF[l,,] = (estf(wt=wt1FF, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta1FF[k,,])^2 boot2CC[l,,] = (estf(wt=wt2CC, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta2CC[k,,])^2 boot2FF[l,,] = (estf(wt=wt2FF, X, T, grids, link, U, x, L, R, zeta, ZETA)-estbeta2FF[k,,])^2 } ESEmat1IPW[k,,] = sqrt(apply(boot1IPW, 2:3, mean, na.rm=TRUE)) ESEmat1EEP[k,,] = sqrt(apply(boot1EEP, 2:3, mean, na.rm=TRUE)) ESEmat2IPW[k,,] = sqrt(apply(boot2IPW, 2:3, mean, na.rm=TRUE)) ESEmat2EEP[k,,] = sqrt(apply(boot2EEP, 2:3, mean, na.rm=TRUE)) ESEmat1CC[k,,] = sqrt(apply(boot1CC, 2:3, mean, na.rm=TRUE)) ESEmat1FF[k,,] = sqrt(apply(boot1FF, 2:3, mean, na.rm=TRUE)) ESEmat2CC[k,,] = sqrt(apply(boot2CC, 2:3, mean, na.rm=TRUE)) ESEmat2FF[k,,] = sqrt(apply(boot2FF, 2:3, mean, na.rm=TRUE)) ESEsam1IPW[k,,] = VARf(estbeta1IPW[k,,], wt1IPW, wtVAR1IPW, X, T, grids, U, x, L, R, id, ID) ESEsam1EEP[k,,] = VARf(estbeta1EEP[k,,], wt1EEP, wtVAR1EEP, X, T, grids, U, x, L, R, id, ID) ESEsam1CC[k,,] = VARf(estbeta1CC[k,,], wt1CC, wtVAR1CC, X, T, grids, U, x, L, R, id, ID) ESEsam1FF[k,,] = VARf(estbeta1FF[k,,], wt1FF, wtVAR1FF, X, T, grids, U, x, L, R, id, ID) ESEsam2IPW[k,,] = VARf(estbeta2IPW[k,,], wt2IPW, wtVAR2IPW, X, T, grids, U, x, L, R, id, ID) ESEsam2EEP[k,,] = VARf(estbeta2EEP[k,,], wt2EEP, wtVAR2EEP, X, T, grids, U, x, L, R, id, ID) ESEsam2CC[k,,] = VARf(estbeta2CC[k,,], wt2CC, wtVAR2CC, X, T, grids, U, x, L, R, id, ID) ESEsam2FF[k,,] = VARf(estbeta2FF[k,,], wt2FF, wtVAR2FF, X, T, grids, U, x, L, R, id, ID) print(k) write.csv( cbind(estbeta1IPW[k,,], ESEmat1IPW[k,,], ESEsam1IPW[k,,]), file = paste("est1IPWc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta1EEP[k,,],ESEmat1EEP[k,,], ESEsam1EEP[k,,]), file = paste("est1EEPc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta1CC[k,,],ESEmat1CC[k,,], ESEsam1CC[k,,]), file = paste("est1CCc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta1FF[k,,],ESEmat1FF[k,,], ESEsam1FF[k,,]), file = paste("est1FFc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta2IPW[k,,],ESEmat2IPW[k,,], ESEsam2IPW[k,,]), file = paste("est2IPWc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta2EEP[k,,],ESEmat2EEP[k,,], ESEsam2EEP[k,,]), file = paste("est2EEPc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta2CC[k,,],ESEmat2CC[k,,], ESEsam2CC[k,,]), file = paste("est2CCc_1", k, ".csv", sep="" )) write.csv( cbind(estbeta2FF[k,,],ESEmat2FF[k,,], ESEsam2FF[k,,]), file = paste("est2FFc_1", k, ".csv", sep="" )) } ########################################################################################