Difficulties fitting a Cox PH model with categorical interactions to complex survey data

I'm attempting to fit a Cox proportional hazards model to a set of NHANES data; the code to load and clean the data is here, and the resulting dataset is here.

The difficulty I'm having is illustrated by the following code:

library(dplyr)
library(survey)
library(survival)

datafile = '2016-05-24.Rdata'
if (!file.exists(datafile)) {
datafile)
}

df = df[complete.cases(df),]

## ----------
## W/o sampling weights
df.working = df %>% filter(!smoker, age.months >= 50*12, age.months < 85*12)
coxfit.unweighted = coxph(Surv(age.months, mort.status == 'deceased') ~
bmi.cat * bmi.max.cat,
data = df.working)
summary(coxfit.unweighted)

## ----------
## W/ sampling weights
## cf the svydesign call in the docs for survey::nhanes
design = svydesign(id = ~ psu, strata = ~ stratum, weights = ~ sample.weight,
nest = TRUE,
data = {df %>% filter(!smoker, age.months >= 50*12, age.months < 85*12)})
coxfit.weighted = svycoxph(Surv(age.months, mort.status == 'deceased') ~
bmi.cat * bmi.max.cat,
design = design)


The variable bmi.cat is a factor variable for body mass index categories (underweight, normal weight, overweight, etc.). The variable bmi.max.cat is a factor variable for the same categories for the respondent's maximum recalled weight. Since bmi.cat is necessarily less than or equal to bmi.max.cat, there are lots of empty cells in the crosstab of these two variables.

The chunk of code fitting coxfit.unweighted using survival::coxph raises a warning about singularities (because of all the empty cells), but does fit a model. It's unweighted, and so I need to use survey::svycoxph to get unbiased estimates. But the call fitting coxfit.weighted returns

Error in solve.default(gvar, coef(g)) : Lapack routine dgesv: system is exactly singular: U[9,9] = 0  Background: I'm working towards replicating the analysis in this paper, table 3. I've omitted the other covariates from the MWE for the sake of simplicity. I have no previous experience with either survival analysis or analyzing data from complex surveys like NHANES, so radical course corrections to my analysis approach would be fine. • hi, it might be easier for others to help you troubleshoot if you include the download in your reproducible example so we can work with the same df object as you. examples of download automation here github.com/ajdamico/asdfree/tree/master/… – Anthony Damico May 25 '16 at 2:46 • Okay, there you go – Dan Hicks May 25 '16 at 10:33 • yo, not sure if your creation of the design aligns with cdc estimates, you should review the replication scripts before you go further.. dplyr+survey together often create bugs, try to avoid unless you're using the srvyr package. naming an object df conflicts with the function, always avoid that.thanks – Anthony Damico May 25 '16 at 12:03 • Thanks, I hadn't realized there was a conflict with df – Dan Hicks May 25 '16 at 12:18 3 Answers This is an old post, but there is a flaw both in the code and the solution that one should be aware of. When you define design = svydesign(id = ~ psu, strata = ~ stratum, weights = ~ sample.weight, nest = TRUE, data = {df %>% filter(!smoker, age.months >= 50*12, age.months < 85*12)})  you are getting rid of all the population you filter out when defining the data= argument. CDC suggests always to define the design by including the whole dataset and then using subset from the survey package like: design <- svydesign(id = ~ psu, strata = ~ stratum, weights = ~ sample.weight, nest = TRUE, data = df) bmi_design <- subset(design, all(!smoker, age.months >= 50*12, age.months < 85*12)))  is it caused by the zeroes in this table? not sure if those will ever converge? table( xbmi.cat , x$bmi.max.cat )  if you use fake data x$bmi.max.cat <- sample(1:5,nrow(x),replace=T)

design <- svydesign( ~ psu , strat= ~ stratum , weights = ~ sample.weight , nest = TRUE , data = subset( x , sample.weight > 0 ) )
w <- subset( design , !smoker & age.months >= 50*12 & age.months < 85*12 )
svycoxph(Surv(age.months, mort.status == 'deceased') ~ bmi.cat * bmi.max.cat, design = w)

• The study I'm trying to reproduce estimated proportional hazards at each combination of bmi.cat and bmi.max.cat. This svycoxph specification doesn't do that. Without coercing the factors to numerics, the same error still appears. – Dan Hicks May 25 '16 at 12:22
• you're right, sorry, i'm not sure of the answer. i don't think those zero cells will ever work (but could be wrong) – Anthony Damico May 25 '16 at 12:35

The solution is to use the interaction function with drop = TRUE to remove the empty levels of the interaction term. Here's how it looks in context:

## W/ sampling weights
## cf the svydesign call in the docs for survey::nhanes
dataf$bmi.cat.inter = interaction(dataf$bmi.cat, dataf\$bmi.max.cat, drop = TRUE)
design = svydesign(id = ~ psu, strata = ~ stratum, weights = ~ sample.weight,
nest = TRUE,
data = {dataf %>% filter(!smoker, age.months >= 50*12, age.months < 85*12)})
coxfit.weighted = svycoxph(Surv(age.months, mort.status == 'deceased') ~
bmi.cat.inter,
design = design)
summary(coxfit.weighted)