Quantile regression in R

Hi everyone,
I have a question :). I am trying to do a quantile regression, comparing median number of counselling sessions (no_couns_sess) between different age categories (age_cat3_l). The data set I am using is called MH_Data_Clean2.
The code I am using from the quantreg package:

model1 <- rq(no_couns_sess ~ age_cat3_l, data = MH_Data_Clean2, tau = 0.5)
summary (model1)

When I run this, I only see:

model1 ← rq(no_couns_sess ~ age_cat3_l, data = MH_Data_Clean2, tau = 0.5)
summary (model1)

Call: rq(formula = no_couns_sess ~ age_cat3_l, tau = 0.5, data = MH_Data_Clean2)

tau: [1] 0.5

Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 1.00000 0.11900 8.40312 0.00000
age_cat3_l 0.00000 0.05356 0.00000 1.00000

Question: how can I get the output to also include 95% confidence intervals?

Many thanks!

Thanks for your post @elburgvb.
Perhaps @cmaronga @mumbua_mutunga @arranh or @cjarvis11 can assist?

Hi,

You could use tbl_regression() from the gtsummary package to display your 95% CI’s in a table.

I’ve written some code to demonstrate this, with the relevant packages. Let me know if you run into difficulties extrapolating from this to your example.

#Install pacman if missing
if(!require(pacman)) install.packages("pacman")

#Load packages
pacman::p_load(
  quantreg,
  janitor,
  broom,
  gtsummary,
  tidyverse
)

#Load data
data("iris")

#Set up regression
model <- iris %>%
  clean_names() %>%
  rq(formula = sepal_length ~ petal_length + sepal_width,
            tau = 0.5)

#Display in gtsummary
tbl_regression(model)

Hi Arran,

Thanks so much for this! I am indeed still running into issues here…
``
#Set up regression
model1 ← MH_Data_Clean2 %>%
clean_names()%>%
rq(formula = no_couns_sess ~ age_cat3_l, tau = 0.5)

