Failed to obtain standard errors of predictors using svyglm() in R package "survey"

Question

I want to fit models to estimate the odds ratios and 95% confidence intervals using svyglm() in the R package "survey". The models are failed to obtain Std. Error of each predictors, showing as Inf, in my data. Similar issue has also been mentioned before (https://stackoverflow.com/questions/42698981/svyglm-in-package-survey-in-r-not-returning-std-errors and https://stat.ethz.ch/pipermail/r-help/2016-November/442870.html), but remains to be addressed. Any suggestions and comments are appreciated.

# --------- data structure

load("D:/excisedata/data1.RData")

str(mydata)

'data.frame':   6508 obs. of  13 variables:
  $ wt       :Class 'labelled' num  8987 5587 26771 35316 5921 ...
.. .. LABEL: design:sample weight 
$ psu      :Class 'labelled' int  1 1 1 2 1 2 1 1 1 1 ...
.. .. LABEL: design:PSU 
$ strat    :Class 'labelled' int  52 51 48 52 51 51 50 44 44 44 ...
.. .. LABEL: design:stratum 
$ age      :Class 'labelled' int  11 15 44 70 16 14 11 19 10 7 ...
.. .. LABEL: Age (years) 
$ sex      :Class 'labelled' Factor w/ 2 levels "Male","Female": 2 1 2 1 2 2 2 1 2 2 ...
.. .. LABEL: Gender 
$ race     :Class 'labelled' Factor w/ 4 levels "Non-Hispanic White",..: 2 2 2 1 2 2 1 1 1 2 ...
.. .. LABEL: Race/ethnicity 
$ edu      :Class 'labelled' Factor w/ 3 levels "Less than high school",..: 1 3 3 3 1 1 2 3 1 2 ...
.. .. LABEL: Education 
$ sala     :Class 'labelled' Factor w/ 2 levels "<= 1","> 1": 1 2 2 2 2 1 2 2 2 1 ...
.. .. LABEL: salary 
$ bmi_cat  :Class 'labelled' Factor w/ 3 levels "Normal","Overweight",..: 1 2 3 1 1 2 1 1 1 2 ...
.. .. LABEL: BMI categories 
$ cotin_cat:Class 'labelled' Factor w/ 3 levels "Low","Medium",..: 2 1 1 1 1 2 2 2 2 2 ...
.. .. LABEL: Serum cotinine categories 
$ cal      :Class 'labelled' int  1402 4110 1458 2168 1688 2866 1040 2232 2134 903 ...
.. .. LABEL: Dietary calories (kcal) 
$ treat    : Factor w/ 3 levels "1","2","3": 3 1 1 3 2 2 3 3 2 1 ...
..- attr(*, "label")= chr "1-low,2-medium,3-high"
$ disease  : num  0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "0-negative,1-positive"

# --------- survey design

library("survey")
sampdesign <- svydesign(id=~psu,
                        strata=~strat,
                        weights=~wt,
                        nest=TRUE,
                        data=mydata)

# --------- model 1, failed

fit<-svyglm(disease~treat+age+sex+race+edu+sala+bmi_cat+cotin_cat+cal,family="binomial",design=sampdesign)
summary(fit)


Call:
  svyglm(formula = disease ~ treat + age + sex + race + edu + sala + 
           bmi_cat + cotin_cat + cal, design = sampdesign, family = "binomial")

Survey design:
  svydesign(id = ~psu, strata = ~strat, weights = ~wt, nest = TRUE, 
            data = mydata)

Coefficients:
                            Estimate           Std. Error
(Intercept)                  -8.910e-01        Inf
treat2                       -7.455e-02        Inf
treat3                       -7.125e-02        Inf
age                          -6.100e-03        Inf
sexFemale                    -5.867e-01        Inf
raceNon-Hispanic Black        8.896e-01        Inf
raceHispanic                  5.157e-01        Inf
raceOthers                    6.365e-01        Inf
eduHigh school or equivalent -1.848e-02        Inf
eduAbove high school         -7.933e-02        Inf
sala> 1                      -1.803e-01        Inf
bmi_catOverweight             1.036e-01        Inf
bmi_catObese                  1.944e-01        Inf
cotin_catMedium               5.323e-02        Inf
cotin_catHigh                 1.998e-01        Inf
cal                          -4.999e-05        Inf

(Dispersion parameter for binomial family taken to be 1.000021)

Number of Fisher Scoring iterations: 4

Answer 1

I think (though you don't show enough to be sure) that you have zero or negative residual degrees of freedom in your model. The residual df is the design df (the number of PSUs minus the number of strata) minus the number of predictors, which can easily get negative when you have two large clusters per stratum. Having zero or negative df makes inference tricky -- there is an argument that this df estimate is conservative, but there isn't a good solution.

You can extract the standard errors with

SE(fit)

and if you want to use a different residual degrees of freedom, you can specify that to summary and get $p$-values. In particular, if none of your covariates are at the cluster level, you may be able to use

summary(fit, df=degf(sampdesign))

In the forthcoming version 4.1 the package will report standard errors in this situation (but not $p$-values unless a different df= is specified)