PAA 2025 models

#Loading libraries

Code

# Loading libraries

librarian::shelf(knitr, car, foreign, dplyr, janitor, ggplot2, ggpubr, usmapdata, usmap, gtsummary, gt, lme4, emmeans, summarytools, data.table, cobalt, WeightIt, margins, marginaleffects, readxl, nlme, multilevel, MASS, survey)

  The 'cran_repo' argument in shelf() was not set, so it will use
  cran_repo = 'https://cran.r-project.org' by default.

  To avoid this message, set the 'cran_repo' argument to a CRAN
  mirror URL (see https://cran.r-project.org/mirrors.html) or set
  'quiet = TRUE'.

Code

opts_chunk$set(echo = TRUE)

setwd("/Users/briansurratt/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/PAA 2025/R files")

options(scipen = 999) # This turns off scientific notation.

Reading in the data

Importing the cleaned data for late renter households.

Code

latehh <- read.csv("./data/latehh.csv")

nrow(latehh)

[1] 16660

PAA: Probability graphs

Eviction risk by state score (PAA: Do the probability graphs need to be weighted?)

Code

# Calculating statistics for state score so that I can compare actual values to model-based values.

scoreprobs <- latehh %>%
  group_by(score) %>%
  summarize(p=mean(highrisk), n=n())
              
scoreprobs$num <- scoreprobs$p*(1-scoreprobs$p)

scoreprobs$sep <- sqrt(scoreprobs$num/scoreprobs$n)

scoreprobs$me <- 2*scoreprobs$sep
  
ggplot(data = scoreprobs,
       aes(x = score, y = p, ymin=p-me, ymax=p+me)
       ) +
  geom_line() +
  ylim(0,1) +
  geom_ribbon(alpha=0.3,
              aes(color=NULL)
              ) +
labs(title="Eviction Risk by State Score", 
       subtitle="Respondents in Sample", 
       caption="Source: Household Pulse Survey"
     ) +
  xlab(label="State Score") +
  ylab(label="Proportion of Respondents Who Are At High Risk Of Eviction") 

Eviction risk by state score for black vs. non-black

Code

scoreprobs <- latehh %>%
  group_by(nh_black, score) %>%
  summarize(p=mean(highrisk), n=n())

`summarise()` has grouped output by 'nh_black'. You can override using the
`.groups` argument.

Code

scoreprobs$num <- scoreprobs$p*(1-scoreprobs$p)

scoreprobs$sep <- sqrt(scoreprobs$num/scoreprobs$n)

scoreprobs$me <- 2*scoreprobs$sep

scoreprobs <- scoreprobs %>% 
  mutate(black = case_when(nh_black == 1 ~ 'Black',
                           nh_black == 0 ~ 'Non-Black'
                           )
         )

ggplot(scoreprobs, aes(x = score, y = p, group=black)) +
  geom_point() +
  geom_line(aes(color = black)) + 
  labs(x = "Score", y = "Predicted Probability of High Eviction Risk", title = "Predicted Probability of High Eviction Risk for Black versus Non-Black") +
  scale_color_discrete(name="Probability of High Eviction Risk",
                       breaks=c("Black", "Non-Black"),
                       labels=c("Black", "Non-Black")
                       )

Eviction risk by state score for children vs. no children

Code

childprobs <- latehh %>%
  group_by(children, score) %>%
  summarize(p=mean(highrisk), n=n())

`summarise()` has grouped output by 'children'. You can override using the
`.groups` argument.

Code

childprobs$num <- childprobs$p*(1-childprobs$p)

childprobs$sep <- sqrt(childprobs$num/childprobs$n)

childprobs$me <- 2*childprobs$sep

childprobs <- childprobs %>% 
  mutate(child = case_when(children == 1 ~ 'Children',
                           children == 0 ~ 'No Children'
                           )
         )

ggplot(childprobs, aes(x = score, y = p, group = child)) +
  geom_point() +
  geom_line(aes(color = child)) + 
  labs(x = "Score", y = "Predicted Probability of High Eviction Risk", title = "Predicted Probability of High Eviction Risk Children vs. No children") +
  scale_color_discrete(name="Probability of High Eviction Risk",
                       breaks=c("Children", "No Children"),
                       labels=c("Children", "No Children")
                       )

Eviction risk by age

Code

# Calculating statistics for specific ages so that I can compare actual values to model-based values.

ageprobs<- latehh %>%
    group_by(age) %>%
    summarize(p=mean(highrisk),n=n())
              
ageprobs$num <- ageprobs$p*(1-ageprobs$p)

ageprobs$sep <- sqrt(ageprobs$num/ageprobs$n)

ageprobs$me <- 2*ageprobs$sep
  
ggplot(data =ageprobs, aes(x = age, y = p, ymin=p-me, ymax=p+me)) +
     geom_line() + ylim(0,1) + geom_ribbon(alpha=0.3,aes(color=NULL)) +
labs(title="Eviction Risk by Age", 
       subtitle="Respondents in Sample", 
       caption="Source: Household Pulse Survey") +
       xlab(label="Age at Survey") +
  ylab(label="Proportion of Respondents Who Are At High Risk Of Eviction") 

Eviction risk by race/eth

Code

# Calculating statistics for races so that I can compare actual values to model-based values.

raceprobs<- latehh %>%
    group_by(race_eth) %>%
    summarize(p=mean(highrisk), n=n())
              
raceprobs$num <- raceprobs$p*(1-raceprobs$p)

raceprobs$sep <- sqrt(raceprobs$num/raceprobs$n)

raceprobs$me <- 2*raceprobs$sep
  
ggplot(data =raceprobs, aes(x = race_eth, y = p, ymin=p-me, ymax=p+me)) +
  geom_boxplot() + 
  ylim(0,1) + 
  geom_ribbon(alpha=0.3,aes(color=NULL)) +
labs(title="Eviction Risk by Race/Ethnicity", 
       subtitle="Respondents in Sample", 
       caption="Source: Household Pulse Survey") +
       xlab(label="Race/Ethnicity") +
  ylab(label="Proportion of Respondents Who Are At High Risk Of Eviction") 

Eviction risk by income level

Not using this code since I formatted the inclvl variable differently

Calculating statistics for inclvl so that I can compare actual values to model-based values.

latehhinclvl <- factor(latehhinclvl, levels=c(“12500”, “30000”, “42500”, “62500”, “87500”, “125000”, “175000”, “300000”))

latehhinclvl <- as.numeric(latehhinclvl)

incprobs<- latehh %>% group_by(inclvl) %>% summarize(p=mean(highrisk), n=n())

incprobsnum <- incprobsp*(1-incprobs$p)

incprobssep <- sqrt(incprobsnum/incprobs$n)

incprobsme <- 2*incprobssep

ggplot(data = incprobs, aes(x = inclvl, y = p, ymin=p-me, ymax=p+me)) + geom_line() + ylim(0,1) + geom_ribbon(alpha=0.3,aes(color=NULL)) + labs(title=“Eviction Risk by Income Level”, subtitle=“Respondents in Sample”, caption=“Source: Household Pulse Survey”) + xlab(label=“Income Level”) + ylab(label=“Proportion of Respondents Who Are At High Risk Of Eviction”)

Eviction risk by number of children

Code

# Calculating statistics for number of children so that I can compare actual values to model-based values.

kidsprobs<- latehh %>%
    group_by(thhld_numkid) %>%
    summarize(p=mean(highrisk), n=n())
              
kidsprobs$num <- kidsprobs$p*(1-kidsprobs$p)

kidsprobs$sep <- sqrt(kidsprobs$num/kidsprobs$n)

kidsprobs$me <- 2*kidsprobs$sep
  
ggplot(data =kidsprobs, aes(x = thhld_numkid, y = p, ymin=p-me, ymax=p+me)) +
  geom_line() + ylim(0,1) + geom_ribbon(alpha=0.3,aes(color=NULL)) +
  labs(title="Eviction Risk by Number of Children",
       subtitle="Respondents in Sample",
       caption="Source: Household Pulse Survey") +
  xlab(label="Number of Children") +
  ylab(label="Proportion of Respondents Who Are At High Risk Of Eviction") 

Scatterplot of score vs highrisk

I need to make a table that groups by state score and finds a mean or median of the evictrisk variable.

Making a table grouped by state score with mean evictrisk.

“Evictrisk” is defined as “How likely is it that your household will have to leave this home or apartment within the next two months because of eviction? Select only one answer.” 1 is not likely at all, 2 is not very likely, 3 is somewhat likely, and 4 is very likely.

Code

scoreplot <- latehh %>%
  group_by(state) %>%
  summarise(score = mean(score), highrisk = mean(highrisk))

scoreplot

# A tibble: 51 × 3
   state                score highrisk
   <chr>                <dbl>    <dbl>
 1 Alabama               0.83    0.5  
 2 Alaska                1.23    0.441
 3 Arizona               0.23    0.476
 4 Arkansas              0       0.483
 5 California            1.33    0.406
 6 Colorado              3.38    0.429
 7 Connecticut           3.78    0.365
 8 Delaware              3.88    0.375
 9 District of Columbia  4.63    0.314
10 Florida               1.58    0.439
# ℹ 41 more rows

ggplot2 scatter plot

Code

risk_by_score_plot <- ggplot(scoreplot, aes(x=score, y=highrisk)) +
  geom_point() +
  geom_smooth(method="lm", se=TRUE, fullrange=FALSE, level=0.95) +
  ylim(0,1) +
  labs(title = "Mean Perceived Likelihood of Eviction by State Score",
       subtitle = "Renters Who Are Not Current On Rent", 
       x="State Score", 
       y = "Mean Perceived Likelihood of Eviction",
       caption = "Source: HH Pulse and COVID State Housing Policy Score") +
   stat_cor(method = "pearson", label.x = 2, label.y = 0.65) +
   theme(
    plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 15),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15)
  )
   #theme_classic()

  #labs(title="Perceived Eviction Risk by State Score", x="State Score", y = "Perceived Eviction Risk") +
  #stat_regline_equation(label.x=2, label.y=3.5) + # This adds the regression equation
  #stat_cor(aes(label=after_stat(rr.label)), label.x=2, label.y=3.2) + # This adds the r squared
 

risk_by_score_plot

`geom_smooth()` using formula = 'y ~ x'

Code

ggsave("risk_by_score_plot.png", path = "./figures/")

Saving 7 x 5 in image
`geom_smooth()` using formula = 'y ~ x'

Matching

Create treatment categories

Let’s check the distribution of scores

Code

descr(latehh$score)

Descriptive Statistics  
latehh$score  
N: 16660  

                       score
----------------- ----------
             Mean       1.91
          Std.Dev       1.32
              Min       0.00
               Q1       0.83
           Median       1.43
               Q3       3.13
              Max       4.63
              MAD       1.56
              IQR       2.30
               CV       0.69
         Skewness       0.28
      SE.Skewness       0.02
         Kurtosis      -1.13
          N.Valid   16660.00
                N   16660.00
        Pct.Valid     100.00

Code

boxplot(latehh$score)

Let’s make 3 treatment groups based on terciles.

Code

quantile(latehh$score, probs = 0:3/3)

       0% 33.33333% 66.66667%      100% 
     0.00      1.28      2.73      4.63 

Let’s create 3 score groups. Low protection is 0 to <1.28, medium protection is 1.28 to < 2.73, and high protection is 2.73 to 5.

Code

