######################################################################################## ######## R code for "Quantile Regression Modeling of Latent Trajectory Features ###### ######## with Longitudinal Data" ###### ######## For academic research only ###### ######################################################################################## ## Data Notation ## N: sample size ## mm: the number of longitudinal observations. Here we assume mm is bounded above by 10. ## XX: baseline covariates, include 1, X1, X2 ## YY: the longitudinal response observations Yi1,..., Yi10. ## TT: the longitudinal time points ti1,..., ti10. ## Note: This program is for fixed bandwidth h=0.8. The longitudinal data is generated ## from a linear trajectory model, Yij = Ai + Bi tij + epsilon_ij. ## One reference is Wang et al. (2012, Biometrika)'s code from Huixia Wang's website rm(list=ls(all=TRUE)) ####################################################################################### ################# the package used in this paper ################################## ####################################################################################### library(quantreg) ####################################################################################### ######################################################################################## ############## First Part: R functions ################################# ######################################################################################## ######################################################################################## ###### rlaplace(): generate random variable from Laplace distribution ############## ######################################################################################## rlaplace <- function(N, mu, sigma) { ## generate N iid laplace random variables with mean mu, variance sigma^2 ## following the definition of laplace on wikipedia, laplace(mu, sigma)=laplace(mu, b) ## variance = 2*b^2 = sigma^2, so b = sigma/sqrt(2) b <- sigma/sqrt(2) u <- runif(N, 0, 1) x <- mu - b*sign(u-0.5)*log(1-2*abs(u-0.5)) return(x) } ######################################################################################## ######################################################################################## ############# qlaplace(): calculate quantile for Laplace distribution ############## ######################################################################################## qlaplace <- function(p, mu, sigma) { ## calculate p-th quantile for a laplace distribution with mean mu, variance sigma^2 ## sigma^2 = 2*b^2 is the variance b <- sigma/sqrt(2) mu - b*sign(p-0.5)*log(1-2*abs(p-0.5)) } ######################################################################################## ######################################################################################## ############# XgenerateB(): using quantile regression model to generate the ####### ############# random slope B, given the covariate X ####### ######################################################################################## XgenerateB <- function(c1, c2, c3, X1, X2) { ## Given X, generate random slope B eb <- rnorm(N, 0, 1) c1 + c2*X1 + c3*X2 + (0.1+X1+X2)*eb } ######################################################################################## ######################################################################################## ############# ABgenerateY(): generate YY, given A, B and TT ################# ######################################################################################## ABgenerateY <- function(A, B, TT, m, sigma) { ## Given random A and B, generate longitudinal data YY, this is a N*m matrix ## A: the random intercept ## B: the random slope ## TT: matrix induced by longitudinal time tij ## m is the upper bound for mi ## sigma^2: variance for error epsilon_ij AA <- matrix( rep(A, times=m), ncol=m ) BB <- matrix( rep(B, times=m), ncol=m ) #### epsilonij <- matrix( rnorm(N*m, 0, sigma), ncol=m) epsilonij <- matrix( rlaplace(N*m, 0, sigma), ncol=m ) #### epsilonij <- matrix( runif(N*m, -sqrt(3)/2, sqrt(3)/2), ncol=m ) AA + BB*TT + epsilonij } ######################################################################################## ######################################################################################## ############# estsigma(): the function to estimate sigma^2 ################# ######################################################################################## estsigma <- function(YY, TT, mm, ww) { ## estimate the variance sigma^2 of epsilon_ij ## YY: longitudinal response, N*m matrix ## TT: longitudinal time points, N*m matrix ## mm: number of longitudinal observations, N*1 vector ## ww: the perturbed variables, 1 or from Exp(1). RSS <- rep(0, times = N ) for (i in 1:N) { ni <- mm[i] Yi <- YY[i, 1:ni] Zi <- matrix( c( rep(1, times=ni), TT[i, 1:ni]), ncol=2 ) Pi <- Zi%*%solve( t(Zi)%*%Zi )%*%t(Zi) RSS[i] <- t(Yi)%*%( diag(ni) - Pi )%*%Yi } sqrt( 1/(sum(mm)-2*N)*sum(ww*RSS)/mean(ww) ) } ######################################################################################## ######################################################################################## ############# estsubjcoef(): the function to calculate D_i ############### ######################################################################################## estsubjcoef <- function(TT, mm) { # mm is a N*1 vector, denote the number in each subject subjcoef <- rep(0, times=N) for (i in 1:N) { ni <- mm[i] Zi <- matrix( c( rep(1, times=ni), TT[i, 1:ni]), ncol=2 ) SZZi <- solve( t(Zi)%*%Zi ) subjcoef[i] <- sqrt( SZZi[2, 2] ) } return( subjcoef ) } ######################################################################################## ############# estB(): least square naive estimator \hat{B}_i ############### ######################################################################################## estB <- function(YY, TT, mm) { # mm is a N*1 vector, denote the number in each subject hatB <- rep(0, times = N ) for (i in 1:N) { ni <- mm[i] Yi <- YY[i, 1:ni] Zi <- cbind(1, TT[i, 1:ni]) hatcoef <- solve( t(Zi)%*%Zi )%*%t(Zi)%*%Yi hatB[i] <- hatcoef[2] } return( hatB ) } ######################################################################################## ######################################################################################## ############# K(): the smooth function ############### ######################################################################################## K <- function(u, type) { if (type=="norm") pnorm(u) # normal CDF else if (type=="poly") 1*( abs(u)<=1 )*( 0.5+105/64*(u-5*u^3/3+7*u^5/5-3*u^7/7) ) + 1*(u > 1) else if (type=="exp") exp(u)/( 1+exp(u) ) # the kernel function used: page 1337 of Horowitz (1996) } ######################################################################################## ######################################################################################## ############# K1(): the first derivative function of K() ############### ######################################################################################## K1 <- function(u, type) { # first derivative of K if (type=="norm") dnorm(u) else if (type=="poly") 1*( abs(u) <= 1 )*105/64*(1-5*u^2+7*u^4-3*u^6) else if (type=="exp") exp(u)/( 1+exp(u) )^2 } ######################################################################################## ######################################################################################## ############# K2(): the second derivative function of K() ############### ######################################################################################## K2 <- function(u, type) { # second derivative of K if (type=="norm") -u*dnorm(u) else if (type=="poly") 1*( abs(u)<=1 )*105/64*(-10*u+28*u^3-18*u^5) else if (type=="exp") exp(u)*( 1-exp(2*u) )/( 1+exp(u) )^4 } ######################################################################################## ######################################################################################## ############# CSobj.CL(): the objective function ############### ######################################################################################## CSobj.CL <- function(beta, hatB, sigma, XX, tau, h, ww, subjnorm) { u <- (hatB - XX %*% beta)/subjnorm part1 <- u*( tau - 1 + K(u/h, type="norm") ) part2 <- sigma^2/2*( 2/h*K1(u/h, type="norm") + u/(h^2)*K2(u/h, type="norm") ) out <- sum( ww*(part1-part2) ) return( out ) } ######################################################################################## ######################################################################################## ############# estbeta.CN0(): minimizer of the objective function ############### ######################################################################################## estbeta.CN0 <- function(betaini, hatB, sigma, XX, tau, h, ww, subjnorm) { betahat <- optim(par=betaini, fn = CSobj.CL, hatB=hatB, sigma=sigma, XX=XX, tau=tau, h=h, ww=ww, subjnorm=subjnorm)$par return( betahat ) } ######################################################################################## ######################################################################################## ######### Second Part: The data Setup in Simulation Studies ################# ######### basic set-up when the sample size N=200 ################# ######################################################################################## set.seed(100) N <- 200 ###### sample size NS <- 1000 ###### number of simulations c1 <- 2 ###### parameter in quantile regression model c2 <- 1 ###### parameter in quantile regression model c3 <- 1 ###### parameter in quantile regression model esttau <- seq(0.1, 0.9, by=0.1) ###### estimated quantile levels tau ntau <- length( esttau ) ###### number of estimated quantile levels ######################################################################################## beta.true <- matrix(0, nrow=ntau, ncol=3) ### true quantile regression parameters naivematrix <- estmatrix0 <- array(0, dim=c(ntau, NS, 3) ) ESEmatrix <- ESEnaivematrix <- array(0, dim=c(ntau, NS, 3) ) CP.est <- se.est <- CP.naive <- se.naive <- matrix(0, nrow=ntau, ncol=3) for (j in 1:NS) { print(j) ## the j-th simulation X1 <- runif(N, 0, 0.5) ## baseline covariate X2 <- rbinom(N, 1, 0.5) ## baseline covariate XX <- cbind(1, X1, X2) A <- rexp(N, 0.8) ## random intercept B <- XgenerateB(c1=c1, c2=c2, c3=c3, X1=X1, X2=X2) ## quantile regression model generate slope m <- 10 ## the maximum longitudinal observations tij <- matrix( rexp(N*m, 0.8), ncol=m) TT <- t( apply(tij, 1, cumsum) ) ## longitudinal time points YY <- ABgenerateY(A=A, B=B, TT=TT, m=m, sigma=1) ## longitudinal response observations mm <- 4 + floor( runif(N, 0, m-4) ) ## number of longitudinal observations estsig <- estsigma(YY, TT, mm, ww=1) ## hat{sigma} subjnorm <- estsubjcoef(TT, mm) ## D_i in paper hatB <- estB(YY, TT, mm) ## hat{B} h <- 0.8 ## fixed bandwidth for (i in 1:ntau) { betaini <- rq( hatB~XX-1, tau=esttau[i] )$coeff ## naive estimator naivematrix[i, j, ] <- betaini beta.true[i, ] <- c( c1+0.1*qnorm(esttau[i]), c2+qnorm(esttau[i]), c3+qnorm(esttau[i]) ) estmatrix0[i, j, ] <- estbeta.CN0(betaini=betaini, hatB=hatB, sigma=estsig, XX=XX, tau=esttau[i], h=h, ww=1, subjnorm=subjnorm) ## proposed estimator NSBW <- 200 ### number of perturbation boot.est <- boot.naive <- matrix(0, NSBW, 3) for (k in 1:NSBW) ## This loop is for standard error estimation { ww <- rexp(N, 1) ## perturbed random variables from Exp(1) estsigk <- estsigma(YY, TT, mm, ww=ww) ## sigma* boot.est[k, ] <- estbeta.CN0(betaini=betaini, hatB=hatB, sigma=estsigk, XX=XX, tau=esttau[i], h=h, ww=ww, subjnorm=subjnorm) ## proposed beta* boot.naive[k, ] <- rq( hatB~XX-1, tau=esttau[i], weights=ww )$coeff ## naive beta* } ESEmatrix[i,j, ] <- apply(boot.est, 2, sd) ## standard error for proposed estimator ESEnaivematrix[i,j, ] <- apply(boot.naive, 2, sd) ## standard error for naive estimator } } beta.naive <- cbind( apply(naivematrix[, , 1], 1, mean), apply(naivematrix[, , 2], 1, mean), apply(naivematrix[, , 3], 1, mean) ) sd.naive <- cbind( apply(naivematrix[, , 1], 1, sd), apply(naivematrix[, , 2], 1, sd), apply(naivematrix[, , 3], 1, sd) ) beta.est0 <- cbind( apply(estmatrix0[, , 1], 1, mean), apply(estmatrix0[, , 2], 1, mean), apply(estmatrix0[, , 3], 1, mean) ) sd.est0 <- cbind( apply(estmatrix0[, , 1], 1, sd), apply(estmatrix0[, , 2], 1, sd), apply(estmatrix0[, , 3], 1, sd) ) for (i in 1:ntau) { truematrix <- matrix( rep(beta.true[i, ], each=NS), nrow=NS ) CPmatrix <- ( truematrix - 1.96*ESEmatrix[i, ,] <= estmatrix0[i, ,] )*( estmatrix0[i, ,] <= truematrix + 1.96*ESEmatrix[i, ,] ) CP.est[i, ] <- apply(CPmatrix, 2, mean)*100 se.est[i, ] <- apply(ESEmatrix[i, ,], 2, mean) CPnaivematrix <- ( truematrix - 1.96*ESEnaivematrix[i, ,] <= naivematrix[i, ,] )*( naivematrix[i, ,] <= truematrix + 1.96*ESEnaivematrix[i, ,] ) CP.naive[i, ] <- apply(CPnaivematrix, 2, mean)*100 se.naive[i, ] <- apply(ESEnaivematrix[i, ,], 2, mean) } beta.est0-beta.true #### bias of the proposed estimator beta.naive-beta.true #### bias of the naive estimator se.est #### estimated empirical error of the proposed estimator sd.est0 #### empirical derivation of the proposed estimator se.naive #### estimated empirical error of the naive estimator sd.naive #### empirical derivation of the naive estimator CP.est #### empirical coverage probability of the proposed estimator CP.naive #### empirical coverage probability of the naive estimator par( mfrow=c(3, 3) ) plot( esttau, beta.est0[,1]-beta.true[,1], xlim=c(0,1), ylim=c(-0.2, 0.2), type="l", lty=1, col="red", xlab="tau", ylab="Bias" ) lines(esttau, beta.naive[,1]-beta.true[,1] , type="l", lty=2, col="black" ) legend("topleft", expression("Proposed", "Naive"), lty=c(1, 2), col=c("red", "black") ) plot( esttau, beta.est0[,2]-beta.true[,2], xlim=c(0,1), ylim=c(-0.2, 0.2), type="l", lty=1, col="red", xlab="tau", ylab="Bias" ) lines( esttau, beta.naive[,2]-beta.true[,2] , type="l", lty=2, col="black" ) legend("topright", expression("Proposed", "Naive"), lty=c(1, 2), col=c("red", "black") ) plot( esttau, beta.est0[,3]-beta.true[,3], xlim=c(0,1), ylim=c(-0.1, 0.1), type="l", lty=1, col="red", xlab="tau", ylab="Bias" ) lines( esttau, beta.naive[,3]-beta.true[,3] , type="l", lty=2, col="black" ) legend("topright", expression("Proposed", "Naive"), lty=c(1, 2), col=c("red", "black") ) plot( esttau, se.est[,1], xlim=c(0,1), ylim=c(0, 0.2), type="l", lty=1, col="red", xlab="tau", ylab="Empirical Standard Error" ) plot( esttau, se.est[,2], xlim=c(0,1), ylim=c(0, 1), type="l", lty=1, col="red", xlab="tau", ylab="Empirical Standard Error" ) plot( esttau, se.est[,3], xlim=c(0,1), ylim=c(0, 0.4), type="l", lty=1, col="red", xlab="tau", ylab="Empirical Standard Error" ) plot( esttau, CP.est[, 1], xlim=c(0,1), ylim=c(0, 100), type="l", lty=1, col="red", xlab="tau", ylab="Coverage Probability" ) plot( esttau, CP.est[, 2], xlim=c(0,1), ylim=c(0, 100), type="l", lty=1, col="red", xlab="tau", ylab="Coverage Probability" ) plot( esttau, CP.est[, 3], xlim=c(0,1), ylim=c(0, 100), type="l", lty=1, col="red", xlab="tau", ylab="Coverage Probability" )