#Display in gtsummary
tbl_regression(model1)`
``

This is the code I used.

But now this is the error I get:
! broom::tidy() failed to tidy the model.
! tidy_parameters() also failed.
:heavy_multiplication_x: Sorry, model_parameters() failed with the following error (possible class rq not supported):

singular matrix in ‘backsolve’. First zero in diagonal [1]
! broom::tidy() failed to tidy the model.
! tidy_parameters() also failed.
:heavy_multiplication_x: Sorry, model_parameters() failed with the following error (possible class rq not supported):

singular matrix in ‘backsolve’. First zero in diagonal [1]
:heavy_multiplication_x: There was an error calling tidy_fun(). Most likely, this is because the
function supplied in tidy_fun= was misspelled, does not exist, is not
compatible with your object, or was missing necessary arguments (e.g. conf.level= or conf.int=). See error message below.
Error:
! Error: ! Unable to tidy x.
Run rlang::last_trace() to see where the error occurred.
Warning messages:
1: The exponentiate argument is not supported in the tidy() method for rq objects and will be ignored.
2: In summary.rq(x) : 9028 non-positive fis
3: In summary.rq(x, covariance = TRUE) : 9028 non-positive fis
4: In summary.rq(x) : 9028 non-positive fis
5: In summary.rq(x, covariance = TRUE) : 9028 non-positive fis

I have installed all the packages you mentioned. I checked the class of both variables (no_couns_sess and age_cat3_l, they are both numeric)

This is what the frequency table of no_couns_sess (number of counselling sessions) looks like:
no_couns_sess n percent valid_percent
1 5227 5.032736e-01 0.5789765175
2 1070 1.030233e-01 0.1185201595
3 924 8.896592e-02 0.1023482499
4 679 6.537647e-02 0.0752104564
5 561 5.401502e-02 0.0621400089
6 233 2.243405e-02 0.0258085955
7 135 1.299827e-02 0.0149534781
8 81 7.798960e-03 0.0089720868
9 42 4.043905e-03 0.0046521932
10 32 3.081071e-03 0.0035445281
11 16 1.540535e-03 0.0017722641
12 9 8.665511e-04 0.0009968985
13 8 7.702677e-04 0.0008861320
14 4 3.851338e-04 0.0004430660
15 2 1.925669e-04 0.0002215330
16 1 9.628346e-05 0.0001107665
18 1 9.628346e-05 0.0001107665
22 1 9.628346e-05 0.0001107665
24 1 9.628346e-05 0.0001107665
50 1 9.628346e-05 0.0001107665
NA 1358 1.307529e-01 NA

I am grateful for any advice you can give!!!
Warm wishes,
Elburg

Does the code that I put in the thread work for you?

It should be creating a table that appears in the plot window.

Hi Arran,

Yes, your code works perfectly.

When I adapt your code to my dataset/variables looking at number of counselling sessions and age groups, it runs smoothly as well:

#Table 5: Compare distribution of total number of counselling sessions between groups: age categories
#Set up regression
model1 ← MH_Data_Clean2 %>%
clean_names()%>%
rq(formula = no_couns_sess ~ age_cat3_l, tau = 0.5)

#Display in gtsummary
tbl_regression(model1)

However, when I run the same code but with the variable ‘sex’ instead of age group, I get an error. Code:
#Table 5: Compare distribution of total number of counselling sessions between groups: sex
#Set up regression
model2 ← MH_Data_Clean2 %>%
clean_names()%>%
rq(formula = no_couns_sess ~ sex,
tau = 0.5)

#Display in gtsummary
tbl_regression(model2)

Error:

#Table 5: Compare distribution of total number of counselling sessions between groups: sex
#Set up regression
model2 ← MH_Data_Clean2 %>%

  • clean_names()%>%
  • rq(formula = no_couns_sess ~ sex,
  •  tau = 0.5)
    

Warning message:
In rq.fit.br(x, y, tau = tau, …) : Solution may be nonunique

#Display in gtsummary
tbl_regression(model2)
! broom::tidy() failed to tidy the model.
! tidy_parameters() also failed.
:heavy_multiplication_x: Sorry, model_parameters() failed with the following error (possible class rq not supported):

singular matrix in ‘backsolve’. First zero in diagonal [1]
! broom::tidy() failed to tidy the model.
! tidy_parameters() also failed.
:heavy_multiplication_x: Sorry, model_parameters() failed with the following error (possible class rq not supported):

singular matrix in ‘backsolve’. First zero in diagonal [1]
:heavy_multiplication_x: There was an error calling tidy_fun(). Most likely, this is because the
function supplied in tidy_fun= was misspelled, does not exist, is not
compatible with your object, or was missing necessary arguments (e.g. conf.level= or conf.int=). See error message below.
Error:
! Error: ! Unable to tidy x.
Run rlang::last_trace() to see where the error occurred.
Warning messages:
1: The exponentiate argument is not supported in the tidy() method for rq objects and will be ignored.
2: In summary.rq(x) : 9544 non-positive fis
3: In summary.rq(x, covariance = TRUE) : 9544 non-positive fis
4: In summary.rq(x) : 9544 non-positive fis
5: In summary.rq(x, covariance = TRUE) : 9544 non-positive fis

I think your error is to do with your variable ‘sex’ in the regression.

It could be because one of variables in ‘sex’ is exclusively (or almost exclusively) associated with your outcome (concurvity {mgcv} - unknown error - General - Posit Community), or ‘sex’ is very highly correlated with another variable (https://www.reddit.com/r/rstats/comments/p7hxnn/quantile_regression_backsolver_errors/).

Could you carry out a correlation on your variables, and also look at splitting the dataset by ‘sex’ to see if there are any strong associations?

I’ve added some example code below, with the same iris dataset

#Install pacman if missing
if(!require(pacman)) install.packages("pacman")

#Load packages
pacman::p_load(
  quantreg,
  janitor,
  broom,
  gtsummary,
  
  corrr,   #New package added
  
  tidyverse
)

#Load data
data("iris")

#Look at the correlations
iris %>%
  corrr::correlate()

#Stratify by species 
iris %>%
  tbl_summary(
    by = Species
  )

What is the output of you running the below code?

MH_Data_Clean2 %>%
  clean_names() %>%
  select( **include here all variables you are interested in** ) %>%
  tbl_summary(
    by = sex
  )

Sorry for the late reply :slight_smile:


The only thing I can think of here is that the values of ‘no_couns_sess’ looks the same for both values of ‘sex’. To investigate that further can you plot what the distribution of ‘no_couns_sess’ by ‘sex’?

ggplot(data = MH_Data_Clean2 ,
       mapping = aes(
         x = sex,
         y = no_couns_sess
       )) +
  geom_boxplot() +
  theme_bw()

Hi Arran!

Apologies for the delayed reply!

I tried it out and indeed the values look the same for both groups, so that must have been the problem.

Thank you so much for your patience and support in figuring this out and I apologize for the rabbit hole I pulled you into!!!

Warm wishes,
Elburg

1 Like

Glad I could help!