QRegVCM (version 1.2)

AHeVT: AHe V(t)-approach

Description

The adapted He (1997) approach considering a simple hetersocedastic varying-coefficient model, V(t). $$Y(t)=\sum_{k=0}^{p}\beta_{k}(t)X^{(k)}(t)+V(t)\varepsilon(t)$$.

Usage

AHeVT(VecX, times, subj, X, y, d, tau, kn, degree, lambda, gam)

Arguments

VecX

The representative values for each covariate used to estimate the desired conditional quantile curves.

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.

tau

The quantiles of interest.

kn

The number of knots for each covariate.

degree

The degree of B-spline basis for each covariate.

lambda

The grid of smoothing parameter to control the trade of 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).

Value

hat_bt50

The median coefficients estimators.

hat_VT

The variability estimator.

C

The estimators of the tau-th quantile of the estimated residuals.

qhat

The conditional quantile curves estimator.

References

Andriyana, Y., Gijbels, I., and Verhasselt, A. P-splines quantile regression estimation in varying coefficient models. Test 23, 1 (2014a),153--194.

Andriyana, Y., Gijbels, I. and Verhasselt, A. (2014b). Quantile regression in varying coefficient models: non-crossingness and heteroscedasticity. Manuscript.

He, X. (1997). Quantile curves without crossing. The American Statistician, 51, 186--192.

See Also

rq.fit.sfn as.matrix.csr truncSP

Examples

# NOT RUN {
<!-- %library(truncSP) -->
# }
# NOT RUN {
data(PM10)

PM10 = PM10[order(PM10$day,PM10$hour,decreasing=FALSE),]

y = PM10$PM10[1:200]
times = PM10$hour[1:200]
subj = PM10$day[1:200]
dim = length(y)
x0 = rep(1,200)
x1 = PM10$cars[1:200]
x2 = PM10$wind.speed[1:200]

X = cbind(x0, x1, x2)

VecX = c(1, max(x1), max(x2))

##########################
#### Input parameters ####
##########################
kn = c(10, 10, 10)
degree = c(3, 3, 3)
taus = seq(0.1,0.9,0.1)
lambdas = c(1,1.5,2)
d = c(1, 1, 1)
gam = 1/2
##########################


AHe = AHeVT(VecX=VecX, times=times, subj=subj, X=X, y=y, d=d, tau=taus,
      kn=kn, degree=degree, lambda=lambdas, gam=gam)
hat_bt50 = AHe$hat_bt50
hat_VT = AHe$hat_Vt
C = AHe$C
qhat = AHe$qhat

qhat1 = qhat[,1]
qhat2 = qhat[,2]
qhat3 = qhat[,3]
qhat4 = qhat[,4]
qhat5 = qhat[,5]
qhat6 = qhat[,6]
qhat7 = qhat[,7]
qhat8 = qhat[,8]
qhat9 = qhat[,9]

hat_bt0 = hat_bt50[seq(1,dim)]
hat_bt1 = hat_bt50[seq((dim+1),(2*dim))]
hat_bt2 = hat_bt50[seq((2*dim+1),(3*dim))]

i = order(times, hat_VT, qhat1, qhat2, qhat3, qhat4, qhat5, qhat6, qhat7,
    qhat8, qhat9, hat_bt0, hat_bt1, hat_bt2);
times = times[i]; hat_VT=hat_VT[i]; qhat1 = qhat1[i]; qhat2=qhat2[i];
qhat3=qhat3[i]; qhat4=qhat4[i]; qhat5=qhat5[i]; qhat6=qhat6[i];
qhat7=qhat7[i]; qhat8=qhat8[i]; qhat9=qhat9[i];
hat_bt0=hat_bt0[i];  hat_bt1=hat_bt1[i];  hat_bt2=hat_bt2[i];


### Plot coefficients

plot(hat_bt0~times, lwd=2, type="l", xlab="hour", ylab="baseline PM10");
plot(hat_bt1~times, lwd=2, type="l", xlab="hour",
    ylab="coefficient of cars");
plot(hat_bt2~times, lwd=2, type="l", xlab="hour",
    ylab="coefficient of wind");


### Plot variability V(t)

plot(hat_VT~times, ylim=c(min(hat_VT), max(hat_VT)), xlab="hour",
    ylab="", type="l", lwd=2);
mtext(expression(hat(V)(t)), side=2, cex=1, line=3)

### Plot conditional quantiles estimators

ylim = c(-4, 6)
plot(qhat1~times, col="magenta", cex=0.2, lty=5, lwd=2, type="l",
    ylim=ylim, xlab="hour", ylab="PM10");
lines(qhat2~times, col="aquamarine4", cex=0.2, lty=4, lwd=2);
lines(qhat3~times, col="blue", cex=0.2, lty=3, lwd=3);
lines(qhat4~times, col="brown", cex=0.2, lty=2, lwd=2);
lines(qhat5~times, col="black", cex=0.2, lty=1, lwd=2);
lines(qhat6~times, col="orange", cex=0.2, lty=2, lwd=2)
lines(qhat7~times, col="darkcyan", cex=0.2, lty=3, lwd=3);
lines(qhat8~times, col="green", cex=0.2, lty=4, lwd=2);
lines(qhat9~times, col="red", cex=0.2, lty=5, lwd=3)

legend("bottom", c(expression(tau==0.9), expression(tau==0.8),
    expression(tau==0.7), expression(tau==0.6), expression(tau==0.5),
    expression(tau==0.4), expression(tau==0.3), expression(tau==0.2),
    expression(tau==0.1)), ncol=1, col=c("red","green","darkcyan",
    "orange","black","brown","blue","aquamarine4","magenta"),
    lwd=c(2,2,3,2,2,2,3,2,2), lty=c(5,4,3,2,1,2,3,4,5))

# }