latehh <- latehh %>% 
  mutate(treatment = case_when(score < 1.28 ~ 'low',
                              score >= 1.28 & score < 2.73 ~ 'medium',
                              score >=2.73 ~ 'high'
                              )
         )

latehh$treatment <- factor(latehh$treatment, levels=c('low', 'medium', 'high'))

Code

tabyl(latehh$treatment)

 latehh$treatment    n   percent
              low 5136 0.3082833
           medium 5713 0.3429172
             high 5811 0.3487995

PAA adjust income and education before matching

What is the median income level in latehh?

Code

tabyl(latehh$income)

 latehh$income    n     percent
             1 6894 0.413805522
             2 3551 0.213145258
             3 2607 0.156482593
             4 1999 0.119987995
             5  829 0.049759904
             6  521 0.031272509
             7  126 0.007563025
             8  133 0.007983193

Code

class(latehh$income)

[1] "integer"

Code

hist(latehh$income)

Create dummy variables for three income categories

INCOME: Total household income (before taxes). inclvl1) Less than $25,000 inclvl2) $25,000 - $49,999 inclvl3) $50,000 and above

Code

latehh <- latehh %>% 
  mutate(inclvl1 = case_when(income == 1 ~ 1,
                               TRUE ~ 0))

tabyl(latehh$inclvl1)

 latehh$inclvl1    n   percent
              0 9766 0.5861945
              1 6894 0.4138055

Code

latehh <- latehh %>% 
  mutate(inclvl2 = case_when(income %in% c(2,3) ~ 1,
                               TRUE ~ 0))

tabyl(latehh$inclvl2)

 latehh$inclvl2     n   percent
              0 10502 0.6303721
              1  6158 0.3696279

Code

latehh <- latehh %>% 
  mutate(inclvl3 = case_when(income %in% c(4:8) ~ 1,
                               TRUE ~ 0))

tabyl(latehh$inclvl3)

 latehh$inclvl3     n   percent
              0 13052 0.7834334
              1  3608 0.2165666

Three education categories:

Code

tabyl(latehh$hsorless)

 latehh$hsorless     n   percent
               0 12001 0.7203481
               1  4659 0.2796519

Code

tabyl(latehh$somecollege)

 latehh$somecollege    n   percent
                  0 9304 0.5584634
                  1 7356 0.4415366

Code

tabyl(latehh$baorhigher)

 latehh$baorhigher     n   percent
                 0 12015 0.7211885
                 1  4645 0.2788115

Code

nrow(latehh)

[1] 16660

Balance before matching

Code

bal.tab(treatment ~ age +
          female + 
          nh_white +
          hispanic +
          nh_black +
          nh_asian +
          other +
          inclvl1 +
          inclvl2 +
          inclvl3 +
          hsorless +
          somecollege +
          baorhigher +
          unmarried +
          children,
        data = latehh,
        focal = "low",
        estimand = "ATT",
        thresholds = c(m = .05)
        )

Note: `s.d.denom` not specified; assuming "pooled".

Balance summary across all treatment pairs
               Type Max.Diff.Un      M.Threshold.Un
age         Contin.      0.0937 Not Balanced, >0.05
female       Binary      0.0547 Not Balanced, >0.05
nh_white     Binary      0.1158 Not Balanced, >0.05
hispanic     Binary      0.1223 Not Balanced, >0.05
nh_black     Binary      0.0371     Balanced, <0.05
nh_asian     Binary      0.0507 Not Balanced, >0.05
other        Binary      0.0360     Balanced, <0.05
inclvl1      Binary      0.0642 Not Balanced, >0.05
inclvl2      Binary      0.0053     Balanced, <0.05
inclvl3      Binary      0.0605 Not Balanced, >0.05
hsorless     Binary      0.0215     Balanced, <0.05
somecollege  Binary      0.0686 Not Balanced, >0.05
baorhigher   Binary      0.0901 Not Balanced, >0.05
unmarried    Binary      0.0569 Not Balanced, >0.05
children     Binary      0.0722 Not Balanced, >0.05

Balance tally for mean differences
                    count
Balanced, <0.05         4
Not Balanced, >0.05    11

Variable with the greatest mean difference
 Variable Max.Diff.Un      M.Threshold.Un
 hispanic      0.1223 Not Balanced, >0.05

Sample sizes
    medium high  low
All   5713 5811 5136

Weighting

Code

w.out <- weightit(treatment ~ age +
                    female +
                    nh_white +
                    hispanic +
                    nh_black +
                    nh_asian +
                    other +
                    inclvl1 +
                    inclvl2 +
                    inclvl3 +
                    hsorless +
                    somecollege +
                    baorhigher +
                    unmarried +
                    children,
                  data = latehh,
                  focal = "low",
                  estimand = "ATT",
                  method="GLM")

Code

summary(w.out)

                  Summary of weights

- Weight ranges:

          Min                                  Max
low    1.0000            ||                 1.0000
medium 0.1972 |---------------------------| 2.2091
high   0.2583  |--------------------|       1.8098

- Units with the 5 most extreme weights by group:
                                          
             5      4      3      2      1
    low      1      1      1      1      1
         13343   7217  14257   7078   7390
 medium 2.1634 2.1634 2.1785 2.1938 2.2091
          5484   8350   6016   8377   9494
   high 1.7549  1.761 1.7678 1.8029 1.8098

- Weight statistics:

       Coef of Var   MAD Entropy # Zeros
low          0.000 0.000   0.000       0
medium       0.463 0.396   0.111       0
high         0.332 0.270   0.056       0

- Effective Sample Sizes:

            low  medium    high
Unweighted 5136 5713.   5811.  
Weighted   5136 4706.12 5235.66

Balance after matching

Code

bal.tab(w.out, stats = c("m", "v"), thresholds = c(m = .05))

Balance summary across all treatment pairs
               Type Max.Diff.Adj     M.Threshold Max.V.Ratio.Adj
age         Contin.       0.0155 Balanced, <0.05          1.0429
female       Binary       0.0040 Balanced, <0.05               .
nh_white     Binary       0.0006 Balanced, <0.05               .
hispanic     Binary       0.0008 Balanced, <0.05               .
nh_black     Binary       0.0014 Balanced, <0.05               .
nh_asian     Binary       0.0003 Balanced, <0.05               .
other        Binary       0.0014 Balanced, <0.05               .
inclvl1      Binary       0.0049 Balanced, <0.05               .
inclvl2      Binary       0.0052 Balanced, <0.05               .
inclvl3      Binary       0.0005 Balanced, <0.05               .
hsorless     Binary       0.0033 Balanced, <0.05               .
somecollege  Binary       0.0084 Balanced, <0.05               .
baorhigher   Binary       0.0051 Balanced, <0.05               .
unmarried    Binary       0.0015 Balanced, <0.05               .
children     Binary       0.0011 Balanced, <0.05               .

Balance tally for mean differences
                    count
Balanced, <0.05        15
Not Balanced, >0.05     0

Variable with the greatest mean difference
 Variable Max.Diff.Adj     M.Threshold
      age       0.0155 Balanced, <0.05

Effective sample sizes
            medium    high  low
Unadjusted 5713.   5811.   5136
Adjusted   4706.12 5235.66 5136

Balance plots

Age

Code

bal.plot(w.out, var.name = "age", which = "both")

Female

Code

bal.plot(w.out, var.name = "female", which = "both", type = "histogram")

nh_white

Code

bal.plot(w.out, var.name = "nh_white", which = "both", type = "histogram")

nh_black

Code

bal.plot(w.out, var.name = "nh_black", which = "both", type = "histogram")

hispanic

Code

bal.plot(w.out, var.name = "hispanic", which = "both", type = "histogram")

nh_asian

Code

bal.plot(w.out, var.name = "nh_asian", which = "both", type = "histogram")

other

Code

bal.plot(w.out, var.name = "other", which = "both", type = "histogram")

inclvl1

Code

bal.plot(w.out, var.name = "inclvl1",
         which = "both",
         type = "histogram",
         mirror = TRUE)

###inclvl2

Code

bal.plot(w.out, var.name = "inclvl2",
         which = "both",
         type = "histogram",
         mirror = TRUE)

inclvl3

Code

bal.plot(w.out, var.name = "inclvl3",
         which = "both",
         type = "histogram",
         mirror = TRUE)

hsorless

Code

bal.plot(w.out, var.name = "hsorless",
         which = "both",
         type = "histogram",
         mirror = TRUE)

somecollege

Code

bal.plot(w.out, var.name = "somecollege",
         which = "both",
         type = "histogram",
         mirror = TRUE)

baorhigher

Code

bal.plot(w.out, var.name = "baorhigher",
         which = "both",
         type = "histogram",
         mirror = TRUE)

unmarried

Code

bal.plot(w.out, var.name = "unmarried", which = "both", type = "histogram")

children

Code

bal.plot(w.out, var.name = "children", which = "both", type = "histogram")

Adding weights to latehh

Code

latehh$weights <- w.out$weights

Code

# Export to csv files

fwrite(latehh, "./data/latehh.csv")

Filter and write csv files for low, medium, and high state scores

Filter for low state score late households

Code

low_latehh <- latehh %>% 
  filter(treatment == "low")

nrow(low_latehh)

[1] 5136

Write out csv for low state score late households

Code

# Export to csv files

fwrite(low_latehh, "./data/low_latehh.csv")

Filter for medium state score late households

Code

medium_latehh <- latehh %>% 
  filter(treatment == "medium")

nrow(medium_latehh)

[1] 5713

Write out csv for medium state score late households

Code

# Export to csv files

fwrite(medium_latehh, "./data/medium_latehh.csv")

Filter for high state score late households

Code

high_latehh <- latehh %>% 
  filter(treatment == "high")

nrow(high_latehh)

[1] 5811

Write out csv for high state score late households

Code

# Export to csv files

fwrite(high_latehh, "./data/high_latehh.csv")

G-computation for the treatment effect

No weight

Model 1: glmer, no controls (glmer, latehh, highrisk, random intercept + score)

Code

m1 <- glmer(highrisk ~ score +
              (1 | abbv),
            data = latehh,
            family = binomial,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10)

summary(m1)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + (1 | abbv)
   Data: latehh
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  22584.0   22607.1  -11289.0   22578.0     16657 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.0746 -0.8437 -0.7505  1.1238  1.4348 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03886  0.1971  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept) -0.11875    0.05485  -2.165    0.0304 *  
score       -0.09851    0.02363  -4.169 0.0000306 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
      (Intr)
score -0.802

Model 2: gmler with controls (glmer, latehh, highrisk, random intercept + all covariates)

Code

latehh$race_eth <- factor(latehh$race_eth, levels=c("nh_white",  "nh_black",  "hispanic", "nh_asian", "other"))
latehh$inclvl <- factor(latehh$inclvl, levels=c("25to49",  "lessthan25",  "50andover"))
latehh$educ <- factor(latehh$educ, levels=c("somecollege",  "HSorless",  "BAorhigher"))

Code

