#Objective: This study aimed to investigate the impact of insurance status on wait times for otolaryngology care, comparing Medicaid-insured patients to commercially insured patients.
#
# Study Design: The study utilized an audit methodology, known as a "mystery caller" study, to assess appointment availability and patient experiences regarding access to care in otolaryngology.
#
# Setting: The study included physicians representing various otolaryngology subspecialties across the United States, excluding military medical practices.
#
# Methods: Physicians were selected from patient-facing directories and stratified by region. Mystery callers, posing as patients with either Medicaid or commercial insurance, made two separate calls to each physician's office to obtain the earliest possible appointment time. The calls were standardized and completed within one week. Data on appointment availability, earliest appointment dates, and additional information were collected using a secure electronic data capture tool.
#
# Results: Out of 612 physicians contacted, 301 physicians accepting new patients were included in the analysis. The median wait time across all subspecialties and insurance types was 34.3 business days. The study found a statistically significant difference in wait times based on insurance type, with Medicaid-insured patients experiencing a 13.8% longer wait time than commercially insured patients. The model-estimated average wait times were 32.4 days for the commercially insured group and 36.8 days for the Medicaid group.
#
# Conclusions: The study revealed that Medicaid-insured patients in otolaryngology care faced longer wait times compared to commercially insured patients. These findings contribute to the existing literature on access to care and highlight the need to address potential disparities in wait times to promote equitable access to otolaryngology services.| npi | name | N |
|---|---|---|
## CSV file saved successfully!| npi | name | Reason for exclusions | insurance | business_days_until_appointment |
|---|---|---|---|---|
## [1] 4 5 28 80 130 166 205 272 297 320 343 345 396 412 495 521 550 586 603
## [20] 668 691 815 915 942| x |
|---|
| 4 |
| 5 |
| 28 |
| 80 |
| 130 |
| 166 |
| 205 |
| 272 |
| 297 |
| 320 |
| 343 |
| 345 |
| 396 |
| 412 |
| 495 |
| 521 |
| 550 |
| 586 |
| 603 |
| 668 |
| 691 |
| 815 |
| 915 |
| 942 |
## # A tibble: 0 × 2
## # ℹ 2 variables: npi <dbl>, calls_count <int>| npi | calls_count |
|---|---|
A total of 379 unique orthopedic surgeons were identified in the dataset and were successfully contacted (i.e., with a recorded wait time for an appointment) in 47 states including the District of Columbia. The excluded states include North Dakota, West Virginia, Wyoming and District of Columbia.
| median_wait_time | Q1 | Q3 |
|---|---|---|
| 13 | 6 | 31 |
The median wait time across all subspecialties and insurance types was 13 business days, with an interquartile range (IQR) of 6 to 31.
Each physician received 2, phone calls with identical clinical scenarios. In this process, 4 ” physicians were excluded after two unsuccessful attempts to reach them. Out of the 978 phone calls made, 34 calls 3.4764826%) were on hold for more than five minutes, 75 7.6687117 %) went to voicemail, 130 13.2924335%) physicians did not accept Medicaid insurance, 40 4.0899796%) required a referral before the appointment, and 8 0.8179959%) resulted in contacting a personal phone number of the physician.
## The data contains 978 observations, grouped by insurance, of the following 6
## variables:
##
## - Blue Cross/Blue Shield (n = 489):
## - business_days_until_appointment: Mean = 19.59, SD = 22.93, range: [0, 182],
## 29.86% missing
## - Med_sch: 2 levels, namely International Medical Graduate (n = 51), US Senior
## Medical Student (n = 399) and missing (n = 39)
## - gender: 2 levels, namely Male (n = 454) and Female (n = 35)
## - ACOG_District: 9 entries, such as Middle Atlantic (12.68%); West South
## Central (12.47%); New England (12.27%) and 6 others
## - academic_affiliation: 2 levels, namely Not Academic (n = 386), Academic (n =
## 71) and missing (n = 32)
##
## - Medicaid (n = 489):
## - business_days_until_appointment: Mean = 24.55, SD = 23.80, range: [0, 129],
## 54.60% missing
## - Med_sch: 2 levels, namely International Medical Graduate (n = 51), US Senior
## Medical Student (n = 399) and missing (n = 39)
## - gender: 2 levels, namely Male (n = 454) and Female (n = 35)
## - ACOG_District: 9 entries, such as Middle Atlantic (12.68%); West South
## Central (12.47%); New England (12.27%) and 6 others
## - academic_affiliation: 2 levels, namely Not Academic (n = 386), Academic (n =
## 71) and missing (n = 32)Graph each variable
We seem to have many physicians who don’t accept MEDICAID. Should them enter the analysis? It seems that these physicians would not be part of the eligible sample if the goal is to see if individuals who accept both leaves you waiting longer if you are on Medicaid. The research question would be: ** Do providers make patients to wait longer if they are on Medicaid, as compared to Blue Cross? I understand that here the provider must accept both. **
If we include all physicians, then the research question is more like: ** If I have Medicaid, do I wait longer when I try to book an appointment? **
Here we add a variable to the dataset that flag providers who have at least one occurrence of “Physician does not accept MEDICAID”.
df3 datasetThere are 978 providers in the dataset. Many of these providers could not be contacted. In the paper we probably want to report characteristics of contacted providers who have at least one waiting time and are therefore included in the analysis.
Total column is at the provider level. Not meaningful with insurance specific variables, like day of the week.
Insurance Type columns are at the Insurance Type level. Not meaningful with provider level variables.
MEDICAID Acceptance columns are at the provider level. Not meaningful with insurance specific variables.
In the Analysis columns are at the provider level. Not meaningful with insurance specific variables. Here, In the Analysis means those providers that have data for number of days until appointment, and therefore will be used in the analysis. A comparison between “In the Analsyis” and “Not in the Analysis” may give an idea of bias in the data, that is, what kind of providers were easier to reach and because of that are more represented in the sample than they should.
The total number of excluded people is 307.
| Reason for exclusions | n | Percent |
|---|---|---|
| Physician does not accept MEDICAID | 130 | 42.3% |
| Went to voicemail | 61 | 19.9% |
| Number contacted did not correspond to expected office/specialty | 43 | 14% |
| Physician referral required before scheduling appointment | 36 | 11.7% |
| Greater than 5 minutes on hold | 27 | 8.8% |
| Phone not answered or busy signal on repeat calls | 25 | 8.1% |
| Not accepting new patients | 20 | 6.5% |
| Must see midlevel before seeing physician | 10 | 3.3% |
| Closed medical system (e.g. Kaiser or military hospital) | 5 | 1.6% |
| Physician’s personal phone | 4 | 1.3% |
Demographics of all physicians called
| Overall (N=489) | |
|---|---|
| Gender | |
| - Male | 454 (92.8%) |
| - Female | 35 (7.2%) |
| Medical School Training | |
| - US Senior | 399 (88.7%) |
| - International Medical Graduate | 51 (11.3%) |
| Medical School Location | |
| - Allopathic training | 468 (95.7%) |
| - Osteopathic training | 21 (4.3%) |
| Medical School Graduation Year | |
| - 1970 - 1979 | 29 (5.9%) |
| - 1980 - 1989 | 142 (29.0%) |
| - 1990 - 1999 | 196 (40.1%) |
| - 2000 - 2009 | 115 (23.5%) |
| - 2010 to Present | 7 (1.4%) |
| Academic Affiliation | |
| - Not Academic | 416 (85.1%) |
| - Academic | 73 (14.9%) |
| American College of Obstetricians and Gynecologists Districts | |
| - East North Central | 58 (11.9%) |
| - East South Central | 46 (9.4%) |
| - Middle Atlantic | 62 (12.7%) |
| - Mountain | 50 (10.2%) |
| - New England | 60 (12.3%) |
| - Pacific | 56 (11.5%) |
| - South Atlantic | 54 (11.0%) |
| - West North Central | 42 (8.6%) |
| - West South Central | 61 (12.5%) |
| Rurality | |
| - Metropolitan area | 435 (89.1%) |
| - Rural area | 53 (10.9%) |
| Centeral Scheduling | |
| - No | 293 (59.9%) |
| - Yes, central scheduling number | 196 (40.1%) |
| Number of Phone Transfers | |
| - No transfers | 298 (60.9%) |
| - One transfer | 141 (28.8%) |
| - Two transfers | 38 (7.8%) |
| - More than two transfers | 12 (2.5%) |
| Insurance | |
| - Blue Cross/Blue Shield | 245 (50.1%) |
| - Medicaid | 244 (49.9%) |
Here we compare the providers in the analysis because they have some data available and the ones that are excluded from the analysis. Assuming that the “Total” column is representative, we can have an idea if the analyzed providers are skewed.
Caution - Some variables like “Day of the Week” varies with Insurance type and should not be looked at. For this table we selected only the Insurance Type Medicaid.
Waiting time in Days (Log Scale) for Blue Cross/Blue Shield versus Medicaid. The code you provided will create a scatter plot with points representing the relationship between the insurance variable (x-axis) and the days variable (y-axis). Additionally, it includes a line plot that connects points with the same wait time value.
Since the lines connect points with the same npi, if they are horizontal, it suggests similar waiting times for the same npi across both insurance types. If the lines have a slope, it indicates a difference in waiting times.
Here we show a scatterplot that compares the Private and Medicaid times. Notice that the graph is in logarithmic scale. Points above the diagonal line are providers for whom the Medicaid waiting time was longer than the private insurance waiting time.
We also see a strong linear association, indicating that providers with longer waiting time for private insurance tend to also have longer waiting times for Medicaid.
From this information, the scatterplot likely explores the relationship between the time to appointment for patients with Medicaid versus those with Blue Cross/Blue Shield. The use of a logarithmic scale suggests a wide range of days to appointment, and the linear regression line (best fitting line) helps understand if there’s a linear relationship between these two variables on the log scale. The scatterplot might be used to analyze if there’s a systematic difference in waiting times for appointments between the two types of insurance.
If the x-axis represents the days to appointment for Blue Cross/Blue Shield and the y-axis represents the days for Medicaid, a slope less than 45 degrees suggests that for patients with Medicaid, the increase in waiting time is generally less steep compared to those with Blue Cross/Blue Shield for the same increase in waiting time. This could mean that, on average, waiting times increase more slowly for Medicaid patients than for Blue Cross/Blue Shield patients.
The output from the code indicates that the best-fitting line (linear regression line) intersects with the 45-degree line at 33.6 business days until the new patient appointment. This finding means that for your all orthopedists, at approximately 33.6 business days to an appointment, the waiting times for both Blue Cross/Blue Shield and Medicaid are equal. Beyond this point, the relationship between the waiting times for the two types of insurance changes.After this inflection point (>33.6 days) Medicaid patients start to experience longer waiting times compared to Blue Cross/Blue Shield patients for the same period.
## [1] 33.56613## [1] 33.56613## [1] 0The models need to be able to deal with NA in the days
outcome variable (413) and also non-parametric data.
poissonGiven that the “days” variable represents the count of days until a new patient appointment and is a count variable, the Poisson regression model is appropriate for your data. It will model the relationship between the predictor variables and the count of days until a new patient appointment.
$$
\[\begin{align*}\\ P(\text{{Days until New Patient Appointment}} = x) = \frac{{e^{-\lambda} \cdot \lambda^x}}{{x!}}\\ \\where\\ \log(\lambda) = & \beta_0 + \beta_1 \cdot \text{{Patient Insurance}} \\ & + \beta_2 \cdot \text{{Physician Age}}\\ & + \beta_3 \cdot \text{{Physician Academic Affiliation}} \\ & + \beta_4 \cdot \text{{US Census Bureau Divisions}}\\ & + \beta_5 \cdot \text{{Physician Medical Training}} \\ & + \beta_6 \cdot \text{{Physician Gender}} \\ & + \beta_7 \cdot \text{{Central Appointment Phone Number}} \end{align*}\]
$$
Summary of the Poisson model called poisson. According
to the model, Medicaid is associated with 6% longer waiting time in
terms of number of business days. The ICC is quite high (93%) indicating
that waiting times are correlated (similar) within providers.
## GVIF Df GVIF^(1/(2*Df))
## insurance 1.147395 1 1.071165
## gender 1.173654 1 1.083353
## Call_time_minutes 1.146364 1 1.070684
## ACOG_District 1.833029 8 1.038599
## title 1.163310 1 1.078568
## central 1.083690 1 1.041004
## Grd_yr 1.563310 4 1.057440
## Med_sch 1.134416 1 1.065090
## ntransf 1.178679 3 1.027778
## specialty 1.624522 6 1.041263
## cbsatype10 1.130996 1 1.063483
## academic_affiliation 1.146201 1 1.070608## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: poisson ( log )
## Formula: days ~ insurance + gender + Call_time_minutes + ACOG_District +
## title + central + Grd_yr + Med_sch + ntransf + specialty +
## cbsatype10 + academic_affiliation + (1 | name)
## Data: df3
##
## AIC BIC logLik deviance df.resid
## 3273.0 3395.1 -1605.5 3211.0 348
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7739 -0.5589 -0.0460 0.2350 5.2316
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.8864 0.9415
## Number of obs: 379, groups: name, 251
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 1.829938 0.405500 4.513
## insuranceMedicaid 0.084865 0.028297 2.999
## genderFemale -0.009242 0.247119 -0.037
## Call_time_minutes 0.088333 0.017068 5.175
## ACOG_DistrictEast South Central 0.087372 0.253683 0.344
## ACOG_DistrictMiddle Atlantic 0.099713 0.253593 0.393
## ACOG_DistrictMountain 0.324385 0.254090 1.277
## ACOG_DistrictNew England -1.005776 0.742554 -1.354
## ACOG_DistrictPacific 0.657082 0.250808 2.620
## ACOG_DistrictSouth Atlantic 0.167263 0.241946 0.691
## ACOG_DistrictWest North Central 0.484703 0.282408 1.716
## ACOG_DistrictWest South Central -0.001755 0.233461 -0.008
## titleDO -0.131915 0.283699 -0.465
## centralYes 0.098193 0.041602 2.360
## Grd_yr1980 - 1989 0.148390 0.330457 0.449
## Grd_yr1990 - 1999 0.198590 0.327955 0.606
## Grd_yr2000 - 2009 0.278981 0.343898 0.811
## Grd_yr2010 to Present -0.361452 0.685760 -0.527
## Med_schUS Senior Medical Student -0.104906 0.202104 -0.519
## ntransfOne transfer -0.148349 0.041564 -3.569
## ntransfTwo transfers 0.071141 0.078057 0.911
## ntransfMore than two transfers -0.047344 0.164562 -0.288
## specialtyFoot & Ankle Orthopaedic Surgery -0.077104 0.273759 -0.282
## specialtyGeneral Orthopaedic Surgery 0.114956 0.251474 0.457
## specialtyOrthopaedic Surgery of the Spine 0.315155 0.277886 1.134
## specialtyPediatric Orthopaedic Surgery -0.058999 0.349341 -0.169
## specialtySports Medicine Orthopaedics 0.009865 0.270773 0.036
## specialtySurgery of the Hand 0.456144 0.271160 1.682
## cbsatype10Micro -0.152323 0.204791 -0.744
## academic_affiliationAcademic 0.453439 0.222169 2.041
## Pr(>|z|)
## (Intercept) 0.000006398 ***
## insuranceMedicaid 0.002707 **
## genderFemale 0.970167
## Call_time_minutes 0.000000228 ***
## ACOG_DistrictEast South Central 0.730536
## ACOG_DistrictMiddle Atlantic 0.694170
## ACOG_DistrictMountain 0.201725
## ACOG_DistrictNew England 0.175583
## ACOG_DistrictPacific 0.008796 **
## ACOG_DistrictSouth Atlantic 0.489363
## ACOG_DistrictWest North Central 0.086103 .
## ACOG_DistrictWest South Central 0.994002
## titleDO 0.641944
## centralYes 0.018261 *
## Grd_yr1980 - 1989 0.653400
## Grd_yr1990 - 1999 0.544821
## Grd_yr2000 - 2009 0.417232
## Grd_yr2010 to Present 0.598136
## Med_schUS Senior Medical Student 0.603714
## ntransfOne transfer 0.000358 ***
## ntransfTwo transfers 0.362084
## ntransfMore than two transfers 0.773580
## specialtyFoot & Ankle Orthopaedic Surgery 0.778212
## specialtyGeneral Orthopaedic Surgery 0.647579
## specialtyOrthopaedic Surgery of the Spine 0.256745
## specialtyPediatric Orthopaedic Surgery 0.865885
## specialtySports Medicine Orthopaedics 0.970937
## specialtySurgery of the Hand 0.092531 .
## cbsatype10Micro 0.456998
## academic_affiliationAcademic 0.041254 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| days | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 6.23 | 2.82 – 13.80 | <0.001 |
| insurance [Medicaid] | 1.09 | 1.03 – 1.15 | 0.003 |
| genderFemale | 0.99 | 0.61 – 1.61 | 0.970 |
| Call time minutes | 1.09 | 1.06 – 1.13 | <0.001 |
|
ACOG District [East South Central] |
1.09 | 0.66 – 1.79 | 0.731 |
|
ACOG District [Middle Atlantic] |
1.10 | 0.67 – 1.82 | 0.694 |
| ACOG District [Mountain] | 1.38 | 0.84 – 2.28 | 0.202 |
|
ACOG District [New England] |
0.37 | 0.09 – 1.57 | 0.176 |
| ACOG District [Pacific] | 1.93 | 1.18 – 3.15 | 0.009 |
|
ACOG District [South Atlantic] |
1.18 | 0.74 – 1.90 | 0.489 |
|
ACOG District [West North Central] |
1.62 | 0.93 – 2.82 | 0.086 |
|
ACOG District [West South Central] |
1.00 | 0.63 – 1.58 | 0.994 |
| title [DO] | 0.88 | 0.50 – 1.53 | 0.642 |
| central [Yes] | 1.10 | 1.02 – 1.20 | 0.018 |
| Grd yr [1980 - 1989] | 1.16 | 0.61 – 2.22 | 0.653 |
| Grd yr [1990 - 1999] | 1.22 | 0.64 – 2.32 | 0.545 |
| Grd yr [2000 - 2009] | 1.32 | 0.67 – 2.59 | 0.417 |
| Grd yr [2010 to Present] | 0.70 | 0.18 – 2.67 | 0.598 |
|
Med sch [US Senior Medical Student] |
0.90 | 0.61 – 1.34 | 0.604 |
| ntransf [One transfer] | 0.86 | 0.79 – 0.94 | <0.001 |
| ntransf [Two transfers] | 1.07 | 0.92 – 1.25 | 0.362 |
|
ntransf [More than two transfers] |
0.95 | 0.69 – 1.32 | 0.774 |
|
specialty [Foot & Ankle Orthopaedic Surgery] |
0.93 | 0.54 – 1.58 | 0.778 |
|
specialty [General Orthopaedic Surgery] |
1.12 | 0.69 – 1.84 | 0.648 |
|
specialty [Orthopaedic Surgery of the Spine] |
1.37 | 0.79 – 2.36 | 0.257 |
|
specialty [Pediatric Orthopaedic Surgery] |
0.94 | 0.48 – 1.87 | 0.866 |
|
specialty [Sports Medicine Orthopaedics] |
1.01 | 0.59 – 1.72 | 0.971 |
|
specialty [Surgery of the Hand] |
1.58 | 0.93 – 2.68 | 0.093 |
| cbsatype10 [Micro] | 0.86 | 0.57 – 1.28 | 0.457 |
|
academic affiliation [Academic] |
1.57 | 1.02 – 2.43 | 0.041 |
| Random Effects | |||
| σ2 | 0.07 | ||
| τ00 name | 0.89 | ||
| ICC | 0.93 | ||
| N name | 251 | ||
| Observations | 379 | ||
| Marginal R2 / Conditional R2 | 0.133 / 0.936 | ||
## # Indices of model performance
##
## AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma | Score_log | Score_spherical
## --------------------------------------------------------------------------------------------------------------
## 3273.041 | 3278.758 | 3395.104 | 0.936 | 0.133 | 0.926 | 5.868 | 1.000 | -2.998 | 0.042## We fitted a poisson mixed model (estimated using ML and BOBYQA optimizer) to
## predict days with insurance, gender, Call_time_minutes, ACOG_District, title,
## central, Grd_yr, Med_sch, ntransf, specialty, cbsatype10 and
## academic_affiliation (formula: days ~ insurance + gender + Call_time_minutes +
## ACOG_District + title + central + Grd_yr + Med_sch + ntransf + specialty +
## cbsatype10 + academic_affiliation). The model included name as random effect
## (formula: ~1 | name). The model's total explanatory power is substantial
## (conditional R2 = 0.94) and the part related to the fixed effects alone
## (marginal R2) is of 0.13. The model's intercept, corresponding to insurance =
## Blue Cross/Blue Shield, gender = Male, Call_time_minutes = 0, ACOG_District =
## East North Central, title = MD, central = No, Grd_yr = 1970 - 1979, Med_sch =
## International Medical Graduate, ntransf = No transfers, specialty = Adult
## Reconstructive Orthopaedic Surgery, cbsatype10 = Metro and academic_affiliation
## = Not Academic, is at 1.83 (95% CI [1.04, 2.62], p < .001). Within this model:
##
## - The effect of insurance [Medicaid] is statistically significant and positive
## (beta = 0.08, 95% CI [0.03, 0.14], p = 0.003; Std. beta = 0.08, 95% CI [0.03,
## 0.14])
## - The effect of genderFemale is statistically non-significant and negative
## (beta = -9.24e-03, 95% CI [-0.49, 0.48], p = 0.970; Std. beta = -9.24e-03, 95%
## CI [-0.49, 0.48])
## - The effect of Call time minutes is statistically significant and positive
## (beta = 0.09, 95% CI [0.05, 0.12], p < .001; Std. beta = 0.12, 95% CI [0.07,
## 0.17])
## - The effect of ACOG District [East South Central] is statistically
## non-significant and positive (beta = 0.09, 95% CI [-0.41, 0.58], p = 0.731;
## Std. beta = 0.09, 95% CI [-0.41, 0.58])
## - The effect of ACOG District [Middle Atlantic] is statistically
## non-significant and positive (beta = 0.10, 95% CI [-0.40, 0.60], p = 0.694;
## Std. beta = 0.10, 95% CI [-0.40, 0.60])
## - The effect of ACOG District [Mountain] is statistically non-significant and
## positive (beta = 0.32, 95% CI [-0.17, 0.82], p = 0.202; Std. beta = 0.32, 95%
## CI [-0.17, 0.82])
## - The effect of ACOG District [New England] is statistically non-significant
## and negative (beta = -1.01, 95% CI [-2.46, 0.45], p = 0.176; Std. beta = -1.01,
## 95% CI [-2.46, 0.45])
## - The effect of ACOG District [Pacific] is statistically significant and
## positive (beta = 0.66, 95% CI [0.17, 1.15], p = 0.009; Std. beta = 0.66, 95% CI
## [0.17, 1.15])
## - The effect of ACOG District [South Atlantic] is statistically non-significant
## and positive (beta = 0.17, 95% CI [-0.31, 0.64], p = 0.489; Std. beta = 0.17,
## 95% CI [-0.31, 0.64])
## - The effect of ACOG District [West North Central] is statistically
## non-significant and positive (beta = 0.48, 95% CI [-0.07, 1.04], p = 0.086;
## Std. beta = 0.48, 95% CI [-0.07, 1.04])
## - The effect of ACOG District [West South Central] is statistically
## non-significant and negative (beta = -1.75e-03, 95% CI [-0.46, 0.46], p =
## 0.994; Std. beta = -1.75e-03, 95% CI [-0.46, 0.46])
## - The effect of title [DO] is statistically non-significant and negative (beta
## = -0.13, 95% CI [-0.69, 0.42], p = 0.642; Std. beta = -0.13, 95% CI [-0.69,
## 0.42])
## - The effect of central [Yes] is statistically significant and positive (beta =
## 0.10, 95% CI [0.02, 0.18], p = 0.018; Std. beta = 0.10, 95% CI [0.02, 0.18])
## - The effect of Grd yr [1980 - 1989] is statistically non-significant and
## positive (beta = 0.15, 95% CI [-0.50, 0.80], p = 0.653; Std. beta = 0.15, 95%
## CI [-0.50, 0.80])
## - The effect of Grd yr [1990 - 1999] is statistically non-significant and
## positive (beta = 0.20, 95% CI [-0.44, 0.84], p = 0.545; Std. beta = 0.20, 95%
## CI [-0.44, 0.84])
## - The effect of Grd yr [2000 - 2009] is statistically non-significant and
## positive (beta = 0.28, 95% CI [-0.40, 0.95], p = 0.417; Std. beta = 0.28, 95%
## CI [-0.40, 0.95])
## - The effect of Grd yr [2010 to Present] is statistically non-significant and
## negative (beta = -0.36, 95% CI [-1.71, 0.98], p = 0.598; Std. beta = -0.36, 95%
## CI [-1.71, 0.98])
## - The effect of Med sch [US Senior Medical Student] is statistically
## non-significant and negative (beta = -0.10, 95% CI [-0.50, 0.29], p = 0.604;
## Std. beta = -0.10, 95% CI [-0.50, 0.29])
## - The effect of ntransf [One transfer] is statistically significant and
## negative (beta = -0.15, 95% CI [-0.23, -0.07], p < .001; Std. beta = -0.15, 95%
## CI [-0.23, -0.07])
## - The effect of ntransf [Two transfers] is statistically non-significant and
## positive (beta = 0.07, 95% CI [-0.08, 0.22], p = 0.362; Std. beta = 0.07, 95%
## CI [-0.08, 0.22])
## - The effect of ntransf [More than two transfers] is statistically
## non-significant and negative (beta = -0.05, 95% CI [-0.37, 0.28], p = 0.774;
## Std. beta = -0.05, 95% CI [-0.37, 0.28])
## - The effect of specialty [Foot & Ankle Orthopaedic Surgery] is statistically
## non-significant and negative (beta = -0.08, 95% CI [-0.61, 0.46], p = 0.778;
## Std. beta = -0.08, 95% CI [-0.61, 0.46])
## - The effect of specialty [General Orthopaedic Surgery] is statistically
## non-significant and positive (beta = 0.11, 95% CI [-0.38, 0.61], p = 0.648;
## Std. beta = 0.11, 95% CI [-0.38, 0.61])
## - The effect of specialty [Orthopaedic Surgery of the Spine] is statistically
## non-significant and positive (beta = 0.32, 95% CI [-0.23, 0.86], p = 0.257;
## Std. beta = 0.32, 95% CI [-0.23, 0.86])
## - The effect of specialty [Pediatric Orthopaedic Surgery] is statistically
## non-significant and negative (beta = -0.06, 95% CI [-0.74, 0.63], p = 0.866;
## Std. beta = -0.06, 95% CI [-0.74, 0.63])
## - The effect of specialty [Sports Medicine Orthopaedics] is statistically
## non-significant and positive (beta = 9.87e-03, 95% CI [-0.52, 0.54], p = 0.971;
## Std. beta = 9.87e-03, 95% CI [-0.52, 0.54])
## - The effect of specialty [Surgery of the Hand] is statistically
## non-significant and positive (beta = 0.46, 95% CI [-0.08, 0.99], p = 0.093;
## Std. beta = 0.46, 95% CI [-0.08, 0.99])
## - The effect of cbsatype10 [Micro] is statistically non-significant and
## negative (beta = -0.15, 95% CI [-0.55, 0.25], p = 0.457; Std. beta = -0.15, 95%
## CI [-0.55, 0.25])
## - The effect of academic affiliation [Academic] is statistically significant
## and positive (beta = 0.45, 95% CI [0.02, 0.89], p = 0.041; Std. beta = 0.45,
## 95% CI [0.02, 0.89])
##
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald z-distribution approximation.Checking the binned residuals but because the data is non-parametric the residuals will not be normally distributed. Collinearity was tested.
Here we see that the Normal model is quite reasonable for this data, as the residuals looks normally distributed.
## Warning: Probably bad model fit. Only about 74% of the residuals are inside the error bounds.## GVIF Df GVIF^(1/(2*Df))
## insurance 1.147395 1 1.071165
## gender 1.173654 1 1.083353
## Call_time_minutes 1.146364 1 1.070684
## ACOG_District 1.833029 8 1.038599
## title 1.163310 1 1.078568
## central 1.083690 1 1.041004
## Grd_yr 1.563310 4 1.057440
## Med_sch 1.134416 1 1.065090
## ntransf 1.178679 3 1.027778
## specialty 1.624522 6 1.041263
## cbsatype10 1.130996 1 1.063483
## academic_affiliation 1.146201 1 1.070608## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## insurance 1.15 [1.07, 1.33] 1.07 0.87 [0.75, 0.94]
## gender 1.17 [1.09, 1.35] 1.08 0.85 [0.74, 0.92]
## Call_time_minutes 1.15 [1.06, 1.33] 1.07 0.87 [0.75, 0.94]
## ACOG_District 1.83 [1.62, 2.12] 1.35 0.55 [0.47, 0.62]
## title 1.16 [1.08, 1.35] 1.08 0.86 [0.74, 0.93]
## central 1.08 [1.02, 1.30] 1.04 0.92 [0.77, 0.98]
## Grd_yr 1.56 [1.40, 1.80] 1.25 0.64 [0.56, 0.72]
## Med_sch 1.13 [1.06, 1.32] 1.07 0.88 [0.76, 0.95]
## ntransf 1.18 [1.09, 1.36] 1.09 0.85 [0.74, 0.92]
## specialty 1.62 [1.45, 1.87] 1.27 0.62 [0.53, 0.69]
## cbsatype10 1.13 [1.05, 1.32] 1.06 0.88 [0.76, 0.95]
## academic_affiliation 1.15 [1.06, 1.33] 1.07 0.87 [0.75, 0.94]## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.7).
## - For variable: (Whole model)In order to have an idea if there is over-dispersion we divide the Pearson Chi-square by the degree of freedom of the residuals. This ratio should be around 1, with values larger then 1 indicating over-dispersion and lower than 1 indicating under-dispersion. In our case we get the value 1.488 which indicates some over-dispersion. However, if we have overdispersion, our p-value is going to be too small than it should be, so that a significant p-value will be less significant under over-dispersion.
## # Overdispersion test
##
## dispersion ratio = 1.600
## Pearson's Chi-Squared = 556.824
## p-value = < 0.001## chisq ratio rdf
## 556.823967290310292810 1.600068871523880221 348.000000000000000000
## p
## 0.000000000007349562## Warning: Autocorrelated residuals detected (p < .001).## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## insurance 1.15 [1.07, 1.33] 1.07 0.87 [0.75, 0.94]
## gender 1.17 [1.09, 1.35] 1.08 0.85 [0.74, 0.92]
## Call_time_minutes 1.15 [1.06, 1.33] 1.07 0.87 [0.75, 0.94]
## ACOG_District 1.83 [1.62, 2.12] 1.35 0.55 [0.47, 0.62]
## title 1.16 [1.08, 1.35] 1.08 0.86 [0.74, 0.93]
## central 1.08 [1.02, 1.30] 1.04 0.92 [0.77, 0.98]
## Grd_yr 1.56 [1.40, 1.80] 1.25 0.64 [0.56, 0.72]
## Med_sch 1.13 [1.06, 1.32] 1.07 0.88 [0.76, 0.95]
## ntransf 1.18 [1.09, 1.36] 1.09 0.85 [0.74, 0.92]
## specialty 1.62 [1.45, 1.87] 1.27 0.62 [0.53, 0.69]
## cbsatype10 1.13 [1.05, 1.32] 1.06 0.88 [0.76, 0.95]
## academic_affiliation 1.15 [1.06, 1.33] 1.07 0.87 [0.75, 0.94]## [1] FALSETesting assumptions you can use the logLik function to get the log-likelihood of the model, and calculate the residual deviance as -2 * logLik(model). The residual degrees of freedom can be computed as the number of observations minus the number of parameters estimated (which includes both fixed effects and random effects).
The number of parameters estimated can be calculated as the number of fixed effects plus the number of random effects parameters. The number of fixed effects can be obtained from the length of fixef(model), and the number of random effects parameters can be obtained from the length of VarCorr(model).
If the dispersion parameter is considerably greater than 1, it indicates overdispersion. If it is less than 1, it indicates underdispersion. A value around 1 is considered ideal for Poisson regression.
## 'log Lik.' 3.405133 (df=31)This command will create a residuals plot that can help you check the assumptions of your Poisson regression model. If the plot shows a random scatter, then the assumptions are likely met. If the plot shows a clear pattern or trend, then the assumptions might not be met, and you might need to consider a different modeling approach.
The Poisson regression assumes that the log of the expected count is a linear function of the predictors. One way to check this is to plot the observed counts versus the predicted counts and see if the relationship looks linear.
We will need to check interaction of insurance with the
other significant variables. “significant variables in the model
estimates” refer to predictors that have a significant effect on the
response variable individually, while the “ANOVA” assesses the overall
significance of the model and the joint significance of all
predictors.
| x |
|---|
| insurance |
| Call_time_minutes |
| central |
| ntransf |
| academic_affiliation |
To include interaction terms in a regression model, you can use the : operator or the * operator in the formula. The : operator represents the interaction between two variables, while the * operator represents the interaction and also includes the main effects of the two variables. This will include interactions between insurance and each of the other significant variables (academic_affiliation, ACOG_District, central), in addition to the main effects of these variables.
Please note that interpreting interaction effects can be complex, especially in nonlinear models such as Poisson regression. The coefficients for the interaction terms represent the difference in the log rate of days for a one-unit change in insurance, for different levels of the other variables. However, the actual effects on the rate of days can vary depending on the values of the other variables.
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: poisson ( log )
## Formula: days ~ insurance + gender + Call_time_minutes + ACOG_District +
## title + central + Grd_yr + Med_sch + ntransf + specialty +
## cbsatype10 + academic_affiliation + (1 | name)
## Data: df3
##
## AIC BIC logLik deviance df.resid
## 3273.0 3395.1 -1605.5 3211.0 348
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7739 -0.5589 -0.0460 0.2350 5.2316
##
## Random effects:
## Groups Name Variance Std.Dev.
## name (Intercept) 0.8864 0.9415
## Number of obs: 379, groups: name, 251
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 1.829938 0.405500 4.513
## insuranceMedicaid 0.084865 0.028297 2.999
## genderFemale -0.009242 0.247119 -0.037
## Call_time_minutes 0.088333 0.017068 5.175
## ACOG_DistrictEast South Central 0.087372 0.253683 0.344
## ACOG_DistrictMiddle Atlantic 0.099713 0.253593 0.393
## ACOG_DistrictMountain 0.324385 0.254090 1.277
## ACOG_DistrictNew England -1.005776 0.742554 -1.354
## ACOG_DistrictPacific 0.657082 0.250808 2.620
## ACOG_DistrictSouth Atlantic 0.167263 0.241946 0.691
## ACOG_DistrictWest North Central 0.484703 0.282408 1.716
## ACOG_DistrictWest South Central -0.001755 0.233461 -0.008
## titleDO -0.131915 0.283699 -0.465
## centralYes 0.098193 0.041602 2.360
## Grd_yr1980 - 1989 0.148390 0.330457 0.449
## Grd_yr1990 - 1999 0.198590 0.327955 0.606
## Grd_yr2000 - 2009 0.278981 0.343898 0.811
## Grd_yr2010 to Present -0.361452 0.685760 -0.527
## Med_schUS Senior Medical Student -0.104906 0.202104 -0.519
## ntransfOne transfer -0.148349 0.041564 -3.569
## ntransfTwo transfers 0.071141 0.078057 0.911
## ntransfMore than two transfers -0.047344 0.164562 -0.288
## specialtyFoot & Ankle Orthopaedic Surgery -0.077104 0.273759 -0.282
## specialtyGeneral Orthopaedic Surgery 0.114956 0.251474 0.457
## specialtyOrthopaedic Surgery of the Spine 0.315155 0.277886 1.134
## specialtyPediatric Orthopaedic Surgery -0.058999 0.349341 -0.169
## specialtySports Medicine Orthopaedics 0.009865 0.270773 0.036
## specialtySurgery of the Hand 0.456144 0.271160 1.682
## cbsatype10Micro -0.152323 0.204791 -0.744
## academic_affiliationAcademic 0.453439 0.222169 2.041
## Pr(>|z|)
## (Intercept) 0.000006398 ***
## insuranceMedicaid 0.002707 **
## genderFemale 0.970167
## Call_time_minutes 0.000000228 ***
## ACOG_DistrictEast South Central 0.730536
## ACOG_DistrictMiddle Atlantic 0.694170
## ACOG_DistrictMountain 0.201725
## ACOG_DistrictNew England 0.175583
## ACOG_DistrictPacific 0.008796 **
## ACOG_DistrictSouth Atlantic 0.489363
## ACOG_DistrictWest North Central 0.086103 .
## ACOG_DistrictWest South Central 0.994002
## titleDO 0.641944
## centralYes 0.018261 *
## Grd_yr1980 - 1989 0.653400
## Grd_yr1990 - 1999 0.544821
## Grd_yr2000 - 2009 0.417232
## Grd_yr2010 to Present 0.598136
## Med_schUS Senior Medical Student 0.603714
## ntransfOne transfer 0.000358 ***
## ntransfTwo transfers 0.362084
## ntransfMore than two transfers 0.773580
## specialtyFoot & Ankle Orthopaedic Surgery 0.778212
## specialtyGeneral Orthopaedic Surgery 0.647579
## specialtyOrthopaedic Surgery of the Spine 0.256745
## specialtyPediatric Orthopaedic Surgery 0.865885
## specialtySports Medicine Orthopaedics 0.970937
## specialtySurgery of the Hand 0.092531 .
## cbsatype10Micro 0.456998
## academic_affiliationAcademic 0.041254 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## [1] "insurance" "Call_time_minutes" "central"
## [4] "ntransf" "academic_affiliation"## $data
## ntransf = No transfers:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 16.36658 1.341873 Inf 13.93700 19.21971
## Medicaid 18.11542 1.502420 Inf 15.39760 21.31295
##
## ntransf = One transfer:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 14.72457 1.272135 Inf 12.43090 17.44145
## Medicaid 16.29795 1.421332 Inf 13.73726 19.33595
##
## ntransf = Two transfers:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 16.32577 1.668230 Inf 13.36271 19.94586
## Medicaid 18.07024 1.803732 Inf 14.85932 21.97500
##
## ntransf = More than two transfers:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 17.26735 2.414334 Inf 13.12835 22.71126
## Medicaid 19.11243 2.719718 Inf 14.46070 25.26055
##
## Results are averaged over the levels of: central, academic_affiliation
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $plot## $emmeans
## ntransf = No transfers:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 16.4 1.34 Inf 13.9 19.2
## Medicaid 18.1 1.50 Inf 15.4 21.3
##
## ntransf = One transfer:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 14.7 1.27 Inf 12.4 17.4
## Medicaid 16.3 1.42 Inf 13.7 19.3
##
## ntransf = Two transfers:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 16.3 1.67 Inf 13.4 19.9
## Medicaid 18.1 1.80 Inf 14.9 22.0
##
## ntransf = More than two transfers:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 17.3 2.41 Inf 13.1 22.7
## Medicaid 19.1 2.72 Inf 14.5 25.3
##
## Results are averaged over the levels of: central, academic_affiliation
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## ntransf = No transfers:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## ntransf = One transfer:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## ntransf = Two transfers:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## ntransf = More than two transfers:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## Results are averaged over the levels of: central, academic_affiliation
## Tests are performed on the log scale## $data
## Call_time_minutes = 2.809223:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 16.14451 1.394132 Inf 13.63079 19.12181
## Medicaid 17.86962 1.556476 Inf 15.06517 21.19612
##
## Results are averaged over the levels of: ntransf, central, academic_affiliation
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $plot## $emmeans
## Call_time_minutes = 2.81:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 16.1 1.39 Inf 13.6 19.1
## Medicaid 17.9 1.56 Inf 15.1 21.2
##
## Results are averaged over the levels of: ntransf, central, academic_affiliation
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## Call_time_minutes = 2.8092225241479:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## Results are averaged over the levels of: ntransf, central, academic_affiliation
## Tests are performed on the log scale## [1] "insurance" "Call_time_minutes" "central"
## [4] "ntransf" "academic_affiliation"## $data
## central = No:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 15.09737 1.329513 Inf 12.70405 17.94156
## Medicaid 16.71058 1.479933 Inf 14.04776 19.87816
##
## central = Yes:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 17.26429 1.521281 Inf 14.52590 20.51890
## Medicaid 19.10904 1.702312 Inf 16.04761 22.75451
##
## Results are averaged over the levels of: ntransf, academic_affiliation
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $plot## $emmeans
## central = No:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 15.1 1.33 Inf 12.7 17.9
## Medicaid 16.7 1.48 Inf 14.0 19.9
##
## central = Yes:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 17.3 1.52 Inf 14.5 20.5
## Medicaid 19.1 1.70 Inf 16.0 22.8
##
## Results are averaged over the levels of: ntransf, academic_affiliation
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## central = No:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## central = Yes:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## Results are averaged over the levels of: ntransf, academic_affiliation
## Tests are performed on the log scale## [1] "insurance" "Call_time_minutes" "central"
## [4] "ntransf" "academic_affiliation"## $data
## academic_affiliation = Not Academic:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 12.29590 0.835998 Inf 10.76186 14.04862
## Medicaid 13.60977 0.943823 Inf 11.88012 15.59124
##
## academic_affiliation = Academic:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 21.19773 3.241173 Inf 15.70866 28.60485
## Medicaid 23.46279 3.593102 Inf 17.37906 31.67620
##
## Results are averaged over the levels of: ntransf, central
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $plot## $emmeans
## academic_affiliation = Not Academic:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 12.3 0.836 Inf 10.8 14.0
## Medicaid 13.6 0.944 Inf 11.9 15.6
##
## academic_affiliation = Academic:
## insurance rate SE df asymp.LCL asymp.UCL
## Blue Cross/\nBlue Shield 21.2 3.241 Inf 15.7 28.6
## Medicaid 23.5 3.593 Inf 17.4 31.7
##
## Results are averaged over the levels of: ntransf, central
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## academic_affiliation = Not Academic:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## academic_affiliation = Academic:
## contrast ratio SE df null z.ratio p.value
## (Blue Cross/\nBlue Shield) / Medicaid 0.903 0.021 Inf 1 -4.374 <.0001
##
## Results are averaged over the levels of: ntransf, central
## Tests are performed on the log scale## - Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K, Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2023). _rmarkdown: Dynamic Documents for R_. R package version 2.25, <https://github.com/rstudio/rmarkdown>. Xie Y, Allaire J, Grolemund G (2018). _R Markdown: The Definitive Guide_. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, <https://bookdown.org/yihui/rmarkdown>. Xie Y, Dervieux C, Riederer E (2020). _R Markdown Cookbook_. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 9780367563837, <https://bookdown.org/yihui/rmarkdown-cookbook>.
## - Arel-Bundock V (2022). "modelsummary: Data and Model Summaries in R." _Journal of Statistical Software_, *103*(1), 1-23. doi:10.18637/jss.v103.i01 <https://doi.org/10.18637/jss.v103.i01>.
## - Attali D, Baker C (2023). _ggExtra: Add Marginal Histograms to 'ggplot2', and More 'ggplot2' Enhancements_. R package version 0.10.1, <https://CRAN.R-project.org/package=ggExtra>.
## - Auguie B (2017). _gridExtra: Miscellaneous Functions for "Grid" Graphics_. R package version 2.3, <https://CRAN.R-project.org/package=gridExtra>.
## - Barrett T, Dowle M, Srinivasan A (2023). _data.table: Extension of `data.frame`_. R package version 1.14.10, <https://CRAN.R-project.org/package=data.table>.
## - Bates D, Mächler M, Bolker B, Walker S (2015). "Fitting Linear Mixed-Effects Models Using lme4." _Journal of Statistical Software_, *67*(1), 1-48. doi:10.18637/jss.v067.i01 <https://doi.org/10.18637/jss.v067.i01>.
## - Bates D, Maechler M, Jagan M (2023). _Matrix: Sparse and Dense Matrix Classes and Methods_. R package version 1.6-3, <https://CRAN.R-project.org/package=Matrix>.
## - Ben-Shachar MS, Lüdecke D, Makowski D (2020). "effectsize: Estimation of Effect Size Indices and Standardized Parameters." _Journal of Open Source Software_, *5*(56), 2815. doi:10.21105/joss.02815 <https://doi.org/10.21105/joss.02815>, <https://doi.org/10.21105/joss.02815>.
## - Bengtsson H (2021). "A Unifying Framework for Parallel and Distributed Processing in R using Futures." _The R Journal_, *13*(2), 208-227. doi:10.32614/RJ-2021-048 <https://doi.org/10.32614/RJ-2021-048>, <https://doi.org/10.32614/RJ-2021-048>.
## - Bolker B, Robinson D (2022). _broom.mixed: Tidying Methods for Mixed Models_. R package version 0.2.9.4, <https://CRAN.R-project.org/package=broom.mixed>.
## - Couch SP, Bray AP, Ismay C, Chasnovski E, Baumer BS, Çetinkaya-Rundel M (2021). "infer: An R package for tidyverse-friendly statistical inference." _Journal of Open Source Software_, *6*(65), 3661. doi:10.21105/joss.03661 <https://doi.org/10.21105/joss.03661>.
## - Croissant Y, Millo G (2018). _Panel Data Econometrics with R_. Wiley. Croissant Y, Millo G (2008). "Panel Data Econometrics in R: The plm Package." _Journal of Statistical Software_, *27*(2), 1-43. doi:10.18637/jss.v027.i02 <https://doi.org/10.18637/jss.v027.i02>. Millo G (2017). "Robust Standard Error Estimators for Panel Models: A Unifying Approach." _Journal of Statistical Software_, *82*(3), 1-27. doi:10.18637/jss.v082.i03 <https://doi.org/10.18637/jss.v082.i03>.
## - Daróczi G, Tsegelskyi R (2022). _pander: An R 'Pandoc' Writer_. R package version 0.6.5, <https://CRAN.R-project.org/package=pander>.
## - Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). _Rcpp: Seamless R and C++ Integration_. R package version 1.0.12, <https://CRAN.R-project.org/package=Rcpp>. Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). "Extending R with C++: A Brief Introduction to Rcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
## - Falissard B (2022). _psy: Various Procedures Used in Psychometrics_. R package version 1.2, <https://CRAN.R-project.org/package=psy>.
## - Fox J, Venables B, Damico A, Salverda AP (2021). _english: Translate Integers into English_. R package version 1.2-6, <https://CRAN.R-project.org/package=english>.
## - Fox J, Weisberg S (2019). _An R Companion to Applied Regression_, 3rd edition. Sage, Thousand Oaks CA. <https://socialsciences.mcmaster.ca/jfox/Books/Companion/index.html>. Fox J, Weisberg S (2018). "Visualizing Fit and Lack of Fit in Complex Regression Models with Predictor Effect Plots and Partial Residuals." _Journal of Statistical Software_, *87*(9), 1-27. doi:10.18637/jss.v087.i09 <https://doi.org/10.18637/jss.v087.i09>. Fox J (2003). "Effect Displays in R for Generalised Linear Models." _Journal of Statistical Software_, *8*(15), 1-27. doi:10.18637/jss.v008.i15 <https://doi.org/10.18637/jss.v008.i15>. Fox J, Hong J (2009). "Effect Displays in R for Multinomial and Proportional-Odds Logit Models: Extensions to the effects Package." _Journal of Statistical Software_, *32*(1), 1-24. doi:10.18637/jss.v032.i01 <https://doi.org/10.18637/jss.v032.i01>.
## - Fox J, Weisberg S (2019). _An R Companion to Applied Regression_, Third edition. Sage, Thousand Oaks CA. <https://socialsciences.mcmaster.ca/jfox/Books/Companion/>.
## - Fox J, Weisberg S, Price B (2022). _carData: Companion to Applied Regression Data Sets_. R package version 3.0-5, <https://CRAN.R-project.org/package=carData>.
## - Frick H, Chow F, Kuhn M, Mahoney M, Silge J, Wickham H (2023). _rsample: General Resampling Infrastructure_. R package version 1.2.0, <https://CRAN.R-project.org/package=rsample>.
## - Gohel D (2023). _officer: Manipulation of Microsoft Word and PowerPoint Documents_. R package version 0.6.3, <https://CRAN.R-project.org/package=officer>.
## - Gohel D, Skintzos P (2023). _flextable: Functions for Tabular Reporting_. R package version 0.9.4, <https://CRAN.R-project.org/package=flextable>.
## - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
## - Grothendieck G, Kates L, Petzoldt T (2016). _proto: Prototype Object-Based Programming_. R package version 1.0.0, <https://CRAN.R-project.org/package=proto>.
## - Halekoh U, Højsgaard S, Yan J (2006). "The R Package geepack for Generalized Estimating Equations." _Journal of Statistical Software_, *15/2*, 1-11. Yan J, Fine JP (2004). "Estimating Equations for Association Structures." _Statistics in Medicine_, *23*, 859-880. Yan J (2002). "geepack: Yet Another Package for Generalized Estimating Equations." _R-News_, *2/3*, 12-14.
## - Hayashi H, Kojima H, Nishida K, Saito K, Yasuda Y (2023). _exploratory: R package for Exploratory_. R package version 9.0.3, <https://github.com/exploratory-io/exploratory_func>.
## - Heinzen E, Sinnwell J, Atkinson E, Gunderson T, Dougherty G (2021). _arsenal: An Arsenal of 'R' Functions for Large-Scale Statistical Summaries_. R package version 3.6.3, <https://CRAN.R-project.org/package=arsenal>.
## - Kassambara A (2023). _ggpubr: 'ggplot2' Based Publication Ready Plots_. R package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
## - Kuhn M (2023). _modeldata: Data Sets Useful for Modeling Examples_. R package version 1.2.0, <https://CRAN.R-project.org/package=modeldata>.
## - Kuhn M (2023). _tune: Tidy Tuning Tools_. R package version 1.1.2, <https://CRAN.R-project.org/package=tune>.
## - Kuhn M, Couch S (2023). _workflowsets: Create a Collection of 'tidymodels' Workflows_. R package version 1.0.1, <https://CRAN.R-project.org/package=workflowsets>.
## - Kuhn M, Frick H (2023). _dials: Tools for Creating Tuning Parameter Values_. R package version 1.2.0, <https://CRAN.R-project.org/package=dials>.
## - Kuhn M, Vaughan D (2023). _parsnip: A Common API to Modeling and Analysis Functions_. R package version 1.1.1, <https://CRAN.R-project.org/package=parsnip>.
## - Kuhn M, Vaughan D, Hvitfeldt E (2023). _yardstick: Tidy Characterizations of Model Performance_. R package version 1.2.0, <https://CRAN.R-project.org/package=yardstick>.
## - Kuhn M, Wickham H (2020). _Tidymodels: a collection of packages for modeling and machine learning using tidyverse principles._. <https://www.tidymodels.org>.
## - Kuhn M, Wickham H, Hvitfeldt E (2023). _recipes: Preprocessing and Feature Engineering Steps for Modeling_. R package version 1.0.9, <https://CRAN.R-project.org/package=recipes>.
## - Kuhn, Max (2008). "Building Predictive Models in R Using the caret Package." _Journal of Statistical Software_, *28*(5), 1–26. doi:10.18637/jss.v028.i05 <https://doi.org/10.18637/jss.v028.i05>, <https://www.jstatsoft.org/index.php/jss/article/view/v028i05>.
## - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package: Tests in Linear Mixed Effects Models." _Journal of Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13 <https://doi.org/10.18637/jss.v082.i13>.
## - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version 1.8.9, <https://CRAN.R-project.org/package=emmeans>.
## - Liaw A, Wiener M (2002). "Classification and Regression by randomForest." _R News_, *2*(3), 18-22. <https://CRAN.R-project.org/doc/Rnews/>.
## - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
## - Lüdecke D (2023). _sjPlot: Data Visualization for Statistics in Social Science_. R package version 2.8.15, <https://CRAN.R-project.org/package=sjPlot>.
## - Lüdecke D, Ben-Shachar M, Patil I, Makowski D (2020). "Extracting, Computing and Exploring the Parameters of Statistical Models using R." _Journal of Open Source Software_, *5*(53), 2445. doi:10.21105/joss.02445 <https://doi.org/10.21105/joss.02445>.
## - Lüdecke D, Ben-Shachar M, Patil I, Waggoner P, Makowski D (2021). "performance: An R Package for Assessment, Comparison and Testing of Statistical Models." _Journal of Open Source Software_, *6*(60), 3139. doi:10.21105/joss.03139 <https://doi.org/10.21105/joss.03139>.
## - Lüdecke D, Ben-Shachar M, Patil I, Wiernik B, Bacher E, Thériault R, Makowski D (2022). "easystats: Framework for Easy Statistical Modeling, Visualization, and Reporting." _CRAN_. R package, <https://easystats.github.io/easystats/>.
## - Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021). "see: An R Package for Visualizing Statistical Models." _Journal of Open Source Software_, *6*(64), 3393. doi:10.21105/joss.03393 <https://doi.org/10.21105/joss.03393>.
## - Lüdecke D, Waggoner P, Makowski D (2019). "insight: A Unified Interface to Access Information from Model Objects in R." _Journal of Open Source Software_, *4*(38), 1412. doi:10.21105/joss.01412 <https://doi.org/10.21105/joss.01412>.
## - Makowski D, Ben-Shachar M, Lüdecke D (2019). "bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework." _Journal of Open Source Software_, *4*(40), 1541. doi:10.21105/joss.01541 <https://doi.org/10.21105/joss.01541>, <https://joss.theoj.org/papers/10.21105/joss.01541>.
## - Makowski D, Ben-Shachar M, Patil I, Lüdecke D (2020). "Estimation of Model-Based Predictions, Contrasts and Means." _CRAN_. <https://github.com/easystats/modelbased>.
## - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023). "Automated Results Reporting as a Practical Tool to Improve Reproducibility and Methodological Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
## - Makowski D, Wiernik B, Patil I, Lüdecke D, Ben-Shachar M (2022). "correlation: Methods for Correlation Analysis." Version 0.8.3, <https://CRAN.R-project.org/package=correlation>. Makowski D, Ben-Shachar M, Patil I, Lüdecke D (2020). "Methods and Algorithms for Correlation Analysis in R." _Journal of Open Source Software_, *5*(51), 2306. doi:10.21105/joss.02306 <https://doi.org/10.21105/joss.02306>, <https://joss.theoj.org/papers/10.21105/joss.02306>.
## - Merkle E, You D (2023). _nonnest2: Tests of Non-Nested Models_. R package version 0.5-6, <https://CRAN.R-project.org/package=nonnest2>.
## - Meyer D, Dimitriadou E, Hornik K, Weingessel A, Leisch F (2023). _e1071: Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien_. R package version 1.7-14, <https://CRAN.R-project.org/package=e1071>.
## - Microsoft, Weston S (2022). _foreach: Provides Foreach Looping Construct_. R package version 1.5.2, <https://CRAN.R-project.org/package=foreach>.
## - Müller K (2020). _here: A Simpler Way to Find Your Files_. R package version 1.0.1, <https://CRAN.R-project.org/package=here>.
## - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
## - Ooms J (2023). _writexl: Export Data Frames to Excel 'xlsx' Format_. R package version 1.4.2, <https://CRAN.R-project.org/package=writexl>.
## - Paluszynska A, Biecek P, Jiang Y (2020). _randomForestExplainer: Explaining and Visualizing Random Forests in Terms of Variable Importance_. R package version 0.10.1, <https://CRAN.R-project.org/package=randomForestExplainer>.
## - Patil I, Makowski D, Ben-Shachar M, Wiernik B, Bacher E, Lüdecke D (2022). "datawizard: An R Package for Easy Data Preparation and Statistical Transformations." _Journal of Open Source Software_, *7*(78), 4684. doi:10.21105/joss.04684 <https://doi.org/10.21105/joss.04684>.
## - R Core Team (2023). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
## - Rich B (2023). _table1: Tables of Descriptive Statistics in HTML_. R package version 1.4.3, <https://CRAN.R-project.org/package=table1>.
## - Robinson D, Hayes A, Couch S (2023). _broom: Convert Statistical Objects into Tidy Tibbles_. R package version 1.0.5, <https://CRAN.R-project.org/package=broom>.
## - Rosseel Y (2012). "lavaan: An R Package for Structural Equation Modeling." _Journal of Statistical Software_, *48*(2), 1-36. doi:10.18637/jss.v048.i02 <https://doi.org/10.18637/jss.v048.i02>.
## - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer, New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
## - Schloerke B, Cook D, Larmarange J, Briatte F, Marbach M, Thoen E, Elberg A, Crowley J (2023). _GGally: Extension to 'ggplot2'_. R package version 2.2.0, <https://CRAN.R-project.org/package=GGally>.
## - Tang Y, Horikoshi M, Li W (2016). "ggfortify: Unified Interface to Visualize Statistical Result of Popular R Packages." _The R Journal_, *8*(2), 474-485. doi:10.32614/RJ-2016-060 <https://doi.org/10.32614/RJ-2016-060>, <https://doi.org/10.32614/RJ-2016-060>. Horikoshi M, Tang Y (2018). _ggfortify: Data Visualization Tools for Statistical Analysis Results_. <https://CRAN.R-project.org/package=ggfortify>.
## - Ushey K, Wickham H (2023). _renv: Project Environments_. R package version 1.0.3, <https://CRAN.R-project.org/package=renv>.
## - Vaughan D, Couch S (2023). _workflows: Modeling Workflows_. R package version 1.1.3, <https://CRAN.R-project.org/package=workflows>.
## - Vaughan D, Dancho M (2022). _furrr: Apply Mapping Functions in Parallel using Futures_. R package version 0.3.1, <https://CRAN.R-project.org/package=furrr>.
## - Wang T, Merkle EC (2018). "merDeriv: Derivative Computations for Linear Mixed Effects Models with Application to Robust Standard Errors." _Journal of Statistical Software, Code Snippets_, *87*(1), 1-16. doi:10.18637/jss.v087.c01 <https://doi.org/10.18637/jss.v087.c01>.
## - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
## - Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
## - Wickham H (2023). _stringr: Simple, Consistent Wrappers for Common String Operations_. R package version 1.5.1, <https://CRAN.R-project.org/package=stringr>.
## - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. R package version 1.1.4, <https://CRAN.R-project.org/package=dplyr>.
## - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package version 1.0.2, <https://CRAN.R-project.org/package=purrr>.
## - Wickham H, Hester J, Bryan J (2023). _readr: Read Rectangular Text Data_. R package version 2.1.4, <https://CRAN.R-project.org/package=readr>.
## - Wickham H, Hester J, Csárdi G (2023). _pkgbuild: Find Tools Needed to Build R Packages_. R package version 1.4.3, <https://CRAN.R-project.org/package=pkgbuild>.
## - Wickham H, Pedersen T, Seidel D (2023). _scales: Scale Functions for Visualization_. R package version 1.3.0, <https://CRAN.R-project.org/package=scales>.
## - Wickham H, Vaughan D, Girlich M (2023). _tidyr: Tidy Messy Data_. R package version 1.3.0, <https://CRAN.R-project.org/package=tidyr>.
## - Wilke C (2023). _cowplot: Streamlined Plot Theme and Plot Annotations for 'ggplot2'_. R package version 1.1.2, <https://CRAN.R-project.org/package=cowplot>.
## - Xie Y (2023). _knitr: A General-Purpose Package for Dynamic Report Generation in R_. R package version 1.45, <https://yihui.org/knitr/>. Xie Y (2015). _Dynamic Documents with R and knitr_, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, <https://yihui.org/knitr/>. Xie Y (2014). "knitr: A Comprehensive Tool for Reproducible Research in R." In Stodden V, Leisch F, Peng RD (eds.), _Implementing Reproducible Computational Research_. Chapman and Hall/CRC. ISBN 978-1466561595.
## - Zeileis A, Grothendieck G (2005). "zoo: S3 Infrastructure for Regular and Irregular Time Series." _Journal of Statistical Software_, *14*(6), 1-27. doi:10.18637/jss.v014.i06 <https://doi.org/10.18637/jss.v014.i06>.
## - Zeileis A, Hothorn T (2002). "Diagnostic Checking in Regression Relationships." _R News_, *2*(3), 7-10. <https://CRAN.R-project.org/doc/Rnews/>.
## - Zeileis A, Köll S, Graham N (2020). "Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in R." _Journal of Statistical Software_, *95*(1), 1-36. doi:10.18637/jss.v095.i01 <https://doi.org/10.18637/jss.v095.i01>. Zeileis A (2004). "Econometric Computing with HC and HAC Covariance Matrix Estimators." _Journal of Statistical Software_, *11*(10), 1-17. doi:10.18637/jss.v011.i10 <https://doi.org/10.18637/jss.v011.i10>. Zeileis A (2006). "Object-Oriented Computation of Sandwich Estimators." _Journal of Statistical Software_, *16*(9), 1-16. doi:10.18637/jss.v016.i09 <https://doi.org/10.18637/jss.v016.i09>.
## - Zhu H (2021). _kableExtra: Construct Complex Table with 'kable' and Pipe Syntax_. R package version 1.3.4, <https://CRAN.R-project.org/package=kableExtra>.–>
–>
–>
–>
–>
–>
–>
–>