Questions regarding ideal approach for filtering out invalid dv values from a regression (and discrepancies with R results)

[kburchfiel]

I often work with dependent-variable values that have missing data (e.g., due to participant non-response). These can cause logistic-regression function calls to fail, since svy (and other libraries) expect to see only 0 and 1 as DV values. Therefore, I wanted to check to see what the best way to handle these missing values would be, and why my svy-based methods produce slightly different standard errors than R-based ones.

Here's some Python code that imports a fictional dataset with a DV column ("Enjoy_Driving_Fast_Strongly_Agree") that has three unique values: 0 for False, 1 for True, and -1 for missing data.

import numpy as np
import pandas as pd
import polars as pl
import svy
import svy
pd.set_option('display.max_columns', 1000)

df_survey = pd.read_csv(
'https://raw.githubusercontent.com/ifstudies/carsurveydata/refs/\
heads/main/car_survey_updated.csv')
df_survey['int_vals'] = df_survey.index

And here's my svy sample object:

# The following code is based on:
# https://svylab.com/docs/svy/tutorials/sample_quicktour.html
design = svy.Design(wgt = 'Weight')
sample = svy.Sample(pl.DataFrame(df_survey), design=design)
print(sample)

I initially tried to handle these invalid values by adding a where argument to my regression code: (Thank you for enabling support for the where parameter, by the way!)

unfiltered_logit_model = sample.glm.fit(y = "Enjoy_Driving_Fast_Strongly_Agree",
                                        x = [svy.Cat('Car_Color')],
                                        family = 'binomial',
                                        link = 'logit',
where = (
svy.col("Enjoy_Driving_Fast_Strongly_Agree") != -1)).to_polars().to_pandas()
unfiltered_logit_model

However, this brings up the following svy error message:

