Context

This file holds the output for the logistic regression model of psychology today mental health providers insurance acceptance:

Code book for variables in full model:

c_provider_1k: the supply of all types of mental health providers per 1000k county residents, centered around its mean
- Mean: 0.94
provs: provider type in categories.
- Levels: “Masters level therapists and counselors”, “All other therapists and counselors”, “Psychologists/PhD”, “Advanced practice providers”
- Level names were changed to be better on output
any_child: does a provider also see clients under 18?
- Binary
Rurality: Binary–urban or rural.
- Based on Categories of the 2013 and 2006 NCHS Urban-Rural Classification Scheme for Counties.
- Metro distinctions 1-4 = urban
- metro distinctions 5-6 = rural (micropoliton and non-core)
income_med: Binary
- whether provider is in a county with a median household income that is above or below the median of the sample.
- Data from ACS 5-year survey
group_pr: is the provider part of a group practice?
- Determined by grouping based on exact addresses.
- Binary
frac_white_med: binary
- Data from 5-year ACS
- Binary, whether provider is in a county with a fraction of county population that is white that is above or below the median of the sample

Load Required Packages

Data Processing

Univariate model:

I chose to change the reference category to masters level providers since it is the largest group

data_acs_inc$provs <- contr_code_treatment(data_acs_inc$provs, 
                                                   base = "masters")

uni_prov <- glm(insurance_flag ~ provs, data = data_acs_inc, 
                family = "binomial")
summary(uni_prov)

## 
## Call:
## glm(formula = insurance_flag ~ provs, family = "binomial", data = data_acs_inc)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               0.50325    0.00727  69.222   <2e-16 ***
## provs.advanced-masters   -0.38909    0.04038  -9.635   <2e-16 ***
## provs.counselors-masters  0.44135    0.01134  38.917   <2e-16 ***
## provs.psych-masters      -0.24335    0.01437 -16.938   <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: 226890  on 175082  degrees of freedom
## Residual deviance: 224196  on 175079  degrees of freedom
## AIC: 224204
## 
## Number of Fisher Scoring iterations: 4

round(
  exp(
    cbind(OR = coef(uni_prov), confint(uni_prov))
    ), 
  2)

## Waiting for profiling to be done...

##                            OR 2.5 % 97.5 %
## (Intercept)              1.65  1.63   1.68
## provs.advanced-masters   0.68  0.63   0.73
## provs.counselors-masters 1.55  1.52   1.59
## provs.psych-masters      0.78  0.76   0.81

Multivariate model, no interaction

model1 <- glm(insurance_flag ~ provs + c_provider_1k + rurality + any_child 
              + income_med + group_pr + frac_white_med + telehealth, data = data_acs_inc, 
              family = "binomial")
summary(model1)

## 
## Call:
## glm(formula = insurance_flag ~ provs + c_provider_1k + rurality + 
##     any_child + income_med + group_pr + frac_white_med + telehealth, 
##     family = "binomial", data = data_acs_inc)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               0.233783   0.013126  17.811  < 2e-16 ***
## provs.advanced-masters   -0.402757   0.040745  -9.885  < 2e-16 ***
## provs.counselors-masters  0.471490   0.011467  41.116  < 2e-16 ***
## provs.psych-masters      -0.194728   0.014580 -13.356  < 2e-16 ***
## c_provider_1k            -0.134196   0.008298 -16.173  < 2e-16 ***
## ruralityRural             0.167428   0.028521   5.870 4.35e-09 ***
## any_child                 0.127924   0.010715  11.939  < 2e-16 ***
## income_med               -0.226293   0.011024 -20.528  < 2e-16 ***
## group_pr                  0.228906   0.010577  21.641  < 2e-16 ***
## frac_white_med            0.305216   0.010405  29.334  < 2e-16 ***
## telehealth                0.134975   0.011875  11.366  < 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: 226890  on 175082  degrees of freedom
## Residual deviance: 221376  on 175072  degrees of freedom
## AIC: 221398
## 
## Number of Fisher Scoring iterations: 4

round(
  exp(
    cbind(OR = coef(model1), confint(model1))
    ), 
  2)

## Waiting for profiling to be done...

##                            OR 2.5 % 97.5 %
## (Intercept)              1.26  1.23   1.30
## provs.advanced-masters   0.67  0.62   0.72
## provs.counselors-masters 1.60  1.57   1.64
## provs.psych-masters      0.82  0.80   0.85
## c_provider_1k            0.87  0.86   0.89
## ruralityRural            1.18  1.12   1.25
## any_child                1.14  1.11   1.16
## income_med               0.80  0.78   0.81
## group_pr                 1.26  1.23   1.28
## frac_white_med           1.36  1.33   1.38
## telehealth               1.14  1.12   1.17

Multivariate model, interaction between rurality x providers per 1,000 county population

Interaction term specification

Rurality is a binary variable (rural and non-rural)
providers per 1k county population is continuous and centered

model1int <- glm(insurance_flag ~ provs + c_provider_1k*rurality + any_child 
               + income_med + group_pr + frac_white_med + telehealth, data = data_acs_inc, family = "binomial")
summary(model1int)