m2 <- glmer(highrisk ~ score +
              age +
              female +
              race_eth +
              unmarried +
              children +
              inclvl +
              educ +
              (1 | abbv),
            data = latehh,
            family = binomial,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m2)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + age + female + race_eth + unmarried + children +  
    inclvl + educ + (1 | abbv)
   Data: latehh
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  21654.9   21770.7  -10812.5   21624.9     16645 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5797 -0.8763 -0.5673  1.0066  3.4430 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03012  0.1735  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                  Estimate Std. Error z value             Pr(>|z|)    
(Intercept)      -0.132604   0.095814  -1.384              0.16636    
score            -0.070952   0.021847  -3.248              0.00116 ** 
age              -0.004196   0.001312  -3.198              0.00139 ** 
female           -0.056061   0.035904  -1.561              0.11843    
race_ethnh_black  0.239276   0.042895   5.578     0.00000002430467 ***
race_ethhispanic -0.142971   0.047487  -3.011              0.00261 ** 
race_ethnh_asian -1.154769   0.091585 -12.609 < 0.0000000000000002 ***
race_ethother     0.075254   0.065623   1.147              0.25148    
unmarried         0.224559   0.038067   5.899     0.00000000365597 ***
children          0.177097   0.035179   5.034     0.00000047979001 ***
inclvllessthan25  0.330794   0.036716   9.010 < 0.0000000000000002 ***
inclvl50andover  -0.403349   0.047477  -8.496 < 0.0000000000000002 ***
educHSorless     -0.036374   0.041843  -0.869              0.38469    
educBAorhigher   -0.279035   0.039551  -7.055     0.00000000000173 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 14 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Model 2: Marginal effects

Code

#average marginal effect
mfx2<-margins(m2,variables=c("score",
                "age",
                "female",
                "race_eth",
                "unmarried",
                "children",
                "inclvl",
                "educ"))

summary(mfx2)

           factor     AME     SE        z      p   lower   upper
              age -0.0010 0.0003  -3.2019 0.0014 -0.0015 -0.0004
         children  0.0405 0.0080   5.0492 0.0000  0.0248  0.0563
   educBAorhigher -0.0642 0.0091  -7.0340 0.0000 -0.0821 -0.0463
     educHSorless -0.0085 0.0098  -0.8695 0.3846 -0.0276  0.0106
           female -0.0128 0.0082  -1.5618 0.1183 -0.0289  0.0033
  inclvl50andover -0.0900 0.0104  -8.6555 0.0000 -0.1103 -0.0696
 inclvllessthan25  0.0784 0.0087   9.0364 0.0000  0.0614  0.0954
 race_ethhispanic -0.0333 0.0110  -3.0316 0.0024 -0.0548 -0.0118
 race_ethnh_asian -0.2319 0.0150 -15.4725 0.0000 -0.2613 -0.2025
 race_ethnh_black  0.0569 0.0102   5.5646 0.0000  0.0368  0.0769
    race_ethother  0.0178 0.0155   1.1434 0.2529 -0.0127  0.0482
            score -0.0162 0.0050  -3.2638 0.0011 -0.0260 -0.0065
        unmarried  0.0514 0.0087   5.9209 0.0000  0.0344  0.0684

Model 3: glmer with interactions

Code

m3 <- glmer(highrisk ~
              score*(age +
                       female +
                       race_eth +
                       unmarried +
                       children +
                       inclvl +
                       educ
                     ) +
              (1 | abbv),
            data = latehh,
            family = binomial,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0918918 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m3)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score * (age + female + race_eth + unmarried + children +  
    inclvl + educ) + (1 | abbv)
   Data: latehh
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  21664.2   21872.7  -10805.1   21610.2     16633 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5786 -0.8786 -0.5643  1.0056  3.5169 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.02914  0.1707  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                         Estimate Std. Error z value        Pr(>|z|)    
(Intercept)            -0.0073147  0.1504944  -0.049        0.961235    
score                  -0.1339348  0.0653887  -2.048        0.040532 *  
age                    -0.0053762  0.0023150  -2.322        0.020217 *  
female                 -0.0436282  0.0630840  -0.692        0.489196    
race_ethnh_black        0.3191543  0.0728967   4.378 0.0000119679326 ***
race_ethhispanic       -0.2787083  0.0881251  -3.163        0.001563 ** 
race_ethnh_asian       -1.1728288  0.1725489  -6.797 0.0000000000107 ***
race_ethother           0.0471236  0.1107576   0.425        0.670497    
unmarried               0.2092401  0.0660085   3.170        0.001525 ** 
children                0.1138623  0.0607910   1.873        0.061067 .  
inclvllessthan25        0.2630916  0.0635814   4.138 0.0000350545754 ***
inclvl50andover        -0.4678613  0.0845529  -5.533 0.0000000314163 ***
educHSorless           -0.0823904  0.0722974  -1.140        0.254452    
educBAorhigher         -0.2360506  0.0686642  -3.438        0.000587 ***
score:age               0.0005999  0.0009885   0.607        0.543946    
score:female           -0.0080491  0.0272238  -0.296        0.767487    
score:race_ethnh_black -0.0441875  0.0319096  -1.385        0.166123    
score:race_ethhispanic  0.0670928  0.0378072   1.775        0.075963 .  
score:race_ethnh_asian  0.0079619  0.0733933   0.108        0.913613    
score:race_ethother     0.0154463  0.0487544   0.317        0.751381    
score:unmarried         0.0073283  0.0292434   0.251        0.802127    
score:children          0.0329208  0.0265669   1.239        0.215283    
score:inclvllessthan25  0.0367870  0.0278164   1.322        0.186004    
score:inclvl50andover   0.0347179  0.0361240   0.961        0.336514    
score:educHSorless      0.0245696  0.0315885   0.778        0.436685    
score:educBAorhigher   -0.0230095  0.0300297  -0.766        0.443543    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 26 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0918918 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

With weight

Model 4: glmer, no controls

Code

m4 <- glmer(highrisk ~ score +
              (1 | abbv),
            data = latehh,
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10)

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Code

summary(m4)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  21076.7   21099.9  -10535.4   21070.7     16657 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4796 -0.8801 -0.5908  1.0322  1.7949 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.04122  0.203   
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.10947    0.05610  -1.951  0.05101 . 
score       -0.06791    0.02430  -2.795  0.00519 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
      (Intr)
score -0.801

Code

avg_comparisons(m4, variables = "score",
                vcov = "HC3",
                newdata = subset(latehh, treatment == "low"))

Warning: Unable to extract a variance-covariance matrix using this `vcov`
  argument. Standard errors are computed using the default variance
  instead. Perhaps the model or argument is not supported by the
  `sandwich` ('HC0', 'HC3', ~clusterid, etc.) or `clubSandwich` ('CR0',
  etc.) packages. If you believe that the model is supported by one of
  these two packages, you can open a feature request on Github.

 Estimate Std. Error     z Pr(>|z|)   S   2.5 %   97.5 %
  -0.0167    0.00601 -2.78  0.00536 7.5 -0.0285 -0.00496

Term: score
Type:  response 
Comparison: +1

Code

avg_predictions(m4, variables = "score",
                vcov = "HC3",
                newdata = subset(latehh, treatment == "low"))

Warning: Unable to extract a variance-covariance matrix using this `vcov`
  argument. Standard errors are computed using the default variance
  instead. Perhaps the model or argument is not supported by the
  `sandwich` ('HC0', 'HC3', ~clusterid, etc.) or `clubSandwich` ('CR0',
  etc.) packages. If you believe that the model is supported by one of
  these two packages, you can open a feature request on Github.

 score Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
  0.00    0.475    0.01390 34.2   <0.001 848.9 0.448  0.503
  0.83    0.461    0.01031 44.8   <0.001   Inf 0.441  0.482
  1.43    0.451    0.00864 52.2   <0.001   Inf 0.434  0.468
  3.13    0.423    0.01111 38.1   <0.001   Inf 0.401  0.445
  4.63    0.399    0.01798 22.2   <0.001 359.2 0.363  0.434

Type:  response 

Model 5: glmer with controls

Code

m5 <- glmer(highrisk ~ score +
              age +
              female +
              race_eth +
              unmarried +
              children +
              inclvl +
              educ +
              (1 | abbv),
            data = latehh,
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00334097 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m5)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + age + female + race_eth + unmarried + children +  
    inclvl + educ + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20431.4   20547.2  -10200.7   20401.4     16645 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9531 -0.8487 -0.4304  0.9925  2.7104 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.0351   0.1874  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                  Estimate Std. Error z value             Pr(>|z|)    
(Intercept)      -0.182648   0.099508  -1.836              0.06643 .  
score            -0.071638   0.023116  -3.099              0.00194 ** 
age              -0.003661   0.001378  -2.657              0.00788 ** 
female           -0.050496   0.037766  -1.337              0.18120    
race_ethnh_black  0.228948   0.042727   5.358    0.000000083972952 ***
race_ethhispanic -0.155000   0.055834  -2.776              0.00550 ** 
race_ethnh_asian -1.057335   0.118431  -8.928 < 0.0000000000000002 ***
race_ethother     0.070547   0.061376   1.149              0.25038    
unmarried         0.218908   0.039324   5.567    0.000000025944897 ***
children          0.180676   0.036108   5.004    0.000000562395930 ***
inclvllessthan25  0.342887   0.037269   9.200 < 0.0000000000000002 ***
inclvl50andover  -0.369069   0.051694  -7.139    0.000000000000937 ***
educHSorless     -0.022897   0.041941  -0.546              0.58512    
educBAorhigher   -0.248847   0.040536  -6.139    0.000000000831026 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 14 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00334097 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

avg_comparisons(m5, variables = "score",
                vcov = "HC3",
                newdata = subset(latehh, treatment == "low"))

Warning: Unable to extract a variance-covariance matrix using this `vcov`
  argument. Standard errors are computed using the default variance
  instead. Perhaps the model or argument is not supported by the
  `sandwich` ('HC0', 'HC3', ~clusterid, etc.) or `clubSandwich` ('CR0',
  etc.) packages. If you believe that the model is supported by one of
  these two packages, you can open a feature request on Github.

 Estimate Std. Error     z Pr(>|z|)   S   2.5 %   97.5 %
  -0.0169    0.00547 -3.09    0.002 9.0 -0.0276 -0.00618

Term: score
Type:  response 
Comparison: +1

Code

avg_predictions(m5, variables = "score",
                vcov = "HC3",
                newdata = subset(latehh, treatment == "low"))

Warning: Unable to extract a variance-covariance matrix using this `vcov`
  argument. Standard errors are computed using the default variance
  instead. Perhaps the model or argument is not supported by the
  `sandwich` ('HC0', 'HC3', ~clusterid, etc.) or `clubSandwich` ('CR0',
  etc.) packages. If you believe that the model is supported by one of
  these two packages, you can open a feature request on Github.

 score Estimate Std. Error    z Pr(>|z|)      S 2.5 % 97.5 %
  0.00    0.475    0.01266 37.5   <0.001 1019.9 0.450  0.500
  0.83    0.461    0.00940 49.0   <0.001    Inf 0.442  0.479
  1.43    0.451    0.00788 57.2   <0.001    Inf 0.435  0.466
  3.13    0.422    0.01013 41.7   <0.001    Inf 0.402  0.442
  4.63    0.397    0.01640 24.2   <0.001  428.4 0.365  0.429

Type:  response 

Marginal effects

Code

#average marginal effect
mfxm5<-margins(m5,variables=c("score",
                "children",
                "age",
                "female",
                "race_eth",
                "inclvl",
                "educ",
                "unmarried"))

