Return to Question

Included an automatic download in the MWE, and a link in the question text to download the MWE

Source Link

edited May 25, 2016 at 10:32

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

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

load(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 difficulty I'm having is illustrated by the following code:

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

load('2016-05-24.Rdata')

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

Source Link

asked May 24, 2016 at 20:21

Dan Hicks

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)

load('2016-05-24.Rdata')

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.