## 
## Call:
## glm(formula = insurance_flag ~ provs + c_provider_1k * rurality + 
##     any_child + income_med + group_pr + frac_white_med + telehealth, 
##     family = "binomial", data = data_acs_inc)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  0.233103   0.013129  17.755  < 2e-16 ***
## provs.advanced-masters      -0.400604   0.040744  -9.832  < 2e-16 ***
## provs.counselors-masters     0.472605   0.011471  41.202  < 2e-16 ***
## provs.psych-masters         -0.193717   0.014582 -13.284  < 2e-16 ***
## c_provider_1k               -0.139869   0.008406 -16.640  < 2e-16 ***
## ruralityRural                0.252909   0.034798   7.268 3.65e-13 ***
## any_child                    0.128140   0.010716  11.958  < 2e-16 ***
## income_med                  -0.224733   0.011034 -20.367  < 2e-16 ***
## group_pr                     0.228176   0.010580  21.567  < 2e-16 ***
## frac_white_med               0.304475   0.010405  29.262  < 2e-16 ***
## telehealth                   0.134899   0.011876  11.359  < 2e-16 ***
## c_provider_1k:ruralityRural  0.240498   0.053839   4.467 7.93e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 226890  on 175082  degrees of freedom
## Residual deviance: 221356  on 175071  degrees of freedom
## AIC: 221380
## 
## Number of Fisher Scoring iterations: 4

mod.ORs <- round(
  exp(
    cbind(OR = coef(model1int), confint(model1int))
    ), 
  2)

## Waiting for profiling to be done...

mod.ORs

##                               OR 2.5 % 97.5 %
## (Intercept)                 1.26  1.23   1.30
## provs.advanced-masters      0.67  0.62   0.73
## provs.counselors-masters    1.60  1.57   1.64
## provs.psych-masters         0.82  0.80   0.85
## c_provider_1k               0.87  0.86   0.88
## ruralityRural               1.29  1.20   1.38
## any_child                   1.14  1.11   1.16
## income_med                  0.80  0.78   0.82
## group_pr                    1.26  1.23   1.28
## frac_white_med              1.36  1.33   1.38
## telehealth                  1.14  1.12   1.17
## c_provider_1k:ruralityRural 1.27  1.15   1.41

Contrasts for interaction

est.rural <- gmodels::estimable(model1int,
        c("c_provider_1k" = 1,
          "c_provider_1k:ruralityRural" = 1),
        conf.int = 0.95)

## Registered S3 method overwritten by 'gdata':
##   method         from     
##   reorder.factor DescTools

est.nonrural <- gmodels::estimable(model1int,
        c("c_provider_1k" = 1),
        conf.int = 0.95)

rbind(exp(est.rural[c("Estimate", "Lower.CI", "Upper.CI")]),
      exp(est.nonrural[c("Estimate", "Lower.CI", "Upper.CI")]))

	Estimate	Lower.CI	Upper.CI
(0 0 0 0 1 0 0 0 0 0 0 1)	1.105866	0.9949508	1.2291460
(0 0 0 0 1 0 0 0 0 0 0 0)	0.869472	0.8550921	0.8840936

Interpretation of interaction terms

The first row shows the estimate for rural and the second row shows the estimate for urban

For Rural: For every unit increase from the mean of provider supply per 1000 county population, providers in rural areas have an 11% increase in odds of accepting insurance (95% CI: 1.00, 1.23).
- The mean of providers per 1000 county population is 0.94 providers per 1000 county population

For non-rural: For every unit increase from the mean of provider supply per 1000 county population, providers in non-rural areas have an 13% decrease in odds of accepting insurance (95% CI: 0.86, 0.88).

Visualizing interaction

#non-interaction
plot_model(model1,
           type = "eff",
           terms = c("c_provider_1k", "rurality"),
           title = "No interaction",
           axis.title = c("Provider Supply (per 1000 county population)", 
                          "P(Insurance Acceptance)"),
           legend.title = "Rurality")

# Effect of c_provider_1k at each level of rurality
plot_model(model1int,
           type = "eff",
           terms = c("c_provider_1k", "rurality"),
           title = "Interaction between provider supply and rurality",
           axis.title = c("Provider Supply (per 1000 county population)", 
                          "P(Insurance Acceptance)"),
           legend.title = "Rurality")

Table output of all models

**Insurance Acceptance**

	Dependent variable:

	Insurance Acceptance
	(1)	(2)	(3)

provs.advanced-masters	0.68^***	0.67^***	0.67^***
	(1.82, 2.13)	(1.80, 2.11)	(1.80, 2.12)

provs.counselors-masters	1.55^***	1.60^***	1.60^***
	(4.63, 4.84)	(4.85, 5.08)	(4.86, 5.09)

provs.psych-masters	0.78^***	0.82^***	0.82^***
	(2.13, 2.25)	(2.21, 2.34)	(2.22, 2.35)

c_provider_1k		0.87^***	0.87^***
		(2.36, 2.44)	(2.35, 2.43)

ruralityRural		1.18^***	1.29^***
		(3.08, 3.45)	(3.39, 3.88)

any_child		1.14^***	1.14^***
		(3.05, 3.18)	(3.05, 3.18)

income_med		0.80^***	0.80^***
		(2.17, 2.27)	(2.18, 2.27)

group_pr		1.26^***	1.26^***
		(3.44, 3.59)	(3.44, 3.59)

frac_white_med		1.36^***	1.36^***
		(3.81, 3.96)	(3.80, 3.96)

telehealth		1.14^***	1.14^***
		(3.07, 3.21)	(3.07, 3.21)

c_provider_1k:ruralityRural			1.27^***
			(3.21, 3.96)

Constant	1.65^***	1.26^***	1.26^***
	(5.15, 5.30)	(3.45, 3.63)	(3.44, 3.63)


Observations	175,083	175,083	175,083
Akaike Inf. Crit.	224,203.80	221,398.10	221,379.70

Note:	p<0.1; p<0.05; p<0.01

Clean model output:

Psychology Today

Aine Huntington

Jane Zhu

Friday, 03 May, 2024