LOAD LIBRARIES
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(DHARMa)
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)
## Warning: package 'performance' was built under R version 4.5.2
library(nlme)
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:nlme':
##
## lmList
library(mgcViz)
## Loading required package: mgcv
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
## Loading required package: qgam
## Registered S3 method overwritten by 'mgcViz':
## method from
## +.gg ggplot2
##
## Attaching package: 'mgcViz'
##
## The following objects are masked from 'package:stats':
##
## qqline, qqnorm, qqplot
library(gamlss)
## Loading required package: splines
## Loading required package: gamlss.data
##
## Attaching package: 'gamlss.data'
##
## The following object is masked from 'package:datasets':
##
## sleep
##
## Loading required package: gamlss.dist
## Loading required package: parallel
## ********** GAMLSS Version 5.5-0 **********
## For more on GAMLSS look at https://www.gamlss.com/
## Type gamlssNews() to see new features/changes/bug fixes.
##
##
## Attaching package: 'gamlss'
##
## The following object is masked from 'package:mgcv':
##
## lp
##
## The following object is masked from 'package:lme4':
##
## refit
##
## The following object is masked from 'package:DHARMa':
##
## getQuantile
library(quantreg)
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
##
## The following object is masked from 'package:Matrix':
##
## det
library(corrplot)
## corrplot 0.95 loaded
FIRST LOOK
df <- read.csv("Autism.csv") %>%
as_tibble()
dim(df)
## [1] 704 21
names(df)
## [1] "A1_Score" "A2_Score" "A3_Score" "A4_Score"
## [5] "A5_Score" "A6_Score" "A7_Score" "A8_Score"
## [9] "A9_Score" "A10_Score" "age" "gender"
## [13] "ethnicity" "jundice" "austim" "contry_of_res"
## [17] "used_app_before" "result" "age_desc" "relation"
## [21] "Class.ASD"
head(df)
## # A tibble: 6 × 21
## A1_Score A2_Score A3_Score A4_Score A5_Score A6_Score A7_Score A8_Score
## <int> <int> <int> <int> <int> <int> <int> <int>
## 1 1 1 1 1 0 0 1 1
## 2 1 1 0 1 0 0 0 1
## 3 1 1 0 1 1 0 1 1
## 4 1 1 0 1 0 0 1 1
## 5 1 0 0 0 0 0 0 1
## 6 1 1 1 1 1 0 1 1
## # ℹ 13 more variables: A9_Score <int>, A10_Score <int>, age <chr>,
## # gender <chr>, ethnicity <chr>, jundice <chr>, austim <chr>,
## # contry_of_res <chr>, used_app_before <chr>, result <int>, age_desc <chr>,
## # relation <chr>, Class.ASD <chr>
str(df)
## tibble [704 × 21] (S3: tbl_df/tbl/data.frame)
## $ A1_Score : int [1:704] 1 1 1 1 1 1 0 1 1 1 ...
## $ A2_Score : int [1:704] 1 1 1 1 0 1 1 1 1 1 ...
## $ A3_Score : int [1:704] 1 0 0 0 0 1 0 1 0 1 ...
## $ A4_Score : int [1:704] 1 1 1 1 0 1 0 1 0 1 ...
## $ A5_Score : int [1:704] 0 0 1 0 0 1 0 0 1 0 ...
## $ A6_Score : int [1:704] 0 0 0 0 0 0 0 0 0 1 ...
## $ A7_Score : int [1:704] 1 0 1 1 0 1 0 0 0 1 ...
## $ A8_Score : int [1:704] 1 1 1 1 1 1 1 0 1 1 ...
## $ A9_Score : int [1:704] 0 0 1 0 0 1 0 1 1 1 ...
## $ A10_Score : int [1:704] 0 1 1 1 0 1 0 0 1 0 ...
## $ age : chr [1:704] "26" "24" "27" "35" ...
## $ gender : chr [1:704] "f" "m" "m" "f" ...
## $ ethnicity : chr [1:704] "White-European" "Latino" "Latino" "White-European" ...
## $ jundice : chr [1:704] "no" "no" "yes" "no" ...
## $ austim : chr [1:704] "no" "yes" "yes" "yes" ...
## $ contry_of_res : chr [1:704] "'United States'" "Brazil" "Spain" "'United States'" ...
## $ used_app_before: chr [1:704] "no" "no" "no" "no" ...
## $ result : int [1:704] 6 5 8 6 2 9 2 5 6 8 ...
## $ age_desc : chr [1:704] "'18 and more'" "'18 and more'" "'18 and more'" "'18 and more'" ...
## $ relation : chr [1:704] "Self" "Self" "Parent" "Self" ...
## $ Class.ASD : chr [1:704] "NO" "NO" "YES" "NO" ...
CLEAN DATA
df[df == "?"] <- NA
colSums(is.na(df))
## A1_Score A2_Score A3_Score A4_Score A5_Score
## 0 0 0 0 0
## A6_Score A7_Score A8_Score A9_Score A10_Score
## 0 0 0 0 0
## age gender ethnicity jundice austim
## 2 0 95 0 0
## contry_of_res used_app_before result age_desc relation
## 0 0 0 0 95
## Class.ASD
## 0
colMeans(is.na(df)) * 100
## A1_Score A2_Score A3_Score A4_Score A5_Score
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## A6_Score A7_Score A8_Score A9_Score A10_Score
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## age gender ethnicity jundice austim
## 0.2840909 0.0000000 13.4943182 0.0000000 0.0000000
## contry_of_res used_app_before result age_desc relation
## 0.0000000 0.0000000 0.0000000 0.0000000 13.4943182
## Class.ASD
## 0.0000000
df <- df %>% drop_na()
df$asd <- ifelse(df$Class.ASD == "YES", 1, 0)
df <- df %>%
mutate(
across(A1_Score:A10_Score, as.integer),
age = as.numeric(age),
result = as.numeric(result),
gender = factor(gender),
ethnicity = factor(ethnicity),
jundice = factor(jundice),
austim = factor(austim),
used_app_before = factor(used_app_before),
relation = factor(relation),
Class.ASD = factor(Class.ASD, levels = c("NO", "YES"))
)
str(df)
## tibble [609 × 22] (S3: tbl_df/tbl/data.frame)
## $ A1_Score : int [1:609] 1 1 1 1 1 0 1 1 1 1 ...
## $ A2_Score : int [1:609] 1 1 1 1 1 1 1 1 1 1 ...
## $ A3_Score : int [1:609] 1 0 0 0 1 0 1 0 1 1 ...
## $ A4_Score : int [1:609] 1 1 1 1 1 0 1 0 1 1 ...
## $ A5_Score : int [1:609] 0 0 1 0 1 0 0 1 0 1 ...
## $ A6_Score : int [1:609] 0 0 0 0 0 0 0 0 1 1 ...
## $ A7_Score : int [1:609] 1 0 1 1 1 0 0 0 1 1 ...
## $ A8_Score : int [1:609] 1 1 1 1 1 1 0 1 1 1 ...
## $ A9_Score : int [1:609] 0 0 1 0 1 0 1 1 1 1 ...
## $ A10_Score : int [1:609] 0 1 1 1 1 0 0 1 0 1 ...
## $ age : num [1:609] 26 24 27 35 36 17 64 29 17 33 ...
## $ gender : Factor w/ 2 levels "f","m": 1 2 2 1 2 1 2 2 2 2 ...
## $ ethnicity : Factor w/ 11 levels "'Middle Eastern '",..: 11 6 6 11 8 4 11 11 3 11 ...
## $ jundice : Factor w/ 2 levels "no","yes": 1 1 2 1 2 1 1 1 2 1 ...
## $ austim : Factor w/ 2 levels "no","yes": 1 2 2 2 1 1 1 1 2 1 ...
## $ contry_of_res : chr [1:609] "'United States'" "Brazil" "Spain" "'United States'" ...
## $ used_app_before: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ result : num [1:609] 6 5 8 6 9 2 5 6 8 10 ...
## $ age_desc : chr [1:609] "'18 and more'" "'18 and more'" "'18 and more'" "'18 and more'" ...
## $ relation : Factor w/ 5 levels "'Health care professional'",..: 5 5 3 5 5 5 3 5 1 4 ...
## $ Class.ASD : Factor w/ 2 levels "NO","YES": 1 1 2 1 2 1 1 1 2 2 ...
## $ asd : num [1:609] 0 0 1 0 1 0 0 0 1 1 ...
EXPLORATORY DATA ANALYSIS
ggplot(df, aes(x = Class.ASD, fill = Class.ASD)) +
geom_bar() +
scale_fill_manual(values = c("violet", "brown")) +
labs(title = "ASD Diagnosis Distribution",
x = "Diagnosis", y = "Count") +
theme_minimal()
hist(df$result,
main = "Distribution of Screening Score",
xlab = "Score (0-10)",
col = "lightblue",
breaks = 11)
mean.score <- mean(df$result)
abline(v = mean.score, lwd = 3, col = "red")
mean.score
## [1] 5.077176
ggplot(df, aes(x = age, fill = Class.ASD)) +
geom_density(alpha = 0.5) +
scale_fill_manual(values = c("blue", "orange")) +
labs(title = "Age Distribution by ASD Diagnosis") +
theme_minimal()
# Gender
ggplot(df, aes(x = gender, fill = Class.ASD)) +
geom_bar(position = "fill") +
scale_fill_manual(values = c("midnightblue", "lightblue")) +
labs(title = "ASD Proportion by Gender",
y = "Proportion") +
theme_minimal()
# Jaundice
ggplot(df, aes(x = jundice, fill = Class.ASD)) +
geom_bar(position = "fill") +
scale_fill_manual(values = c("skyblue", "darkgreen")) +
labs(title = "ASD Proportion by Jaundice at Birth",
y = "Proportion") +
theme_minimal()
# Family history of ASD
ggplot(df, aes(x = austim, fill = Class.ASD)) +
geom_bar(position = "fill") +
scale_fill_manual(values = c("gray", "yellow")) +
labs(title = "ASD Proportion by Family History of ASD",
y = "Proportion") +
theme_minimal()
COLLINEARITY CHECK
items <- df %>% select(A1_Score:A10_Score)
cor_matrix <- cor(items)
corrplot(cor_matrix,
method = "color",
type = "upper",
tl.cex = 0.8,
title = "Correlation Among Screening Items")
Fit Baseline Model
fit.1 <- glm(
formula = asd ~ 1,
data = df,
family = binomial()
)
summary(fit.1)
##
## Call:
## glm(formula = asd ~ 1, family = binomial(), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.86850 0.08881 -9.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 739.4 on 608 degrees of freedom
## Residual deviance: 739.4 on 608 degrees of freedom
## AIC: 741.4
##
## Number of Fisher Scoring iterations: 4
exp(coef(fit.1))
## (Intercept)
## 0.4195804
exp(confint(fit.1))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## 0.3517823 0.4983783
Single Predictor Model
fit.2 <- glm(
formula = asd ~ result,
data = df,
family = binomial()
)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(fit.2)
##
## Call:
## glm(formula = asd ~ result, family = binomial(), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -289.37 49899.89 -0.006 0.995
## result 44.52 7700.48 0.006 0.995
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7.3940e+02 on 608 degrees of freedom
## Residual deviance: 4.9937e-08 on 607 degrees of freedom
## AIC: 4
##
## Number of Fisher Scoring iterations: 25
exp(coef(fit.2))
## (Intercept) result
## 2.129857e-126 2.155566e+19
exp(confint(fit.2))
## Waiting for profiling to be done...
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## 2.5 % 97.5 %
## (Intercept) 0 0
## result Inf Inf
sim.2 <- simulateResiduals(fit.2)
plot(sim.2)
Multiple Predictors
fit.3 <- glm(
formula = asd ~ result + age + gender + jundice + austim,
data = df,
family = binomial()
)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(fit.3)
##
## Call:
## glm(formula = asd ~ result + age + gender + jundice + austim,
## family = binomial(), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.894e+02 5.172e+04 -0.006 0.996
## result 4.449e+01 7.710e+03 0.006 0.995
## age 1.147e-02 4.031e+02 0.000 1.000
## genderm -1.893e-01 7.852e+03 0.000 1.000
## jundiceyes 1.271e-01 1.334e+04 0.000 1.000
## austimyes -3.113e-01 1.255e+04 0.000 1.000
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7.3940e+02 on 608 degrees of freedom
## Residual deviance: 4.9917e-08 on 603 degrees of freedom
## AIC: 12
##
## Number of Fisher Scoring iterations: 25
exp(coef(fit.3))
## (Intercept) result age genderm jundiceyes
## 2.013599e-126 2.100232e+19 1.011539e+00 8.275467e-01 1.135522e+00
## austimyes
## 7.324743e-01
sim.3 <- simulateResiduals(fit.3)
plot(sim.3)
r2_mcfadden(fit.3)
## # R2 for Generalized Linear Regression
## R2: 1.000
## adj. R2: 0.997
GENERALIZED LINEAR MIXED MODEL
fit.glmm <- glmer(
formula = asd ~ result + age + gender + jundice + austim + (1 | ethnicity),
data = df,
family = binomial()
)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determined
summary(fit.glmm)
## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX
## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: asd ~ result + age + gender + jundice + austim + (1 | ethnicity)
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 14.0 44.9 0.0 0.0 602
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.49e-08 -1.49e-08 -1.49e-08 1.49e-08 1.49e-08
##
## Random effects:
## Groups Name Variance Std.Dev.
## ethnicity (Intercept) 1.162e-05 0.003408
## Number of obs: 609, groups: ethnicity, 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.840e+02 8.462e+06 0 1
## result 7.445e+01 1.099e+06 0 1
## age -2.043e-02 1.580e+05 0 1
## genderm -4.511e-01 5.481e+06 0 1
## jundiceyes 1.459e+00 9.347e+06 0 1
## austimyes 4.040e-01 8.122e+06 0 1
##
## Correlation of Fixed Effects:
## (Intr) result age gendrm jndcys
## result -0.650
## age -0.581 0.019
## genderm -0.403 0.026 0.052
## jundiceyes -0.010 -0.081 -0.046 0.002
## austimyes -0.003 -0.162 -0.077 0.101 -0.136
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Hessian is numerically singular: parameters are not uniquely determined
exp(fixef(fit.glmm))
## (Intercept) result age genderm jundiceyes
## 6.568994e-211 2.164385e+32 9.797808e-01 6.369536e-01 4.302814e+00
## austimyes
## 1.497844e+00
ranef(fit.glmm)
## $ethnicity
## (Intercept)
## 'Middle Eastern ' -1.960352e-19
## 'South Asian' -7.738994e-20
## Asian -2.347495e-19
## Black -1.699684e-20
## Hispanic -7.738994e-21
## Latino 0.000000e+00
## others -2.579665e-21
## Others -3.095597e-20
## Pasifika -2.579665e-20
## Turkish -1.031866e-20
## White-European -3.707368e-20
##
## with conditional variances for "ethnicity"
fit.gam <- gam(
formula = asd ~ s(age) + s(result) + gender + jundice + austim,
family = binomial(),
data = df
)
summary(fit.gam)
##
## Family: binomial
## Link function: logit
##
## Formula:
## asd ~ s(age) + s(result) + gender + jundice + austim
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.282e+01 3.839e+05 0 1
## genderm -1.990e-01 2.553e+05 0 1
## jundiceyes 1.451e-01 4.349e+05 0 1
## austimyes -2.786e-01 4.085e+05 0 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(age) 1 1 0 1.0
## s(result) 1 1 0 0.5
##
## R-sq.(adj) = 1 Deviance explained = 100%
## UBRE = -0.9803 Scale est. = 1 n = 609
plot(fit.gam, pages = 1, shade = TRUE, shade.col = "lightblue",
main = "GAM Smooth Effects")
sim.gam <- simulateResiduals(fit.gam)
plot(sim.gam)
Question: does family history matter more for people who already
score high?
q <- seq(.1,.9,.1)
q
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
fit.qr <- rq((result ~ age + gender + jundice + austim),
data = df,
tau = q
)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
summary(fit.qr,se = "nid")
## Warning in summary.rq(xi, U = U, ...): 4 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 1 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 2 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 1 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 17 non-positive fis
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.1
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 2.06389 0.24723 8.34805 0.00000
## age -0.00278 0.01257 -0.22106 0.82512
## genderm 0.00000 0.10006 0.00000 1.00000
## jundiceyes -0.92778 0.55557 -1.66996 0.09545
## austimyes 1.02222 0.14273 7.16200 0.00000
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.2
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 3.10440 0.36613 8.47902 0.00000
## age -0.00549 0.01086 -0.50606 0.61300
## genderm 0.02198 0.24214 0.09076 0.92771
## jundiceyes 0.01648 0.41208 0.04000 0.96811
## austimyes 1.02747 0.35582 2.88763 0.00402
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.3
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 3.30030 0.38711 8.52556 0.00000
## age -0.00601 0.01212 -0.49566 0.62032
## genderm 0.00601 0.22686 0.02648 0.97889
## jundiceyes 0.87988 0.68163 1.29084 0.19725
## austimyes 0.92192 0.34963 2.63684 0.00858
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.4
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 4.00000 0.65924 6.06762 0.00000
## age 0.00000 0.02200 0.00000 1.00000
## genderm 0.00000 0.23679 0.00000 1.00000
## jundiceyes 1.00000 0.62736 1.59398 0.11146
## austimyes 1.00000 0.54566 1.83263 0.06735
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 4.00000 0.60806 6.57831 0.00000
## age 0.00000 0.02149 0.00000 1.00000
## genderm 0.00000 0.24654 0.00000 1.00000
## jundiceyes 2.00000 0.63195 3.16481 0.00163
## austimyes 2.00000 0.51905 3.85317 0.00013
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.6
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 5.00000 0.78532 6.36681 0.00000
## age 0.00000 0.02599 0.00000 1.00000
## genderm 0.00000 0.36595 0.00000 1.00000
## jundiceyes 1.00000 0.44122 2.26647 0.02378
## austimyes 2.00000 0.46702 4.28247 0.00002
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.7
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.00000 0.78540 7.63940 0.00000
## age 0.00000 0.02387 0.00000 1.00000
## genderm 0.00000 0.29666 0.00000 1.00000
## jundiceyes 1.00000 0.30855 3.24097 0.00126
## austimyes 2.00000 0.32789 6.09961 0.00000
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.8
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 7.00000 0.67836 10.31893 0.00000
## age 0.00000 0.02205 0.00000 1.00000
## genderm 0.00000 0.31403 0.00000 1.00000
## jundiceyes 2.00000 0.29803 6.71068 0.00000
## austimyes 1.00000 0.32296 3.09632 0.00205
##
## Call: rq(formula = (result ~ age + gender + jundice + austim), tau = q,
## data = df)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 7.13043 0.57605 12.37821 0.00000
## age 0.04348 0.01527 2.84707 0.00456
## genderm 0.00000 0.23884 0.00000 1.00000
## jundiceyes 0.47826 0.30708 1.55742 0.11989
## austimyes 0.65217 0.39141 1.66621 0.09619
plot(fit.qr, ols = F)
fit.all.q <- rq(result ~ age + gender + jundice + austim,
q, data = df)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
plot(summary(fit.all.q),
main = "Quantile Regression Coefficients Across Quantiles")
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
