Title: | Quantile Regression Using Asymmetric Laplace Distribution |
---|---|
Description: | EM algorithm for estimation of parameters and other methods in a quantile regression. |
Authors: | Luis Benites Sanchez, Christian E. Galarza, Victor H. Lachos |
Maintainer: | Luis Benites Sanchez <[email protected]> |
License: | GPL (>= 3.0) |
Version: | 1.0 |
Built: | 2024-11-16 03:03:10 UTC |
Source: | https://github.com/lbenitesanchez/aldqr |
Data on 102 male and 100 female athletes collected at the Australian Institute of Sport.
This data frame contains the following columns:
(0 = male or 1 = female)
height (cm)
weight (kg)
lean body mass
red cell count
white cell count
Hematocrit
Hemoglobin
plasma ferritin concentration
body mass index, weight/(height)**2
sum of skin folds
Percent body fat
Case Labels
Sport
S. Weisberg (2005). Applied Linear Regression, 3rd edition. New York: Wiley, Section 6.4
Return case-deletion estimating the parameters in a quantile regression
diag.qr(y,x,tau,theta)
diag.qr(y,x,tau,theta)
y |
vector of responses |
x |
the design matrix |
tau |
the quantile to be estimated, this is generally a number strictly between 0 and 1. |
theta |
parameter estimated |
Hessian and gradient matrix. Also the generalized cook distance (GDi), approximation of the likelihood distance (QDi)
Luis Benites Sanchez [email protected], Christian E. Galarza [email protected], Victor Hugo Lachos [email protected]
[1] Koenker, R. W. (2005). Quantile Regression, Cambridge U. Press. [2] Yu, K. & Moyeed, R. (2001). Bayesian quantile regression. Statistics & Probability Letters, 54 (4), 437 to 447. [3] Kotz, S., Kozubowski, T. & Podgorski, K. (2001). The laplace distribution and generalizations: A revisit with applications to communications, economics, engineering, and finance. Number 183. Birkhauser.
## Not run: ############################################################## ### Graphic of the generalized Cook distance for data(AIS) ### ############################################################## #Dados data(ais, package="sn") attach(ais) sexInd <- (sex=="female") + 0 x <- cbind(1,LBM,sexInd) y <- BMI #Percentile perc <- 0.5 res <- EM.qr(y,x,perc) diag <- diag.qr(y,x,perc,res$theta) HessianMatrix <- diag$MatrizQ Gradiente <- diag$mdelta GDI <- c() for (i in 1:202) { GDI[i] <- t(Gradiente[,i]) } #Seccion de los graficos par(mfrow = c(1,1)) plot(seq(1:202),GDI,xlab='Index',ylab=expression(paste(GD[i])),main='p=0.1') abline(h=2*(4+1)/202,lty=2) identify(GDI,n=1) plot(seq(1:202),GDI,xlab='Index',ylab=expression(paste(GD[i])),main='p=0.5') abline(h=2*(4+1)/202,lty=2) identify(GDI,n=1) plot(seq(1:202),GDI,xlab='Index',ylab=expression(paste(GD[i])),main='p=0.9') abline(h=2*(4+1)/202,lty=2) identify(GDI,n=4) ############################################################# ### Graphic of the likelihood displacemente for data(AIS) ### ############################################################# #Dados data(ais, package="sn"); attach(ais); sexInd<-(sex=="female")+0; x=cbind(1,LBM,sexInd); y=BMI #Percentile perc <- 0.9 n <- nrow(x) res <- EM.qr(y,x,perc) thetaest <- res$theta sigmaest <- thetaest[4] betaest <- matrix(thetaest[1:3],3,1) taup2 <- (2/(perc*(1-perc))) thep <- (1-2*peGraphic of the generalized Cook distance for data(AIS)rc)/(perc*(1-perc)) diag <- diag.qr(y,x,perc,thetaest) HessianMatrix <- diag$MatrizQ Gradiente <- diag$mdelta sigma <- sigmaest beta <- betaest muc <- (y-x delta2 <- (y-x gamma2 <- (2+thep^2/taup2)/sigma vchpN <- besselK(sqrt(delta2*gamma2), 0.5-1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2))^(-1) vchp1 <- besselK(sqrt(delta2*gamma2), 0.5+1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2)) Q <- -0.5*n*log(sigmaest)-0.5*(sigmaest*taup2)^{-1}* (sum(vchpN*muc^2 - 2*muc*thep + vchp1*(thep^2+2*taup2))) ######################################################## theta_i <- thetaest sigmaest <- theta_i[4,] betaest <- theta_i[1:3,] sigma <- sigmaest beta <- betaest muc <- (y-x delta2 <- (y-x gamma2 <- (2+thep^2/taup2)/sigma vchpN <- besselK(sqrt(delta2*gamma2), 0.5-1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2))^(-1) vchp1 <- besselK(sqrt(delta2*gamma2), 0.5+1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2)) Q1 <- c() for (i in 1:202) { Q1[i] <- -0.5*n*log(sigmaest[i])-sum(vchpN[,i]*muc[,i]^2 - 2*muc[,i]*thep + vchp1[,i]*(thep^2+2*taup2))/(2*(sigmaest[i]*taup2)) } ######################################################## QDi <- 2*(-Q+Q1) #Depois de escolger perc guardamos os valores de QDi QDi0.1 <- QDi QDi0.5 <- QDi QDi0.9 <- QDi #Seccion de los graficos par(mfrow = c(1,1)) plot(seq(1:202),QDi0.1,xlab='Index',ylab=expression(paste(QD[i])),main='p=0.1') abline(h=mean(QDi0.1)+3.5*sd(QDi0.1),lty=2) identify(QDi0.1,n=3) plot(seq(1:202),QDi0.5,xlab='Index',ylab=expression(paste(QD[i])),main='p=0.5') abline(h=mean(QDi0.5)+3.5*sd(QDi0.5),lty=2) identify(QDi0.5,n=3) plot(seq(1:202),QDi0.9,xlab='Index',ylab=expression(paste(QD[i])),main='p=0.9') abline(h=mean(QDi0.9)+3.5*sd(QDi0.9),lty=2) identify(QDi0.9,n=4) ## End(Not run)
## Not run: ############################################################## ### Graphic of the generalized Cook distance for data(AIS) ### ############################################################## #Dados data(ais, package="sn") attach(ais) sexInd <- (sex=="female") + 0 x <- cbind(1,LBM,sexInd) y <- BMI #Percentile perc <- 0.5 res <- EM.qr(y,x,perc) diag <- diag.qr(y,x,perc,res$theta) HessianMatrix <- diag$MatrizQ Gradiente <- diag$mdelta GDI <- c() for (i in 1:202) { GDI[i] <- t(Gradiente[,i]) } #Seccion de los graficos par(mfrow = c(1,1)) plot(seq(1:202),GDI,xlab='Index',ylab=expression(paste(GD[i])),main='p=0.1') abline(h=2*(4+1)/202,lty=2) identify(GDI,n=1) plot(seq(1:202),GDI,xlab='Index',ylab=expression(paste(GD[i])),main='p=0.5') abline(h=2*(4+1)/202,lty=2) identify(GDI,n=1) plot(seq(1:202),GDI,xlab='Index',ylab=expression(paste(GD[i])),main='p=0.9') abline(h=2*(4+1)/202,lty=2) identify(GDI,n=4) ############################################################# ### Graphic of the likelihood displacemente for data(AIS) ### ############################################################# #Dados data(ais, package="sn"); attach(ais); sexInd<-(sex=="female")+0; x=cbind(1,LBM,sexInd); y=BMI #Percentile perc <- 0.9 n <- nrow(x) res <- EM.qr(y,x,perc) thetaest <- res$theta sigmaest <- thetaest[4] betaest <- matrix(thetaest[1:3],3,1) taup2 <- (2/(perc*(1-perc))) thep <- (1-2*peGraphic of the generalized Cook distance for data(AIS)rc)/(perc*(1-perc)) diag <- diag.qr(y,x,perc,thetaest) HessianMatrix <- diag$MatrizQ Gradiente <- diag$mdelta sigma <- sigmaest beta <- betaest muc <- (y-x delta2 <- (y-x gamma2 <- (2+thep^2/taup2)/sigma vchpN <- besselK(sqrt(delta2*gamma2), 0.5-1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2))^(-1) vchp1 <- besselK(sqrt(delta2*gamma2), 0.5+1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2)) Q <- -0.5*n*log(sigmaest)-0.5*(sigmaest*taup2)^{-1}* (sum(vchpN*muc^2 - 2*muc*thep + vchp1*(thep^2+2*taup2))) ######################################################## theta_i <- thetaest sigmaest <- theta_i[4,] betaest <- theta_i[1:3,] sigma <- sigmaest beta <- betaest muc <- (y-x delta2 <- (y-x gamma2 <- (2+thep^2/taup2)/sigma vchpN <- besselK(sqrt(delta2*gamma2), 0.5-1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2))^(-1) vchp1 <- besselK(sqrt(delta2*gamma2), 0.5+1) /(besselK(sqrt(delta2*gamma2), 0.5))*(sqrt(delta2/gamma2)) Q1 <- c() for (i in 1:202) { Q1[i] <- -0.5*n*log(sigmaest[i])-sum(vchpN[,i]*muc[,i]^2 - 2*muc[,i]*thep + vchp1[,i]*(thep^2+2*taup2))/(2*(sigmaest[i]*taup2)) } ######################################################## QDi <- 2*(-Q+Q1) #Depois de escolger perc guardamos os valores de QDi QDi0.1 <- QDi QDi0.5 <- QDi QDi0.9 <- QDi #Seccion de los graficos par(mfrow = c(1,1)) plot(seq(1:202),QDi0.1,xlab='Index',ylab=expression(paste(QD[i])),main='p=0.1') abline(h=mean(QDi0.1)+3.5*sd(QDi0.1),lty=2) identify(QDi0.1,n=3) plot(seq(1:202),QDi0.5,xlab='Index',ylab=expression(paste(QD[i])),main='p=0.5') abline(h=mean(QDi0.5)+3.5*sd(QDi0.5),lty=2) identify(QDi0.5,n=3) plot(seq(1:202),QDi0.9,xlab='Index',ylab=expression(paste(QD[i])),main='p=0.9') abline(h=mean(QDi0.9)+3.5*sd(QDi0.9),lty=2) identify(QDi0.9,n=4) ## End(Not run)
Return estimating the parameters in a quantile regression
EM.qr(y, x = NULL, tau = NULL, error = 0.000001 , iter = 2000, envelope=FALSE)
EM.qr(y, x = NULL, tau = NULL, error = 0.000001 , iter = 2000, envelope=FALSE)
y |
vector of responses |
x |
the design matrix |
tau |
the quantile to be estimated, this is generally a number strictly between 0 and 1. |
error |
the covergence maximum error |
iter |
maximum iterations of the EM algorithm. |
envelope |
confidence envelopes for a curve based on bootstrap replicates |
Estimated parameter for a quantile regression fit,standard error, log-likelihood.
Luis Benites Sanchez [email protected], Christian E. Galarza [email protected], Victor Hugo Lachos [email protected]
[1] Koenker, R. W. (2005). Quantile Regression, Cambridge U. Press.
[2] Yu, K. & Moyeed, R. (2001). Bayesian quantile regression. Statistics & Probability Letters, 54 (4), 437 to 447.
[3] Kotz, S., Kozubowski, T. & Podgorski, K. (2001). The laplace distribution and generalizations: A revisit with applications to communications, economics, engineering, and finance. Number 183. Birkhauser.
data(ais, package="sn") attach(ais) sexInd <- (sex=="female") + 0 x <- cbind(1,LBM,sexInd) y <- BMI tau <- 0.5 ## EM.qr EM.qr(y,x,tau)
data(ais, package="sn") attach(ais) sexInd <- (sex=="female") + 0 x <- cbind(1,LBM,sexInd) y <- BMI tau <- 0.5 ## EM.qr EM.qr(y,x,tau)