diag.qr: Diagnostics for Quantile Regression Using Asymmetric Laplace...

Description Usage Arguments Value Author(s) References Examples

View source: R/diag.qr.r

Description

Return case-deletion estimating the parameters in a quantile regression

Usage

1
  diag.qr(y,x,tau,theta)

Arguments

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

Value

Hessian and gradient matrix. Also the generalized cook distance (GDi), approximation of the likelihood distance (QDi)

Author(s)

Luis Benites Sanchez [email protected], Christian E. Galarza [email protected], Victor Hugo Lachos [email protected]

References

[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.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
## 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) 

ALDqr documentation built on May 30, 2017, 8:26 a.m.