Variability Estimation and Testing
Estimating and Testing the variability function using the following hetersocedastic varying-coefficient model. [REMOVE_ME]
where V(X(t),t) one of the six variability function in Gijbels etal (2017).
Estimating and Testing the variability function using the following hetersocedastic varying-coefficient model.
where V(X(t),t) one of the six variability function in Gijbels etal (2017).
test_variability(times, subj, X, y, d, kn, degree, lambda, gam,tau, nr.bootstrap.samples, seed, test, model)
times: The vector of time variable.subj: The vector of subjects/individuals.X: The covariate containing 1 as its first component (including intercept in the model).y: The response vector.d: The order of differencing operator for each covariate.kn: The number of knot intervals for each covariate.degree: The degree of B-spline basis for each covariate.lambda: The grid of smoothing parameter to control the trade between fidelity and penalty term (use a fine grid of lambda).gam: The power used in estimating the smooting parameter for each covariate (e.g. gam=1 or gam=0.5).tau: The quantiles of interest.nr.bootstrap.samples: The number of bootstrap samples used.seed: The seed for the random generator in the bootstrap resampling.test: To request for testing the specific shape of the variability function ("Y" for test and "N" for only estimation of the parameters, the default is "Y").model: The variability model used to estimate the quantile of errors (the default is 4, model 4).est_median: the median estimator.
hat_bt50: The median coefficients estimators.
qhat5_s2_m0: The variability (model 0) estimator.
qhat5_s2_m1: The variability (model 1) estimator.
qhat5_s2_m2: The variability (model 2) estimator.
qhat5_s2_m3: The variability (model 3) estimator.
qhat5_s2_m4: The variability (model 4) estimator.
qhat5_s2_m5: The variability (model 5) estimator.
hat_btV_0: The variability coefficients (model 0) estimators.
hat_btV_1: The variability coefficients (model 1) estimators.
hat_btV_2: The variability coefficients (model 2) estimators.
hat_btV_3: The variability coefficients (model 3) estimators.
hat_btV_4: The variability coefficients (model 4) estimators.
hat_btV_5: The variability coefficients (model 5) estimators.
C: The estimators of the tau-th quantile of the estimated residuals.
comp: The pairwise comparisons for testing the variabilty function.
P: The p-values.
GR: The test statistics for the given data.
Gb: The bootstrap test statistics.
Andriyana, Y. and Gijbels, I. & Verhasselt, A. (2014). P-splines quantile regression estimation in varying coefficient models. Test, 23, 153-194.
Andriyana, Y., Gijbels, I. and Verhasselt, A. (2017). Quantile regression in varying-coefficient models: non-crossing quantile curves and heteroscedasticity. Statistical Papers, DOI:10.1007/s00362-016-0847-7.
Gijbels, I., Ibrahim, M. A., and Verhasselt, A. (2017). Testing the heteroscedastic error structure in quantile varying coefficient models. The Canadian Journal of Statistics, DOI:10.1002/cjs.11346.
He, X. (1997). Quantile curves without crossing. The American Statistician, 51, 186-192.
Mohammed Abdulkerim Ibrahim
Some warning messages are related to the function rq.fit.sfn.
rq.fit.sfn
as.matrix.csr
############################################################################ ##### The real data Example in Section S3 of the supplementary material ############################################################################ data(PM10) PM10 = PM10[order(PM10$day,PM10$hour,decreasing=FALSE),] y = PM10$PM10 ## the logarithm of the concentration of PM10 times = PM10$hour ## the time in hours subj = PM10$day ## subject indicator (day) dim=length(y) ## number of rows in the data = 500 x0 = rep(1,dim) ## for intercept # the covariates x1 = PM10$cars ## logarithm of number of cars per hour x2 = PM10$wind.speed ## the wind speed (in meters/second) x3 = PM10$temp ## the temperature (in degree Celsius) X = cbind(x0, x1, x2, x3) ## the covariate matrix px=ncol(X) ########################## ### Input parameters #### ######################### lambda = 1 # used 10^seq(-2, 1, 0.1) in Gijbels etal (2017) kn = rep(3,px) # used rep(10,px) in Gijbels etal (2017) degree = rep(3,px) d = rep(1,px) gam=0.25 nr.bootstrap.samples= 4 # used 200 in Gijbels etal (2017) seed=110 taus = 0.1 ######################### test1=test_variability(times=times, subj=subj, X=X, y=y, d=d, kn=kn, degree=degree, lambda=lambda, gam=gam, tau=taus, nr.bootstrap.samples=nr.bootstrap.samples,seed=seed, test="Y",model=4) #### Testing results test1$comp #the comparisons test1$P ## p-values test1$GR ## test statistics ### estimation results qhat5_s2_m4=test1$qhat5_s2_m4 qhat5_s2_m5=test1$qhat5_s2_m5 qhat5_s2_m0=test1$qhat5_s2_m0*rep(1,dim) gamma0=test1$hat_btV_4[1:dim] gamma1=test1$hat_btV_4[(dim+1):(dim*2)] gamma2=test1$hat_btV_4[(dim*2+1):(dim*3)] gamma3=test1$hat_btV_4[(dim*3+1):(dim*4)] i = order(times, qhat5_s2_m4, qhat5_s2_m5, qhat5_s2_m0,gamma0,gamma1, gamma2,gamma3); times_o = times[i]; qhat5_s2_m4_o=qhat5_s2_m4[i]; qhat5_s2_m5_o=qhat5_s2_m5[i]; qhat5_s2_m0_o=qhat5_s2_m0[i]; gamma0_o=gamma0[i]; gamma1_o=gamma1[i]; gamma2_o=gamma2[i];gamma3_o=gamma3[i] ##### variability functions plots plot(qhat5_s2_m4_o~times_o, col="magenta", cex=0.2, xlab="hour", ylab="estimated variability function") lines(qhat5_s2_m5_o~times_o, col="red", cex=0.2, lty=1, lwd=2); lines(qhat5_s2_m0_o~times_o, col="black", cex=0.2, lty=5, lwd=2); legend("topleft", c("Model 4", "Model 5", "Model 0"), ncol=1, col=c("magenta","red","black"), lwd=c(1,2,2), lty=c(3,1,5)) ### Plot of coefficients for variability function plot(gamma0_o~times_o, lwd=2, type="l", xlab="hour", ylab=expression(hat(gamma)(T))); plot(gamma1_o~times_o, lwd=2, type="l", xlab="hour", ylab="coefficient of logarithm of number of cars per hour"); plot(gamma2_o~times_o, lwd=2, type="l", xlab="hour", ylab="coefficient of wind speed"); plot(gamma3_o~times_o, lwd=2, type="l", xlab="hour", ylab="coefficient of temperature") ## Not run: ############################################################################### ############### The real data Example in Section 6 of Gijbels etal (2017) ############################################################################### data(CD4) subj = CD4$subj ## subject indicator (a man) dim = length(subj) ## number of rows in the data = 1817 y = CD4$CD4 ## the CD4 percentage X0 = rep(1,dim) ## the intercept X1 = CD4$Smooking ## the smoking status X2 = CD4$Age ## age at HIV infection X3 = CD4$PreCD4 ## the pre-infection CD4 percentage times = CD4$Time ## the time in years X = cbind(X0, X1, X2, X3) ## the covariate matrix px=ncol(X) lambdas = c(0.01,1,10) # used 10^seq(-2, 1, 0.1) in Gijbels etal (2017) kn = rep(10,px) # the number of internal knots for each covariate degree = rep(3,px) # the degree of splines d = rep(1,px) ## The differencing order in the penalty term for each covariate gam=0.25 ## the smooting parameter for each covariate nr.bootstrap.samples=100 ## used 200 in Gijbels etal (2017) seed=110 ## the seed for the random generator in the bootstrap resampling taus = seq(0.1,0.9,0.2) test2=test_variability(times=times, subj=subj, X=X, y=y, d=d, kn=kn, degree=degree, lambda=lambdas, gam=gam,tau=taus, nr.bootstrap.samples=nr.bootstrap.samples,seed=seed, test="Y",model=4) test2$comp test2$P ## p-values test2$GR ## test statistics ### estimation results qhat5_s2_m4=test2$qhat5_s2_m4 qhat5_s2_m5=test2$qhat5_s2_m5 qhat5_s2_m0=test2$qhat5_s2_m0*rep(1,dim) gamma0=test2$hat_btV_4[1:dim] gamma1=test2$hat_btV_4[(dim+1):(dim*2)] gamma2=test2$hat_btV_4[(dim*2+1):(dim*3)] gamma3=test2$hat_btV_4[(dim*3+1):(dim*4)] i = order(times, qhat5_s2_m4, qhat5_s2_m5, qhat5_s2_m0,gamma0,gamma1, gamma2,gamma3); times_o = times[i]; qhat5_s2_m4_o=qhat5_s2_m4[i]; qhat5_s2_m5_o=qhat5_s2_m5[i] qhat5_s2_m0_o=qhat5_s2_m0[i]; gamma0_o=gamma0[i]; gamma1_o=gamma1[i]; gamma2_o=gamma2[i];gamma3_o=gamma3[i] ##### variability functions plots plot(qhat5_s2_m4_o~times_o, col="black", cex=0.2, xlab="time since infection", ylab="estimated variability function") lines(qhat5_s2_m5_o~times_o, col="red", cex=0.2, lty=5, lwd=2); lines(qhat5_s2_m0_o~times_o, col="magenta", cex=0.2, lty=1, lwd=2); legend("topleft", c("Model 4", "Model 5", "Model 0"), ncol=1, col=c("black","red","magenta"), lwd=c(1,2,2), lty=c(3,5,1)) ### Plot of coefficients for variability function plot(gamma0_o~times_o, lwd=2, type="l", xlab="time since infection", ylab=expression(hat(gamma)(T))); plot(gamma1_o~times_o, lwd=2, type="l", xlab="time since infection", ylab="coefficient of smoking status"); plot(gamma2_o~times_o, lwd=2, type="l", xlab="time since infection", ylab="coefficient of age"); plot(gamma3_o~times_o, lwd=2, type="l", xlab="time since infection", ylab="coefficient of pre-infection CD4") ## End(Not run)
Useful links