Log.lqr: Robust Logistic Linear Quantile Regression

Description

It performs the logistic transformation in Galarza et.al.(2020) (see references) for estimating quantiles for a bounded response. Once the response is transformed, it uses the lqr function.

Usage

Log.lqr(y,x,p=0.5,a=0,b=1,dist = "normal",nu="",gama="",precision = 10^-6,
epsilon = 0.001,CI=0.95)

Arguments

lower bound for the response (default = 0)

upper bound for the response (default = 1)

epsilon

a small quantity $\epsilon>0$ that ensures that the logistic transform is defined for all values of y

the response vector of dimension $n$ where $n$ is the total of observations.

design matrix for the fixed effects of dimension $N x d$ where $d$ represents the number of fixed effects including the intercept, if considered.

An unique quantile or a set of quantiles related to the quantile regression.

dist

represents the distribution to be used for the error term. The values are normal for Normal distribution, t for Student's t distribution, laplace for Laplace distribution, slash for Slash distribution and cont for the Contaminated normal distribution.

It represents the degrees of freedom when dist = t. For the Slash distribution (dist = slash) it is a shape parameter $\nu>0$. For the Contaminated Normal distribution, $\nu$ is the parameter that represents the percentage of outliers. When is not provided, we use the MLE.

gama

It represents a scale factor for the contaminated normal distribution. When is not provided, we use the MLE.

precision

The convergence maximum error permitted. By default is 10^-6.

Confidence to be used for the Confidence Interval when a grid of quantiles is provided. Default = 0.95.

Value

iter

number of iterations.

criteria

attained criteria value.

beta

fixed effects estimates.

sigma

scale parameter estimate for the error term.

Estimate of nu parameter detailed above.

gamma

Estimate of gamma parameter detailed above.

Standard Error estimates.

table

Table containing the inference for the fixed effects parameters.

loglik

Log-likelihood value.

AIC

Akaike information criterion.

BIC

Bayesian information criterion.

Hannan-Quinn information criterion.

fitted.values

vector containing the fitted values.

residuals

vector containing the residuals.

Details

We follow the transformation in Bottai et.al. (2009) defined as

$$h(y)=logit(y)=log(\frac{y-a}{b-y})$$

that implies

$$Q_{y}(p)=\frac{b\,exp(X\beta) + a}{1 + exp(X\beta)}$$

where $Q_{y}(p)$ represents the conditional quantile of the response. Once estimates for the regression coefficients $\beta_p$ are obtained, inference on $Q_{y}(p)$ can then be made through the inverse transform above. This equation (as function) is provided in the output. See example.

The interpretation of the regression coefficients is analogous to the interpretation of the coefficients of a logistic regression for binary outcomes.

For example, let $x_1$ be the gender (male = 0, female=1). Then $exp(\beta_{0.5,1})$ represents the odds ratio of median score in males vs females, where the odds are defined using the score instead of a probability, $(y-a)/(b-y)$. When the covariate is continous, the respective $\beta$ coeficient can be interpretated as the increment (or decrement) over the log(odd ratio) when the covariate increases one unit.

References

Galarza, C.M., Zhang P. and Lachos, V.H. (2020). Logistic Quantile Regression for Bounded Outcomes Using a Family of Heavy-Tailed Distributions. Sankhya B: The Indian Journal of Statistics. 10.1007/s13571-020-00231-0

Galarza, C., Lachos, V. H., Cabral, C. R. B., & Castro, C. L. (2017). Robust quantile regression using a generalized class of skewed distributions. Stat, 6(1), 113-130.

Examples

Run this code

# NOT RUN {
##Load the data
data(resistance)
attach(resistance)

#EXAMPLE 1.1

#Comparing the resistence to death of two types of tumor-cells.
#The response is a score in [0,4].

boxplot(score~type,ylab="score",xlab="type")

#Student't median logistic quantile regression
res = Log.lqr(y = score,x = cbind(1,type),a=0,b=4,dist="t")

# The odds ratio of median score in type B vs type A
exp(res$beta[2])

#Proving that exp(res$beta[2])  is approx median odd ratio
medA  = median(score[type=="A"])
medB  = median(score[type=="B"])
rateA = (medA - 0)/(4 - medA)
rateB = (medB - 0)/(4 - medB)
odd   = rateB/rateA

round(c(exp(res$beta[2]),odd),3)


#EXAMPLE 1.2
############

#Comparing the resistence to death depending of dose.

#descriptive
plot(dose,score,ylim=c(0,4),col="dark gray");abline(h=c(0,4),lty=2)
dosecat<-cut(dose, 6, ordered = TRUE)
boxplot(score~dosecat,ylim=c(0,4))
abline(h=c(0,4),lty=2)