ModelError: 
  ❌ Distribution domain violation [DOMAIN_VIOLATION]
  The target data violates the requirements for family 'binomial': Non-binary target.
  - where: GLM.fit
  - param: family
  - expected: valid range for binomial
  Hint: Must be 0[/1](http://localhost:32885/1)

I found this message surprising, since I did filter out the -1 values for the DV using where. Perhaps there's a way to update the function to accommodate this approach?

Another approach would be to change these -1 values to np.nan, then rerun my logistic-regression code without the where argument. The N/A values appear to get dropped automatically:

# Creating a new sample object in which -1 DV values are replaced with np.nan:
sample_nan = sample.wrangling.mutate({
    "Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries": (
        svy.when(svy.col(
"Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries") == -1).then(
np.nan).otherwise(svy.col(
"Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries")))})

# Running a regression on this new object:
unfiltered_logit_model = sample_nan.glm.fit(
y = "Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries",
x = [svy.Cat('Car_Color')],
family = 'binomial',
link = 'logit'
).to_polars().to_pandas()
unfiltered_logit_model

This produces the following output:

	term	estimate	std_err	conf_low	conf_high	statistic	p_value	df
0	_intercept_	-0.087897	0.120860	-0.325081	0.149288	-0.727259	4.672470e-01	948
1	Car_Color_Red	-0.003292	0.184311	-0.364996	0.358412	-0.017860	9.857540e-01	948
2	Car_Color_White	-1.300545	0.203769	-1.700435	-0.900655	-6.382450	2.723041e-10	948

Interestingly, this appears to produce the same output as the following code (which uses filter_records beforehand to exclude -1 DV values from the dataset):

unfiltered_logit_model = sample.wrangling.filter_records(where = (
svy.col("Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries") != -1)).glm.fit(
y = "Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries",
x = [svy.Cat('Car_Color')],
family = 'binomial',
link = 'logit',
).to_polars().to_pandas()
unfiltered_logit_model

Output:

	term	estimate	std_err	conf_low	conf_high	statistic	p_value	df
0	_intercept_	-0.087897	0.120860	-0.325081	0.149288	-0.727259	4.672470e-01	948
1	Car_Color_Red	-0.003292	0.184311	-0.364996	0.358412	-0.017860	9.857540e-01	948
2	Car_Color_White	-1.300545	0.203769	-1.700435	-0.900655	-6.382450	2.723041e-10	948

However, the standard errors I'm getting here differ very slightly from those obtained using the survey and srvyr libraries in R:

library(tidyverse)
library(survey)
library(srvyr)
library(broom)

df_survey <- read_csv('/home/ifskjb3/Documents/ifsdocs/Scripts/ifs_data_checking/car_survey_updated.csv')

survey_des <- df_survey %>% as_survey_design(
  weights = 'Weight',
  # pps = "brewer",
  # variance = "YG",
  nest = TRUE)

unfiltered_regression_with_minus_1_vals <- survey_des %>% filter(
Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries != -1) %>%  svyglm(design = ., 
formula = Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries ~ Car_Color, family = quasibinomial)
tidy(unfiltered_regression_with_minus_1_vals)

Output of the tidy() call:

term<chr>	estimate<dbl>	std.error<dbl>	statistic<dbl>	p.value<dbl>
(Intercept)	-0.087896649	0.1208537	-0.72729777	4.672231e-01
Car_ColorRed	-0.003291865	0.1843009	-0.01786136	9.857532e-01
Car_ColorWhite	-1.300544807	0.203758	-6.38279275	2.717204e-10

As with svy, I can also replace these -1 values with NA values, then rerun my logistic regression code. I get the same results via this approach that I did within the previous block of R code:

survey_des_nan <- survey_des %>% mutate(Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries = replace_when(
Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries, 
Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries == -1 ~ NA)) 

unfiltered_regression_with_nan_vals <- survey_des_nan %>% svyglm(design = ., 
formula = Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries ~ Car_Color, 
family = quasibinomial) # I think na.action's default value is na.omit, but
# I could be wrong.
tidy(unfiltered_regression_with_nan_vals)

Output:

term<chr>	estimate<dbl>	std.error<dbl>	statistic<dbl>	p.value<dbl>
(Intercept)	-0.087896649	0.1208537	-0.72729777	4.672231e-01
Car_ColorRed	-0.003291865	0.1843009	-0.01786136	9.857532e-01
Car_ColorWhite	-1.300544807	0.203758	-6.38279275	2.717204e-10

Meanwhile, svy's estimation.prop() method appears to produce identical output to R's survey_prop method.

svy code (using where to filter out -1 DV values:

df_svy_props = sample.estimation.prop(y = 'Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries', 
by = 'Car_Color', where = (
svy.col("Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries") != -1)).to_polars().to_pandas()
df_svy_props

Output:

R code:

survey_des %>% filter(Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries != -1) %>% group_by(
Car_Color, Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries) %>% summarize(
p = survey_prop(vartype = c("se", "ci", "var", "cv")))

Output:

This all makes me wonder whether (1) dropping N/A DV values when running logistic regressions in svy might produce slightly inaccurate output (relative to the R-based output I shared above), and (2) whether updating glm.fit to allow its where argument to filter out invalid DV value might allow for logistic regression results that more closely match R-based results.

Thank you as always for your assistance with these requests!

[kburchfiel]

Good morning,

I saw that svy has received a number of updates, so I tried rerunning the following code:

unfiltered_logit_model = sample.glm.fit(y = "Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries",
                                        x = [svy.Cat('Car_Color')],
                                        family = 'binomial',
                                        link = 'logit',
where = (
svy.col("Enjoy_Driving_Fast_Strongly_Agree_with_invalid_entries") != -1)).to_polars().to_pandas()
unfiltered_logit_model

This code now runs without an error and produces point estimates, standard errors, and p-values that appear equivalent to those I'm getting in R. So perhaps the initial issue got resolved?

Thank you for all these updates--they're very helpful!