summary(mfxm5)

           factor     AME     SE        z      p   lower   upper
              age -0.0008 0.0003  -2.6604 0.0078 -0.0015 -0.0002
         children  0.0415 0.0083   5.0223 0.0000  0.0253  0.0577
   educBAorhigher -0.0575 0.0094  -6.1277 0.0000 -0.0758 -0.0391
     educHSorless -0.0054 0.0098  -0.5460 0.5850 -0.0246  0.0139
           female -0.0116 0.0087  -1.3372 0.1812 -0.0286  0.0054
  inclvl50andover -0.0827 0.0114  -7.2779 0.0000 -0.1050 -0.0605
 inclvllessthan25  0.0815 0.0088   9.2355 0.0000  0.0642  0.0988
 race_ethhispanic -0.0361 0.0129  -2.8020 0.0051 -0.0614 -0.0109
 race_ethnh_asian -0.2171 0.0200 -10.8445 0.0000 -0.2564 -0.1779
 race_ethnh_black  0.0545 0.0102   5.3441 0.0000  0.0345  0.0745
    race_ethother  0.0167 0.0146   1.1464 0.2516 -0.0118  0.0452
            score -0.0165 0.0053  -3.1169 0.0018 -0.0268 -0.0061
        unmarried  0.0503 0.0090   5.5838 0.0000  0.0327  0.0680

Code

exp(-0.0165*5)

[1] 0.9208114

Model 6: glmer with interactions

Code

m6 <- glmer(highrisk ~ 
             score*(age +
                      female +
                      race_eth +
                      unmarried +
                      children +
                      inclvl +
                      educ
                    ) +
              (1 | abbv),
            data = latehh,
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0222058 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m6)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score * (age + female + race_eth + unmarried + children +  
    inclvl + educ) + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20443.0   20651.5  -10194.5   20389.0     16633 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9298 -0.8521 -0.4294  0.9889  2.7100 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03401  0.1844  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                         Estimate Std. Error z value    Pr(>|z|)    
(Intercept)            -0.1176147  0.1536464  -0.765     0.44398    
score                  -0.1049843  0.0678267  -1.548     0.12166    
age                    -0.0034417  0.0023730  -1.450     0.14696    
female                 -0.0480480  0.0647979  -0.742     0.45839    
race_ethnh_black        0.3295545  0.0729087   4.520 0.000006181 ***
race_ethhispanic       -0.2564879  0.0968912  -2.647     0.00812 ** 
race_ethnh_asian       -1.0527542  0.1991058  -5.287 0.000000124 ***
race_ethother           0.0531927  0.1071426   0.496     0.61957    
unmarried               0.2140133  0.0674982   3.171     0.00152 ** 
children                0.1314814  0.0617562   2.129     0.03325 *  
inclvllessthan25        0.2616437  0.0642639   4.071 0.000046732 ***
inclvl50andover        -0.4412270  0.0882632  -4.999 0.000000576 ***
educHSorless           -0.0775051  0.0726711  -1.067     0.28619    
educBAorhigher         -0.2391265  0.0696852  -3.432     0.00060 ***
score:age              -0.0001088  0.0010383  -0.105     0.91653    
score:female           -0.0028230  0.0286096  -0.099     0.92140    
score:race_ethnh_black -0.0546687  0.0319535  -1.711     0.08710 .  
score:race_ethhispanic  0.0534208  0.0428324   1.247     0.21232    
score:race_ethnh_asian -0.0033738  0.0907346  -0.037     0.97034    
score:race_ethother     0.0095139  0.0476226   0.200     0.84166    
score:unmarried         0.0016829  0.0298261   0.056     0.95501    
score:children          0.0262735  0.0272918   0.963     0.33570    
score:inclvllessthan25  0.0444820  0.0282809   1.573     0.11575    
score:inclvl50andover   0.0400459  0.0391690   1.022     0.30660    
score:educHSorless      0.0295389  0.0318422   0.928     0.35358    
score:educBAorhigher   -0.0060617  0.0308375  -0.197     0.84416    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 26 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0222058 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

avg_comparisons(m6, variables = "score",
                vcov = "HC3",
                newdata = subset(latehh, treatment == "low"))

Warning: Unable to extract a variance-covariance matrix using this `vcov`
  argument. Standard errors are computed using the default variance
  instead. Perhaps the model or argument is not supported by the
  `sandwich` ('HC0', 'HC3', ~clusterid, etc.) or `clubSandwich` ('CR0',
  etc.) packages. If you believe that the model is supported by one of
  these two packages, you can open a feature request on Github.

 Estimate Std. Error     z Pr(>|z|)   S   2.5 %   97.5 %
  -0.0166    0.00541 -3.06  0.00222 8.8 -0.0272 -0.00595

Term: score
Type:  response 
Comparison: +1

Code

avg_predictions(m6, variables = "score",
                vcov = "HC3",
                newdata = subset(latehh, treatment == "low"))

Warning: Unable to extract a variance-covariance matrix using this `vcov`
  argument. Standard errors are computed using the default variance
  instead. Perhaps the model or argument is not supported by the
  `sandwich` ('HC0', 'HC3', ~clusterid, etc.) or `clubSandwich` ('CR0',
  etc.) packages. If you believe that the model is supported by one of
  these two packages, you can open a feature request on Github.

 score Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
  0.00    0.474     0.0125 37.8   <0.001   Inf 0.450  0.499
  0.83    0.460     0.0093 49.5   <0.001   Inf 0.442  0.479
  1.43    0.450     0.0078 57.7   <0.001   Inf 0.435  0.466
  3.13    0.422     0.0100 42.2   <0.001   Inf 0.403  0.442
  4.63    0.399     0.0162 24.6   <0.001 442.2 0.367  0.430

Type:  response 

Predicted probability graphs

Plots for Model 2: With controls, without weight

Plot of Model 2: High Eviction Risk by State Score

This model has controls, but no weights, no interaction

Code