#Slash (Non logistic) quantile regression for quantiles 0.05, 0.50 and 0.95
xx1  = dose
xx2  = dose^2
xx3  = dose^3
res3 = lqr(y = score,x = cbind(1,xx1,xx2,xx3),p = c(0.05,0.50,0.95),dist="slash")
seqq=seq(min(dose),max(dose),length.out = 1000)
dd = matrix(data = NA,nrow = 1000,ncol =3)
for(i in 1:3)
{
  dd[,i] = rep(res3[[i]]$beta[1],1000) + res3[[i]]$beta[2]*seqq +
    res3[[i]]$beta[3]*seqq^2 + res3[[i]]$beta[4]*seqq^3
}

plot(dose,score,ylim=c(-1,5),col="gray");abline(h=c(0,4),lty=2)
lines(seqq,dd[,1],lwd=1,col=2)
lines(seqq,dd[,2],lwd=1,col=1)
lines(seqq,dd[,3],lwd=1,col=2)

#Using logistic quantile regression for obtaining predictions inside bounds

res4 = Log.lqr(y = score,x = cbind(1,xx1,xx2,xx3),a = 0,b = 4,p = c(0.05,0.50,0.95),dist="slash")
dd = matrix(data = NA,nrow = 1000,ncol =3)
for(i in 1:3)
{
  dd[,i] = rep(res4[[i]]$beta[1],1000) + res4[[i]]$beta[2]*seqq +
    res4[[i]]$beta[3]*seqq^2 + res4[[i]]$beta[4]*seqq^3
}

#Computing quantiles for the original response (Inverse trnasformation)

pred = function(predlog,a,b)
{
  return((b*exp(predlog)+a)/(1+exp(predlog)))
}

for(i in 1:3)
{
  dd[,i] = pred(dd[,i],a=0,b=4)
}           

#No more prediction curves outof bounds
plot(dose,score,ylim=c(0,4),col="gray");abline(h=c(0,4),lty=2)
lines(seqq,dd[,1],lwd=1,col=2)
lines(seqq,dd[,2],lwd=1,col=1)
lines(seqq,dd[,3],lwd=1,col=2)

#EXAMPLE 1.3
############

#A full model using dose and type for a grid of quantiles

typeB = 1*(type=="B")
res5 = Log.lqr(y = score,x = cbind(1,xx1,xx2,xx3,typeB,typeB*xx1),a = 0,b = 4,
               p = seq(from = 0.05,to = 0.95,by = 0.05),dist = "t")
ddA = ddB = matrix(data = NA,nrow = 1000,ncol = 5)
for(i in 1:5)
{
  k = c(2,5,10,15,18)[i]
  ddA[,i] = rep(res5[[k]]$beta[1],1000) + res5[[k]]$beta[2]*seqq + res5[[k]]$beta[3]*
    seqq^2 + res5[[k]]$beta[4]*seqq^3
  ddB[,i] = rep(res5[[k]]$beta[1],1000) + (res5[[k]]$beta[2] + res5[[k]]$beta[6])*
    seqq + res5[[k]]$beta[3]*seqq^2 + res5[[k]]$beta[4]*seqq^3  + res5[[k]]$beta[5]
}

#Computing quantiles for the original response (Inverse transformation)

for(i in 1:5)
{
  ddA[,i] = pred(ddA[,i],a=0,b=4)
  ddB[,i] = pred(ddB[,i],a=0,b=4)
} 

#Such a beautiful plot
par(mfrow=c(1,2))
plot(dose,score,ylim=c(0,4),col=c((type == "B")*8+(type == "A")*1),main="Type A")
abline(h=c(0,4),lty=2)
lines(seqq,ddA[,1],lwd=2,col=2)
lines(seqq,ddA[,2],lwd=1,col=4)
lines(seqq,ddA[,3],lwd=2,col=1)
lines(seqq,ddA[,4],lwd=1,col=4)
lines(seqq,ddA[,5],lwd=2,col=2)

legend(x = 0,y=4,legend = c("p=0.10","p=0.25","p=0.50","p=0.75","p=0.90")
       ,col=c(2,4,1,4,2),lwd=c(2,1,2,1,2),bty = "n",cex=0.65)

plot(dose,score,ylim=c(0,4),col=c((type == "B")*1 + (type == "A")*8),
     main="Type B");abline(h=c(0,4),lty=2)
lines(seqq,ddB[,1],lwd=2,col=2)
lines(seqq,ddB[,2],lwd=1,col=4)
lines(seqq,ddB[,3],lwd=2,col=1)
lines(seqq,ddB[,4],lwd=1,col=4)
lines(seqq,ddB[,5],lwd=2,col=2)
# }

Run the code above in your browser using DataLab