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

Question

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)) {
    download.file('https://github.com/dhicks/obesity/blob/mwe/2016-05-24.Rdata?raw=true',
                  datafile)
}
load(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(g$var, 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.

Answer 1

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)))

Answer 2

is it caused by the zeroes in this table? not sure if those will ever converge?

table( x$bmi.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)

Answer 3

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)