plot_m2<- plot_predictions(m2, condition = c("score")) +
  labs(title = "Predicted Probability of High Eviction Risk by State Score",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Score",
                       breaks=c(),
                       labels=c("Score")
                       ) +
  scale_fill_discrete(name="Score",
                       breaks=c(),
                       labels=c("Score")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('plot_m2.png')

Code

plot_m2

Plot of Model 2: Race/Ethnicity

This model has controls, but no weights, no interaction

Code

plot_race_ethm2 <- plot_predictions(m2, condition = c("score", "race_eth")) +
  labs(title = "Predicted Probability of High Eviction Risk, by Race/Ethnicity",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Race/Ethnicity",
                       breaks=c("nh_white", "nh_black", "hispanic", "nh_asian", "other"),
                       labels=c("White", "Black", "Hispanic", "Asian", "Other")
                       ) +
  scale_fill_discrete(name="Race/Ethnicity",
                       breaks=c("nh_white", "nh_black", "hispanic", "nh_asian", "other"),
                       labels=c("White", "Black", "Hispanic", "Asian", "Other")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('race_ethm2.png')

Code

plot_race_ethm2

Plot of Model 2: Income

This model has controls, but no weights, no interaction

Code

tabyl(latehh$inclvl)

 latehh$inclvl    n   percent
        25to49 6158 0.3696279
    lessthan25 6894 0.4138055
     50andover 3608 0.2165666

Code

plot_inclvlm2 <- plot_predictions(m2, condition = c("score", "inclvl")) +
  labs(title = "Predicted Probability of High Eviction Risk, Income Level",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Income Level",
                       breaks=c("lessthan25", "25to49", "50andover"),
                       labels=c("Less than $25K", "$25K to $49K", "$50K and over")
                       ) +
  scale_fill_discrete(name="Income Level",
                      breaks=c("lessthan25", "25to49", "50andover"),
                      labels=c("Less than $25K", "$25K to $49K", "$50K and over")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('plot_inclvlm2.png')

Code

plot_inclvlm2

Plot of Model 2: Education

This model has controls, but no weights, no interaction

Code

tabyl(latehh$educ)

 latehh$educ    n   percent
 somecollege 5230 0.3139256
    HSorless 4659 0.2796519
  BAorhigher 6771 0.4064226

Code

plot_educm2 <- plot_predictions(m2, condition = c("score", "educ")) +
  labs(title = "Predicted Probability of High Eviction Risk, Education Level",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Education Level",
                       breaks=c("HSorless", "somecollege", "BAorhigher"),
                       labels=c("Highschool or less", "Some college", "BA or Higher")
                       ) +
  scale_fill_discrete(name="Education Level",
                      breaks=c("HSorless", "somecollege", "BAorhigher"),
                      labels=c("Highschool or less", "Some college", "BA or Higher")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('educm2.png')

Code

plot_educm2

Plot of Model 2: Marital status

This model has controls, but no weights, no interaction

Code

plot_unmarriedm2 <- plot_predictions(m2, condition = c("score", "unmarried")) +
  labs(title = "Predicted Probability of High Eviction Risk, Marital Status",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Marital Status",
                       breaks=c("1", "0"),
                       labels=c("Unmarried", "Married")
                       ) +
  scale_fill_discrete(name="Marital Status",
                      breaks=c("1", "0"),
                      labels=c("Unmarried", "Married")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('unmarriedm2.png')

Code

plot_unmarriedm2

Plot of Model 2: Presence of Children in Household

Code

plot_childrenm2 <- plot_predictions(m2, condition = c("score", "children")) +
  labs(title = "Predicted Probability of High Eviction Risk, Presence of Children in Household",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Presence of Children in Household",
                       breaks=c("1", "0"),
                       labels=c("Children in Household", "No Children in Household")
                       ) +
  scale_fill_discrete(name="Presence of Children in Household",
                      breaks=c("1", "0"),
                      labels=c("Children in Household", "No Children in Household")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('childrenm2.png')

Code

plot_childrenm2

(I’m not adding model 3s because they aren’t weighted)

Plots for Model 5: With controls and with weight

Plot of Model 5: High Eviction Risk by State Score

With controls, with weight

Code

plot_scorem5 <- plot_predictions(m5, condition = c("score")) +
  labs(title = "Predicted Probability of High Likelihood of Eviction by State Score",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Score",
                       breaks=c(),
                       labels=c("Score")
                       ) +
  scale_fill_discrete(name="Score",
                       breaks=c(),
                       labels=c("Score")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Likelihood of Eviction") + # ylim(0, 0)
  png('plot_scorem5.png')

Code

plot_scorem5

Code

ggsave("plot_scorem5.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 5: Race/Ethnicity

Code

plot_m5race_eth <- plot_predictions(m5, condition = c("score", "race_eth")) +
  labs(title = "Predicted Probability of High Likelihood of Eviction, by Race/Ethnicity",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Race/Ethnicity",
                       breaks=c("nh_black", "other", "nh_white", "hispanic", "nh_asian"),
                       labels=c("Black", "Other", "White",  "Hispanic", "Asian")
                       ) +
  scale_fill_discrete(name="Race/Ethnicity",
                       breaks=c("nh_black", "other", "nh_white", "hispanic", "nh_asian"),
                       labels=c("Black", "Other", "White",  "Hispanic", "Asian")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Likelihood of Eviction") + # ylim(0, 0)
  png('m5race_eth.png')

Code

plot_m5race_eth

Code

ggsave("plot_m5race_eth.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 5: Income

Code

plot_m5inclvl <- plot_predictions(m5, condition = c("score", "inclvl")) +
  labs(title = "Predicted Probability of High Eviction Risk, Income Level",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Income Level",
                       breaks=c("lessthan25", "25to49", "50andover"),
                       labels=c("Less than $25K", "$25K to $49K", "$50K and over")
                       ) +
  scale_fill_discrete(name="Income Level",
                      breaks=c("lessthan25", "25to49", "50andover"),
                      labels=c("Less than $25K", "$25K to $49K", "$50K and over")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('plot_m5inclvl.png')

Code

plot_m5inclvl

Code

ggsave("plot_m5inclvl.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 5: Education

Code

plot_m5educ <- plot_predictions(m5, condition = c("score", "educ")) +
  labs(title = "Predicted Probability of High Eviction Risk, Education Level",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Education Level",
                       breaks=c("HSorless", "somecollege", "BAorhigher"),
                       labels=c("Highschool or less", "Some college", "BA or Higher")
                       ) +
  scale_fill_discrete(name="Education Level",
                      breaks=c("HSorless", "somecollege", "BAorhigher"),
                      labels=c("Highschool or less", "Some college", "BA or Higher")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('plot_m5educ')

Code

plot_m5educ

Code

ggsave("plot_m5educ.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 5: Marital Status

Code

plot_m5unmarried <- plot_predictions(m5, condition = c("score", "unmarried")) +
  labs(title = "Predicted Probability of High Likelihood of Eviction, Marital Status",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Marital Status",
                       breaks=c("1", "0"),
                       labels=c("Unmarried", "Married")
                       ) +
  scale_fill_discrete(name="Marital Status",
                      breaks=c("1", "0"),
                      labels=c("Unmarried", "Married")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Likelihood of Eviction") + # ylim(0, 0)
  png('m5unmarried.png')

Code

plot_m5unmarried

Code

ggsave("plot_m5unmarried.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 5: Presence of Children in Household

Code

plot_m5children <- plot_predictions(m5, condition = c("score", "children")) +
  labs(title = "Predicted Probability of High Likelihood of Eviction, Children in Household",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Presence of\nChildren in Household",
                       breaks=c("1", "0"),
                       labels=c("Children in Household", "No Children in Household")
                       ) +
  scale_fill_discrete(name="Presence of\nChildren in Household",
                      breaks=c("1", "0"),
                      labels=c("Children in Household", "No Children in Household")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Likelihood of Eviction") + # ylim(0, 0)
  png('m5children.png')

Code

plot_m5children

Code

ggsave("plot_m5children.png", path = "./figures/")

Saving 7 x 5 in image

##Plots for Model 6: full interactions with score

Model 6s are score * all covariates interactions.

Plot of Model 6: Race/Ethnicity

Code

plot_m6race_eth <- plot_predictions(m6, condition = c("score", "race_eth")) +
  labs(title = "Predicted Probability of High Likelihood of Eviction, by Race/Ethnicity",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Race/Ethnicity",
                       breaks=c("nh_white", "nh_black", "hispanic", "nh_asian", "other"),
                       labels=c("White", "Black", "Hispanic", "Asian", "Other")
                       ) +
  scale_fill_discrete(name="Race/Ethnicity",
                       breaks=c("nh_white", "nh_black", "hispanic", "nh_asian", "other"),
                       labels=c("White", "Black", "Hispanic", "Asian", "Other")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Likelihood of Eviction") + # ylim(0, 0)
  png('m6race_eth.png')

Code

plot_m6race_eth

Code

ggsave("plot_m6race_eth.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 6: Income

Code

plot_m6inclvl <- plot_predictions(m6, condition = c("score", "inclvl")) +
  labs(title = "Predicted Probability of High Eviction Risk, Income Level",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Income Level",
                       breaks=c("lessthan25", "25to49", "50andover"),
                       labels=c("Less than $25K", "$25K to $49K", "$50K and over")
                       ) +
  scale_fill_discrete(name="Income Level",
                      breaks=c("lessthan25", "25to49", "50andover"),
                      labels=c("Less than $25K", "$25K to $49K", "$50K and over")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('plot_m6inclvl.png')

Code

plot_m6inclvl

Code

ggsave("plot_m6inclvl.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 6: Education

Code

plot_m6educ <- plot_predictions(m6, condition = c("score", "educ")) +
  labs(title = "Predicted Probability of High Eviction Risk, Education Level",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Education Level",
                       breaks=c("HSorless", "somecollege", "BAorhigher"),
                       labels=c("Highschool or less", "Some college", "BA or Higher")
                       ) +
  scale_fill_discrete(name="Education Level",
                      breaks=c("HSorless", "somecollege", "BAorhigher"),
                      labels=c("Highschool or less", "Some college", "BA or Higher")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('plot_m6educ')

Code

plot_m6educ

Code

ggsave("plot_m6educ.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 6: Marital Status

Code

plot_m6unmarried <- plot_predictions(m6, condition = c("score", "unmarried")) +
  labs(title = "Predicted Probability of High Eviction Risk, Marital Status",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Marital Status",
                       breaks=c("1", "0"),
                       labels=c("Unmarried", "Married")
                       ) +
  scale_fill_discrete(name="Marital Status",
                      breaks=c("1", "0"),
                      labels=c("Unmarried", "Married")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('m6unmarried.png')

Code

plot_m6unmarried

Code

ggsave("plot_m5unmarried.png", path = "./figures/")

Saving 7 x 5 in image

Plot of Model 6: Presence of Children in Household

Code

plot_m6children <- plot_predictions(m6, condition = c("score", "children")) +
  labs(title = "Predicted Probability of High Eviction Risk, Presence of Children in Household",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
  scale_color_discrete(name="Presence of Children in Household",
                       breaks=c("1", "0"),
                       labels=c("Children in Household", "No Children in Household")
                       ) +
  scale_fill_discrete(name="Presence of Children in Household",
                      breaks=c("1", "0"),
                      labels=c("Children in Household", "No Children in Household")
                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('m6children.png')

Code

plot_m5children

Code

ggsave("plot_m6children.png", path = "./figures/")

Saving 7 x 5 in image

Formatted tables for PAA

Model 1 and 2

Code

t1 <- tbl_regression(m1, exponentiate = TRUE,
                     label = list(
                       score ~ "State Score"
                       )
                     ) %>% 
  bold_labels()

t1

Characteristic	OR	95% CI	p-value
State Score	0.91	0.87, 0.95	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

PAA: Can I modify how the categorical variables display?

Code

t2 <- tbl_regression(m2,
                     exponentiate = TRUE,
                     include = c(score,
                                 age,
                                 female,
                                 race_eth,
                                 unmarried,
                                 children,
                                 inclvl,
                                 educ),
                     label = list(
                       score ~ "State Score",
                       age ~ "Age",
                       female ~ "Female",
                       race_eth ~ "Race/ethnicity",
                       unmarried ~ "Unmarried",
                       inclvl ~ "Income level",
                       children ~ "Children in Household",
                       educ ~ "Education"
                       )
                     )%>% 
  bold_labels()

t2

Characteristic	OR	95% CI	p-value
State Score	0.93	0.89, 0.97	0.001
Age	1.00	0.99, 1.00	0.001
Female	0.95	0.88, 1.01	0.12
Race/ethnicity
nh_white	—	—
nh_black	1.27	1.17, 1.38	<0.001
hispanic	0.87	0.79, 0.95	0.003
nh_asian	0.32	0.26, 0.38	<0.001
other	1.08	0.95, 1.23	0.3
Unmarried	1.25	1.16, 1.35	<0.001
Children in Household	1.19	1.11, 1.28	<0.001
Income level
25to49	—	—
lessthan25	1.39	1.30, 1.50	<0.001
50andover	0.67	0.61, 0.73	<0.001
Education
somecollege	—	—
HSorless	0.96	0.89, 1.05	0.4
BAorhigher	0.76	0.70, 0.82	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Code

model12 <- tbl_merge(
  tbls = list(t1, t2),
  tab_spanner = c("**Model 1**", "**Model 2**")) #%>% 
  #as_gt() %>% # convert to gt table
  
model12

Characteristic	Model 1			Model 2
	OR	95% CI	p-value	OR	95% CI	p-value
State Score	0.91	0.87, 0.95	<0.001	0.93	0.89, 0.97	0.001
Age				1.00	0.99, 1.00	0.001
Female				0.95	0.88, 1.01	0.12
Race/ethnicity
nh_white				—	—
nh_black				1.27	1.17, 1.38	<0.001
hispanic				0.87	0.79, 0.95	0.003
nh_asian				0.32	0.26, 0.38	<0.001
other				1.08	0.95, 1.23	0.3
Unmarried				1.25	1.16, 1.35	<0.001
Children in Household				1.19	1.11, 1.28	<0.001
Income level
25to49				—	—
lessthan25				1.39	1.30, 1.50	<0.001
50andover				0.67	0.61, 0.73	<0.001
Education
somecollege				—	—
HSorless				0.96	0.89, 1.05	0.4
BAorhigher				0.76	0.70, 0.82	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Code

gtsave(as_gt(model12), filename = "model12.png", path = "./figures/") # save table as image

file:////var/folders/n0/x0mfrjpd1p75y8wj5zn84ztc0000gn/T//Rtmp59cpxX/file25b728d460d5.html screenshot completed

Model 4 and 5

Code

t4 <- tbl_regression(m4, exponentiate = TRUE,
                     label = list(
                       score ~ "State Score"
                       )
                     ) %>% 
  bold_labels()

t4

Characteristic	OR	95% CI	p-value
State Score	0.93	0.89, 0.98	0.005
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Code

t5 <- tbl_regression(m5,
                     exponentiate = TRUE,
                     include = c(score,
                                 age,
                                 female,
                                 race_eth,
                                 unmarried,
                                 children,
                                 inclvl,
                                 educ),
                     label = list(
                       score ~ "State Score",
                       age ~ "Age",
                       female ~ "Female",
                       race_eth ~ "Race/ethnicity",
                       unmarried ~ "Unmarried",
                       inclvl ~ "Income level",
                       children ~ "Children in Household",
                       educ ~ "Education"
                       )
                     )%>% 
  bold_labels()

t5

Characteristic	OR	95% CI	p-value
State Score	0.93	0.89, 0.97	0.002
Age	1.00	0.99, 1.00	0.008
Female	0.95	0.88, 1.02	0.2
Race/ethnicity
nh_white	—	—
nh_black	1.26	1.16, 1.37	<0.001
hispanic	0.86	0.77, 0.96	0.006
nh_asian	0.35	0.28, 0.44	<0.001
other	1.07	0.95, 1.21	0.3
Unmarried	1.24	1.15, 1.34	<0.001
Children in Household	1.20	1.12, 1.29	<0.001
Income level
25to49	—	—
lessthan25	1.41	1.31, 1.52	<0.001
50andover	0.69	0.62, 0.77	<0.001
Education
somecollege	—	—
HSorless	0.98	0.90, 1.06	0.6
BAorhigher	0.78	0.72, 0.84	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Code

model45 <- tbl_merge(
  tbls = list(t4, t5),
  tab_spanner = c("**Model 4**", "**Model 5**")) #%>% 
  #as_gt() %>% # convert to gt table
  
model45

Characteristic	Model 4			Model 5
	OR	95% CI	p-value	OR	95% CI	p-value
State Score	0.93	0.89, 0.98	0.005	0.93	0.89, 0.97	0.002
Age				1.00	0.99, 1.00	0.008
Female				0.95	0.88, 1.02	0.2
Race/ethnicity
nh_white				—	—
nh_black				1.26	1.16, 1.37	<0.001
hispanic				0.86	0.77, 0.96	0.006
nh_asian				0.35	0.28, 0.44	<0.001
other				1.07	0.95, 1.21	0.3
Unmarried				1.24	1.15, 1.34	<0.001
Children in Household				1.20	1.12, 1.29	<0.001
Income level
25to49				—	—
lessthan25				1.41	1.31, 1.52	<0.001
50andover				0.69	0.62, 0.77	<0.001
Education
somecollege				—	—
HSorless				0.98	0.90, 1.06	0.6
BAorhigher				0.78	0.72, 0.84	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Code

gtsave(as_gt(model45), filename = "model45.png", path = "./figures/") # save table as image

file:////var/folders/n0/x0mfrjpd1p75y8wj5zn84ztc0000gn/T//Rtmp59cpxX/file25b73b9735a3.html screenshot completed

Model 2 (unweighted) and 5 (weighted)

Code

model25 <- tbl_merge(
  tbls = list(t2, t5),
  tab_spanner = c("**Model 2 (unweighted)**", "**Model 5 (weighted)**")
  ) # %>%
  # as_gt() # %>% # convert to gt table
  # gt::gtsave(filename = "model45.png") # save table as image

model25

Characteristic	Model 2 (unweighted)			Model 5 (weighted)
	OR	95% CI	p-value	OR	95% CI	p-value
State Score	0.93	0.89, 0.97	0.001	0.93	0.89, 0.97	0.002
Age	1.00	0.99, 1.00	0.001	1.00	0.99, 1.00	0.008
Female	0.95	0.88, 1.01	0.12	0.95	0.88, 1.02	0.2
Race/ethnicity
nh_white	—	—		—	—
nh_black	1.27	1.17, 1.38	<0.001	1.26	1.16, 1.37	<0.001
hispanic	0.87	0.79, 0.95	0.003	0.86	0.77, 0.96	0.006
nh_asian	0.32	0.26, 0.38	<0.001	0.35	0.28, 0.44	<0.001
other	1.08	0.95, 1.23	0.3	1.07	0.95, 1.21	0.3
Unmarried	1.25	1.16, 1.35	<0.001	1.24	1.15, 1.34	<0.001
Children in Household	1.19	1.11, 1.28	<0.001	1.20	1.12, 1.29	<0.001
Income level
25to49	—	—		—	—
lessthan25	1.39	1.30, 1.50	<0.001	1.41	1.31, 1.52	<0.001
50andover	0.67	0.61, 0.73	<0.001	0.69	0.62, 0.77	<0.001
Education
somecollege	—	—		—	—
HSorless	0.96	0.89, 1.05	0.4	0.98	0.90, 1.06	0.6
BAorhigher	0.76	0.70, 0.82	<0.001	0.78	0.72, 0.84	<0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Code

gtsave(as_gt(model25), filename = "model25.png", path = "./figures/") # save table as image

file:////var/folders/n0/x0mfrjpd1p75y8wj5zn84ztc0000gn/T//Rtmp59cpxX/file25b73fc1e60c.html screenshot completed

Prediction() code

Code

pre <- predictions(m2, newdata = datagrid(score = unique, highrisk = unique))

pre

 score Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
  0.83    0.494     0.0136 36.4   <0.001 963.2 0.468  0.521
  0.83    0.494     0.0136 36.4   <0.001 963.2 0.468  0.521
  1.23    0.487     0.0128 38.0   <0.001   Inf 0.462  0.512
  1.23    0.487     0.0128 38.0   <0.001   Inf 0.462  0.512
  0.23    0.505     0.0152 33.2   <0.001 800.3 0.475  0.535
--- 72 rows omitted. See ?print.marginaleffects --- 
  3.75    0.443     0.0159 27.8   <0.001 562.6 0.412  0.474
  0.63    0.498     0.0141 35.4   <0.001 910.9 0.470  0.525
  0.63    0.498     0.0141 35.4   <0.001 910.9 0.470  0.525
  3.00    0.456     0.0138 33.1   <0.001 797.5 0.429  0.483
  3.00    0.456     0.0138 33.1   <0.001 797.5 0.429  0.483
Type:  response 

Some exploratory models

Data set with children and marital status groups

Dummy variable for unmarried households with children

Make dummy variable

Code

latehh <- latehh %>%
  mutate(unwch = case_when(unmarried == 1 & children == 1 ~ 1,
                              TRUE ~ 0
                              )
        )

tabyl(latehh$unwch)

 latehh$unwch     n   percent
            0 11658 0.6997599
            1  5002 0.3002401

Model 7, variable for unmarried with children, no controls

Code

m7 <- glmer(highrisk ~ score +
               unwch +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Code

summary(m7)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + unwch + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20943.3   20974.2  -10467.7   20935.3     16656 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.693 -0.858 -0.562  1.034  1.818 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.04059  0.2015  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
            Estimate Std. Error z value             Pr(>|z|)    
(Intercept) -0.24747    0.05712  -4.333            0.0000147 ***
score       -0.06833    0.02419  -2.824              0.00474 ** 
unwch        0.40306    0.03468  11.624 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
      (Intr) score 
score -0.782       
unwch -0.209 -0.004

Model 8, Variable for unmarried with children, with controls

Code

m8 <- glmer(highrisk ~ score +
               unwch +
              age +
              female +
              race_eth +
              inclvl +
              educ +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00453899 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m8)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + unwch + age + female + race_eth + inclvl +  
    educ + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20427.3   20535.4  -10199.7   20399.3     16646 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9947 -0.8455 -0.4304  0.9948  2.8230 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03479  0.1865  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                  Estimate Std. Error z value             Pr(>|z|)    
(Intercept)      -0.014398   0.091908  -0.157              0.87551    
score            -0.071581   0.023046  -3.106              0.00190 ** 
unwch             0.265478   0.037126   7.151    0.000000000000864 ***
age              -0.003770   0.001353  -2.785              0.00535 ** 
female           -0.057200   0.037701  -1.517              0.12922    
race_ethnh_black  0.227960   0.042706   5.338    0.000000094062715 ***
race_ethhispanic -0.161077   0.055715  -2.891              0.00384 ** 
race_ethnh_asian -1.069457   0.118202  -9.048 < 0.0000000000000002 ***
race_ethother     0.072039   0.061367   1.174              0.24043    
inclvllessthan25  0.352946   0.036961   9.549 < 0.0000000000000002 ***
inclvl50andover  -0.378531   0.051365  -7.369    0.000000000000171 ***
educHSorless     -0.027914   0.041888  -0.666              0.50517    
educBAorhigher   -0.245683   0.040545  -6.060    0.000000001365147 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 13 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00453899 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Model 9, Variable for unmarried with children, with controls, with full interaction

Code

m9 <- glmer(highrisk ~ score * (
               unwch +
              age +
              female +
              race_eth +
              inclvl +
              educ
              ) +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0653454 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m9)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score * (unwch + age + female + race_eth + inclvl +  
    educ) + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20436.5   20629.5  -10193.2   20386.5     16635 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9661 -0.8481 -0.4313  0.9905  2.8314 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03388  0.1841  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                         Estimate Std. Error z value     Pr(>|z|)    
(Intercept)             0.0373518  0.1394439   0.268     0.788805    
score                  -0.0987694  0.0615340  -1.605     0.108468    
unwch                   0.2038792  0.0638272   3.194     0.001402 ** 
age                    -0.0034762  0.0023301  -1.492     0.135739    
female                 -0.0527400  0.0647190  -0.815     0.415125    
race_ethnh_black        0.3290438  0.0728961   4.514 0.0000063655 ***
race_ethhispanic       -0.2660531  0.0966788  -2.752     0.005925 ** 
race_ethnh_asian       -1.0648351  0.1989006  -5.354 0.0000000862 ***
race_ethother           0.0540732  0.1071038   0.505     0.613652    
inclvllessthan25        0.2768642  0.0636996   4.346 0.0000138387 ***
inclvl50andover        -0.4605855  0.0875580  -5.260 0.0000001438 ***
educHSorless           -0.0824238  0.0725595  -1.136     0.255979    
educBAorhigher         -0.2344173  0.0696925  -3.364     0.000769 ***
score:unwch             0.0324294  0.0280416   1.156     0.247487    
score:age              -0.0001437  0.0010213  -0.141     0.888077    
score:female           -0.0037007  0.0285332  -0.130     0.896806    
score:race_ethnh_black -0.0546645  0.0319464  -1.711     0.087057 .  
score:race_ethhispanic  0.0557205  0.0427389   1.304     0.192321    
score:race_ethnh_asian -0.0030315  0.0906354  -0.033     0.973318    
score:race_ethother     0.0098417  0.0476341   0.207     0.836314    
score:inclvllessthan25  0.0415685  0.0280391   1.483     0.138203    
score:inclvl50andover   0.0455356  0.0389049   1.170     0.241826    
score:educHSorless      0.0295283  0.0317860   0.929     0.352903    
score:educBAorhigher   -0.0069892  0.0308401  -0.227     0.820714    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0653454 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Plot of Model 9

Code

plot_m9 <- plot_predictions(m9, condition = c("score", "unwch")) +
  labs(title = "Predicted Probability of High Eviction Risk, Unmarried with Children",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
# scale_color_discrete(name="Presence of Children in Household",
#                       breaks=c("1", "0"),
#                       labels=c("Children in Household", "No Children in Household")
#                       ) +
#  scale_fill_discrete(name="Presence of Children in Household",
#                      breaks=c("1", "0"),
#                      labels=c("Children in Household", "No Children in Household")
#                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('m9.png')

Code

plot_m9

Code

ggsave("plot_m9.png", path = "./figures/")

Saving 7 x 5 in image

Categorical with 4 variables for marrital status and child status

Make dummy variable

Code

latehh <- latehh %>%
  mutate(marrchild = case_when(married == 0 & children == 1 ~ "UnWCh",
                               married == 1 & children == 1 ~ "MWCh",
                               married == 0 & children == 0 ~ "UnNoCh",
                               married == 1 & children == 0 ~ "MNoCh"
                              )
        )

Code

latehh$marrchild <- factor(latehh$marrchild, levels=c("UnNoCh",  "UnWCh",  "MNoCh", "MWCh"))

tabyl(latehh$marrchild)

 latehh$marrchild    n   percent
           UnNoCh 6684 0.4012005
            UnWCh 5002 0.3002401
            MNoCh 2111 0.1267107
             MWCh 2863 0.1718487

Model 10, categorical marital status by children, no controls

Code

m10 <- glmer(highrisk ~ score +
               marrchild +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Code

summary(m10)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + marrchild + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20907.3   20953.6  -10447.7   20895.3     16654 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-1.696 -0.850 -0.558  1.023  1.935 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.0417   0.2042  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
               Estimate Std. Error z value           Pr(>|z|)    
(Intercept)    -0.14055    0.06012  -2.338             0.0194 *  
score          -0.06829    0.02444  -2.795             0.0052 ** 
marrchildUnWCh  0.29531    0.03868   7.635 0.0000000000000225 ***
marrchildMNoCh -0.29639    0.05709  -5.191 0.0000002087591866 ***
marrchildMWCh  -0.23051    0.04813  -4.790 0.0000016698233000 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) score  mrrUWC mrrMNC
score       -0.751                     
mrrchldUnWC -0.303 -0.003              
mrrchldMNCh -0.203  0.002  0.318       
mrrchldMWCh -0.241 -0.001  0.378  0.256

Model 11, categorical marital status by children with controls

Code

m11 <- glmer(highrisk ~ score +
               marrchild +
               age +
               female +
               race_eth +
               inclvl +
               educ +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00390339 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m11)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + marrchild + age + female + race_eth + inclvl +  
    educ + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20427.1   20550.6  -10197.6   20395.1     16644 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9863 -0.8458 -0.4325  0.9922  2.7394 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03518  0.1876  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                  Estimate Std. Error z value             Pr(>|z|)    
(Intercept)       0.022730   0.094957   0.239              0.81081    
score            -0.071574   0.023134  -3.094              0.00198 ** 
marrchildUnWCh    0.231189   0.041335   5.593     0.00000002230527 ***
marrchildMNoCh   -0.107625   0.059015  -1.824              0.06820 .  
marrchildMWCh    -0.069858   0.050969  -1.371              0.17050    
age              -0.003828   0.001379  -2.776              0.00550 ** 
female           -0.055347   0.037815  -1.464              0.14329    
race_ethnh_black  0.227050   0.042745   5.312     0.00000010857577 ***
race_ethhispanic -0.156088   0.055848  -2.795              0.00519 ** 
race_ethnh_asian -1.056912   0.118422  -8.925 < 0.0000000000000002 ***
race_ethother     0.071535   0.061403   1.165              0.24401    
inclvllessthan25  0.343537   0.037280   9.215 < 0.0000000000000002 ***
inclvl50andover  -0.367588   0.051698  -7.110     0.00000000000116 ***
educHSorless     -0.023913   0.041957  -0.570              0.56872    
educBAorhigher   -0.246505   0.040553  -6.079     0.00000000121244 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 15 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00390339 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Model 12, categorical marital status by children, full interactions

Code

m12 <- glmer(highrisk ~ score *(
               marrchild +
               age +
               female +
               race_eth +
               inclvl +
               educ
               ) +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.105446 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m12)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score * (marrchild + age + female + race_eth + inclvl +  
    educ) + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20439.8   20663.7  -10190.9   20381.8     16631 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9571 -0.8484 -0.4301  0.9882  2.7295 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03436  0.1854  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                         Estimate Std. Error z value    Pr(>|z|)    
(Intercept)             0.0906206  0.1455483   0.623    0.533537    
score                  -0.1076590  0.0639798  -1.683    0.092433 .  
marrchildUnWCh          0.1547507  0.0710004   2.180    0.029289 *  
marrchildMNoCh         -0.1642321  0.1004976  -1.634    0.102219    
marrchildMWCh          -0.0952141  0.0877388  -1.085    0.277834    
age                    -0.0035254  0.0023763  -1.484    0.137925    
female                 -0.0513630  0.0648807  -0.792    0.428563    
race_ethnh_black        0.3269947  0.0729593   4.482 0.000007399 ***
race_ethhispanic       -0.2578701  0.0969247  -2.661    0.007802 ** 
race_ethnh_asian       -1.0514974  0.1990398  -5.283 0.000000127 ***
race_ethother           0.0535155  0.1071793   0.499    0.617562    
inclvllessthan25        0.2629087  0.0642766   4.090 0.000043087 ***
inclvl50andover        -0.4422272  0.0882651  -5.010 0.000000544 ***
educHSorless           -0.0773641  0.0726902  -1.064    0.287193    
educBAorhigher         -0.2363605  0.0697155  -3.390    0.000698 ***
score:marrchildUnWCh    0.0407562  0.0312495   1.304    0.192160    
score:marrchildMNoCh    0.0314675  0.0442885   0.711    0.477387    
score:marrchildMWCh     0.0143602  0.0390393   0.368    0.712992    
score:age              -0.0001501  0.0010396  -0.144    0.885182    
score:female           -0.0034366  0.0286367  -0.120    0.904477    
score:race_ethnh_black -0.0541645  0.0319820  -1.694    0.090343 .  
score:race_ethhispanic  0.0539764  0.0428494   1.260    0.207786    
score:race_ethnh_asian -0.0036068  0.0907044  -0.040    0.968281    
score:race_ethother     0.0098069  0.0476581   0.206    0.836965    
score:inclvllessthan25  0.0439919  0.0282915   1.555    0.119958    
score:inclvl50andover   0.0417301  0.0391788   1.065    0.286822    
score:educHSorless      0.0288530  0.0318610   0.906    0.365152    
score:educBAorhigher   -0.0065008  0.0308480  -0.211    0.833092    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 28 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.105446 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Plot of Model 12

Code

plot_m12 <- plot_predictions(m12, condition = c("score", "marrchild")) +
  labs(title = "Predicted Probability of High Eviction Risk, Unmarried with Children",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
# scale_color_discrete(name="Presence of Children in Household",
#                       breaks=c("1", "0"),
#                       labels=c("Children in Household", "No Children in Household")
#                       ) +
#  scale_fill_discrete(name="Presence of Children in Household",
#                      breaks=c("1", "0"),
#                      labels=c("Children in Household", "No Children in Household")
#                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('m12.png')

Code

plot_m12

Code

ggsave("plot_m12.png", path = "./figures/")

Saving 7 x 5 in image

Categorical with marital x children x race

Make Dummy Variable

Code

tabyl(latehh, race_eth, unwch)

 race_eth    0    1
 nh_white 5620 1979
 nh_black 2324 1530
 hispanic 1978  982
 nh_asian  979  106
    other  757  405

Code

latehh <- latehh %>%
  mutate(rsp = case_when(race_eth == "nh_white" & unwch == 1 ~ "wsp",
                         race_eth == "nh_white" & unwch == 0 ~ "wnsp",
                         race_eth == "nh_black" & unwch == 1 ~ "bsp",
                         race_eth == "nh_black" & unwch == 0 ~ "bnsp",
                         race_eth == "hispanic" & unwch == 1 ~ "hsp",
                         race_eth == "hispanic" & unwch == 0 ~ "hnsp",
                         race_eth == "nh_asian" & unwch == 1 ~ "asp",
                         race_eth == "nh_asian" & unwch == 0 ~ "ansp",
                         race_eth == "other" & unwch == 1 ~ "osp",
                         race_eth == "other" & unwch == 0 ~ "onsp"
                         )
        )

latehh$rsp <- factor(latehh$rsp, levels=c("wnsp", "wsp", "bnsp", "bsp", "hnsp", "hsp", "ansp", "asp", "onsp", "osp"))

tabyl(latehh$rsp)

 latehh$rsp    n     percent
       wnsp 5620 0.337334934
        wsp 1979 0.118787515
       bnsp 2324 0.139495798
        bsp 1530 0.091836735
       hnsp 1978 0.118727491
        hsp  982 0.058943577
       ansp  979 0.058763505
        asp  106 0.006362545
       onsp  757 0.045438175
        osp  405 0.024309724

Model 13, race by unmarried parent, no controls

Code

m13 <- glmer(highrisk ~ score +
               rsp +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Code

summary(m13)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + rsp + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20745.7   20838.4  -10360.9   20721.7     16648 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6974 -0.8690 -0.5022  1.0204  2.6360 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03539  0.1881  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
            Estimate Std. Error z value             Pr(>|z|)    
(Intercept) -0.23738    0.05738  -4.137        0.00003522332 ***
score       -0.06987    0.02310  -3.025             0.002488 ** 
rspwsp       0.30369    0.05057   6.005        0.00000000191 ***
rspbnsp      0.20120    0.05252   3.831             0.000128 ***
rspbsp       0.63219    0.05841  10.823 < 0.0000000000000002 ***
rsphnsp     -0.09994    0.06780  -1.474             0.140431    
rsphsp       0.19125    0.08315   2.300             0.021450 *  
rspansp     -1.37208    0.13048 -10.516 < 0.0000000000000002 ***
rspasp      -0.23466    0.26731  -0.878             0.380030    
rsponsp      0.10557    0.07530   1.402             0.160907    
rsposp       0.42354    0.09478   4.469        0.00000787438 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        (Intr) score  rspwsp rspbns rspbsp rsphns rsphsp rspansp rspasp rsponsp
score   -0.740                                                                 
rspwsp  -0.259 -0.004                                                          
rspbnsp -0.242 -0.016  0.279                                                   
rspbsp  -0.223 -0.010  0.252  0.298                                            
rsphnsp -0.186 -0.001  0.210  0.220  0.193                                     
rsphsp  -0.147 -0.008  0.172  0.179  0.157  0.153                              
rspansp -0.098  0.000  0.109  0.117  0.102  0.095  0.076                       
rspasp  -0.044 -0.004  0.053  0.055  0.048  0.045  0.037  0.026                
rsponsp -0.167 -0.008  0.193  0.191  0.170  0.155  0.125  0.084   0.042        
rsposp  -0.142 -0.001  0.156  0.153  0.138  0.121  0.098  0.066   0.032  0.121 

Model 14, race by unmarried parent with controls

Code

m14 <- glmer(highrisk ~ score +
               rsp +
               age +
               female +
               race_eth +
               inclvl +
               educ +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m14)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + rsp + age + female + race_eth + inclvl + educ +  
    (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20427.7   20566.6  -10195.8   20391.7     16642 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9746 -0.8442 -0.4321  0.9927  2.9912 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03461  0.186   
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                  Estimate Std. Error z value             Pr(>|z|)    
(Intercept)      -0.005746   0.092486  -0.062             0.950458    
score            -0.071487   0.023007  -3.107             0.001889 ** 
rspwsp            0.228100   0.052186   4.371    0.000012370887321 ***
rspbnsp           0.186413   0.053208   3.503             0.000459 ***
rspbsp            0.523588   0.060480   8.657 < 0.0000000000000002 ***
rsphnsp          -0.146840   0.068901  -2.131             0.033075 *  
rsphsp            0.053099   0.085237   0.623             0.533318    
rspansp          -1.198311   0.132567  -9.039 < 0.0000000000000002 ***
rspasp           -0.224101   0.270003  -0.830             0.406543    
rsponsp           0.073219   0.076203   0.961             0.336630    
rsposp            0.307201   0.096429   3.186             0.001444 ** 
age              -0.003727   0.001356  -2.749             0.005986 ** 
female           -0.057820   0.037708  -1.533             0.125182    
inclvllessthan25  0.352344   0.036977   9.529 < 0.0000000000000002 ***
inclvl50andover  -0.376925   0.051385  -7.335    0.000000000000221 ***
educHSorless     -0.026136   0.041900  -0.624             0.532774    
educBAorhigher   -0.244400   0.040549  -6.027    0.000000001668089 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 17 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
optimizer (bobyqa) convergence code: 0 (OK)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Model 15, race by unmarried parent, full interactions

Code

m15 <- glmer(highrisk ~ score *(
               rsp +
               age +
               female +
               race_eth +
               inclvl +
               educ
               ) +
              (1 | abbv), 
            data = latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

fixed-effect model matrix is rank deficient so dropping 8 columns / coefficients

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
length(par)^2 is not recommended.

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0805244 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

Code

summary(m15)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score * (rsp + age + female + race_eth + inclvl +  
    educ) + (1 | abbv)
   Data: latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
  20439.7   20694.5  -10186.8   20373.7     16627 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9799 -0.8488 -0.4339  0.9899  2.8986 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03383  0.1839  
Number of obs: 16660, groups:  abbv, 51

Fixed effects:
                         Estimate Std. Error z value    Pr(>|z|)    
(Intercept)             0.0641207  0.1407809   0.455    0.648775    
score                  -0.1093682  0.0622251  -1.758    0.078812 .  
rspwsp                  0.0927495  0.0911835   1.017    0.309070    
rspbnsp                 0.2362055  0.0900606   2.623    0.008723 ** 
rspbsp                  0.5833702  0.1028185   5.674 0.000000014 ***
rsphnsp                -0.2802651  0.1191463  -2.352    0.018659 *  
rsphsp                 -0.1247476  0.1488925  -0.838    0.402122    
rspansp                -1.0926048  0.2176661  -5.020 0.000000518 ***
rspasp                 -0.8066620  0.5062529  -1.593    0.111071    
rsponsp                -0.0397922  0.1347117  -0.295    0.767698    
rsposp                  0.3243422  0.1662545   1.951    0.051072 .  
age                    -0.0033886  0.0023346  -1.451    0.146652    
female                 -0.0503728  0.0647237  -0.778    0.436407    
inclvllessthan25        0.2741467  0.0637510   4.300 0.000017059 ***
inclvl50andover        -0.4586747  0.0875706  -5.238 0.000000163 ***
educHSorless           -0.0820173  0.0726155  -1.129    0.258698    
educBAorhigher         -0.2337224  0.0697164  -3.352    0.000801 ***
score:rspwsp            0.0731741  0.0402135   1.820    0.068814 .  
score:rspbnsp          -0.0267827  0.0393473  -0.681    0.496078    
score:rspbsp           -0.0329263  0.0451987  -0.728    0.466320    
score:rsphnsp           0.0729360  0.0535424   1.362    0.173132    
score:rsphsp            0.0929809  0.0645005   1.442    0.149428    
score:rspansp          -0.0632619  0.1038400  -0.609    0.542375    
score:rspasp            0.2971325  0.2117226   1.403    0.160496    
score:rsponsp           0.0613266  0.0595241   1.030    0.302878    
score:rsposp           -0.0110199  0.0741890  -0.149    0.881918    
score:age              -0.0001667  0.0010244  -0.163    0.870742    
score:female           -0.0059548  0.0285555  -0.209    0.834813    
score:inclvllessthan25  0.0433352  0.0280906   1.543    0.122906    
score:inclvl50andover   0.0460008  0.0389310   1.182    0.237365    
score:educHSorless      0.0306657  0.0318038   0.964    0.334938    
score:educBAorhigher   -0.0059675  0.0308566  -0.193    0.846651    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

fit warnings:
fixed-effect model matrix is rank deficient so dropping 8 columns / coefficients
optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0805244 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?
Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?
maxfun < 10 * length(par)^2 is not recommended.

Plot of Model 15

Code

plot_m15 <- plot_predictions(m15, condition = c("score", "rsp")) +
  labs(title = "Predicted Probability of High Eviction Risk",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
# scale_color_discrete(name="Presence of Children in Household",
#                       breaks=c("1", "0"),
#                       labels=c("Children in Household", "No Children in Household")
#                       ) +
#  scale_fill_discrete(name="Presence of Children in Household",
#                      breaks=c("1", "0"),
#                      labels=c("Children in Household", "No Children in Household")
#                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('m15.png')

Code

plot_m15

Code

ggsave("plot_m15.png", path = "./figures/")

Saving 7 x 5 in image

Sample of only unmarried with children, then test by race

Code

tabyl(latehh, unmarried, children)

 unmarried    0    1
         0 2111 2863
         1 6684 5002

Filter for unmarried with children late households

Code

umwc_latehh <- latehh %>% 
  filter(unmarried == 1 & children == 1)

nrow(umwc_latehh)

[1] 5002

Write out csv for unmarried with children late households

Code

# Export to csv files

fwrite(umwc_latehh, "./data/umwc_latehh.csv")

Code

tabyl(umwc_latehh$race_eth)

 umwc_latehh$race_eth    n    percent
             nh_white 1979 0.39564174
             nh_black 1530 0.30587765
             hispanic  982 0.19632147
             nh_asian  106 0.02119152
                other  405 0.08096761

Model 16, race_eth, no controls

Code

m16 <- glmer(highrisk ~ score +
               race_eth +
              (1 | abbv), 
            data = umwc_latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Code

summary(m16)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + race_eth + (1 | abbv)
   Data: umwc_latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
   7118.4    7164.0   -3552.2    7104.4      4995 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6519 -1.0137  0.6641  0.9754  1.5801 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03799  0.1949  
Number of obs: 5002, groups:  abbv, 51

Fixed effects:
                 Estimate Std. Error z value    Pr(>|z|)    
(Intercept)       0.02005    0.07534   0.266      0.7901    
score            -0.05121    0.02959  -1.731      0.0835 .  
race_ethnh_black  0.34981    0.06829   5.123 0.000000302 ***
race_ethhispanic -0.09656    0.09070  -1.065      0.2870    
race_ethnh_asian -0.53061    0.27017  -1.964      0.0495 *  
race_ethother     0.13796    0.10099   1.366      0.1719    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) score  rc_thnh_b rc_thh rc_thnh_s
score       -0.730                                  
rc_thnh_blc -0.357 -0.004                           
rac_thhspnc -0.257 -0.009  0.305                    
rac_thnh_sn -0.084 -0.005  0.102     0.087          
race_eththr -0.249  0.004  0.270     0.209  0.074   

Model 17, race_eth with controls

Code

m17 <- glmer(highrisk ~ score +
               race_eth +
               age +
               female +
               inclvl +
               educ +
              (1 | abbv), 
            data = umwc_latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Code

summary(m17)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score + race_eth + age + female + inclvl + educ +  
    (1 | abbv)
   Data: umwc_latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
   7043.0    7127.7   -3508.5    7017.0      4989 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0736 -0.9739  0.6418  0.9592  1.7104 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03575  0.1891  
Number of obs: 5002, groups:  abbv, 51

Fixed effects:
                  Estimate Std. Error z value        Pr(>|z|)    
(Intercept)       0.525982   0.161205   3.263        0.001103 ** 
score            -0.054359   0.029328  -1.853        0.063813 .  
race_ethnh_black  0.327623   0.068950   4.752 0.0000020182628 ***
race_ethhispanic -0.146581   0.092239  -1.589        0.112026    
race_ethnh_asian -0.503484   0.273345  -1.842        0.065484 .  
race_ethother     0.098867   0.101873   0.970        0.331800    
age              -0.009435   0.002745  -3.437        0.000587 ***
female           -0.242448   0.079239  -3.060        0.002215 ** 
inclvllessthan25  0.403561   0.061745   6.536 0.0000000000632 ***
inclvl50andover  -0.154762   0.100966  -1.533        0.125320    
educHSorless     -0.109409   0.068433  -1.599        0.109868    
educBAorhigher   -0.178576   0.071975  -2.481        0.013098 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) score  rc_thnh_b rc_thh rc_thnh_s rc_tht age    female
score       -0.338                                                       
rc_thnh_blc -0.195 -0.003                                                
rac_thhspnc -0.162 -0.007  0.306                                         
rac_thnh_sn -0.053 -0.007  0.099     0.083                               
race_eththr -0.113  0.005  0.270     0.207  0.074                        
age         -0.722  0.007  0.054     0.055  0.010     0.005              
female      -0.481 -0.006 -0.034     0.037  0.038    -0.001  0.102       
inclvllss25 -0.142 -0.010  0.005    -0.005 -0.012    -0.046 -0.043 -0.062
inclvl50ndv -0.133 -0.002  0.037     0.027 -0.013    -0.002 -0.034  0.079
educHSorlss -0.189 -0.010 -0.007    -0.093  0.010     0.012 -0.002  0.054
educBArhghr -0.093  0.009 -0.007    -0.002 -0.046     0.000 -0.133 -0.007
            incl25 incl50 edcHSr
score                           
rc_thnh_blc                     
rac_thhspnc                     
rac_thnh_sn                     
race_eththr                     
age                             
female                          
inclvllss25                     
inclvl50ndv  0.329              
educHSorlss -0.127  0.021       
educBArhghr  0.012 -0.107  0.434
optimizer (bobyqa) convergence code: 0 (OK)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Model 18, race_eth with interactions

Code

m18 <- glmer(highrisk ~ score*(
               race_eth +
               age +
               female +
               inclvl +
               educ
               ) +
              (1 | abbv), 
            data = umwc_latehh, 
            family = binomial,
            weights = weights,
            control = glmerControl(optimizer = "bobyqa"), 
            nAGQ = 10
            )

Warning in eval(family$initialize, rho): non-integer #successes in a binomial
glm!

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0844054 (tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

Code

summary(m18)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
 Family: binomial  ( logit )
Formula: highrisk ~ score * (race_eth + age + female + inclvl + educ) +  
    (1 | abbv)
   Data: umwc_latehh
Weights: weights
Control: glmerControl(optimizer = "bobyqa")

      AIC       BIC    logLik -2*log(L)  df.resid 
   7045.7    7195.6   -3499.9    6999.7      4979 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0708 -0.9591  0.6367  0.9521  1.8005 

Random effects:
 Groups Name        Variance Std.Dev.
 abbv   (Intercept) 0.03064  0.175   
Number of obs: 5002, groups:  abbv, 51

Fixed effects:
                         Estimate Std. Error z value  Pr(>|z|)    
(Intercept)             0.7432485  0.2605126   2.853   0.00433 ** 
score                  -0.1736950  0.1143315  -1.519   0.12871    
race_ethnh_black        0.5231177  0.1194566   4.379 0.0000119 ***
race_ethhispanic       -0.1931647  0.1619406  -1.193   0.23294    
race_ethnh_asian       -0.9902204  0.5103104  -1.940   0.05233 .  
race_ethother           0.2501404  0.1756727   1.424   0.15448    
age                    -0.0108933  0.0047515  -2.293   0.02187 *  
female                 -0.2859974  0.1357232  -2.107   0.03510 *  
inclvllessthan25        0.1642726  0.1070014   1.535   0.12473    
inclvl50andover        -0.2519375  0.1735705  -1.451   0.14664    
educHSorless           -0.1737181  0.1198154  -1.450   0.14709    
educBAorhigher         -0.2966826  0.1241894  -2.389   0.01690 *  
score:race_ethnh_black -0.1094558  0.0526644  -2.078   0.03768 *  
score:race_ethhispanic  0.0212418  0.0699364   0.304   0.76133    
score:race_ethnh_asian  0.2380784  0.2130540   1.117   0.26380    
score:race_ethother    -0.0894529  0.0785540  -1.139   0.25481    
score:age               0.0008688  0.0020614   0.421   0.67341    
score:female            0.0224642  0.0592969   0.379   0.70481    
score:inclvllessthan25  0.1303695  0.0468075   2.785   0.00535 ** 
score:inclvl50andover   0.0563409  0.0767188   0.734   0.46272    
score:educHSorless      0.0361946  0.0520782   0.695   0.48705    
score:educBAorhigher    0.0667381  0.0546452   1.221   0.22197    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation matrix not shown by default, as p = 22 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0844054 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?
Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

Plot of Model 18

Code

plot_m18 <- plot_predictions(m18, condition = c("score", "race_eth")) +
  labs(title = "Predicted Probability of High Eviction Risk",
       subtitle = "Renters Who Are Not Current On Rent",
       caption = "Source: HH Pulse and COVID State Housing Policy Score"
       ) +
# scale_color_discrete(name="Presence of Children in Household",
#                       breaks=c("1", "0"),
#                       labels=c("Children in Household", "No Children in Household")
#                       ) +
#  scale_fill_discrete(name="Presence of Children in Household",
#                      breaks=c("1", "0"),
#                      labels=c("Children in Household", "No Children in Household")
#                      ) +
  xlab(label = "State Score") +
  ylab(label = "Predicted Probability of High Eviction Risk") + # ylim(0, 0)
  png('m18.png')

Code

plot_m18

Code

ggsave("plot_m18.png", path = "./figures/")

Saving 7 x 5 in image