Survey-weighted negative binomial generalised linear model

svyglm.nb() is an extension to the survey-package to fit survey-weighted negative binomial models. It uses svymle to fit sampling-weighted maximum likelihood estimates, based on starting values provided by glm.nb, as proposed by Lumley (2010, pp249).

svyglm.nb(formula, design, ...)

Arguments

formula

An object of class formula, i.e. a symbolic description of the model to be fitted. See 'Details' in glm.

design

An object of class svydesign, providing a specification of the survey design.

...

Other arguments passed down to glm.nb.

Value

An object of class svymle and svyglm.nb, with some additional information about the model.

Details

For details on the computation method, see Lumley (2010), Appendix E (especially 254ff.)

sjstats implements following S3-methods for svyglm.nb-objects: family(), model.frame(), formula(), print(), predict() and residuals(). However, these functions have some limitations:

• family() simply returns the family-object from the underlying glm.nb-model.

• The predict()-method just re-fits the svyglm.nb-model with glm.nb, overwrites the $coefficients from this model-object with the coefficients from the returned svymle-object and finally calls predict.glm to compute the predicted values.

• residuals() re-fits the svyglm.nb-model with glm.nb and then computes the Pearson-residuals from the glm.nb-object.

References

Lumley T (2010). Complex Surveys: a guide to analysis using R. Wiley

Examples

# ------------------------------------------
# This example reproduces the results from
# Lumley 2010, figure E.7 (Appendix E, p256)
# ------------------------------------------
if (require("survey")) {
  data(nhanes_sample)

  # create survey design
  des <- svydesign(
    id = ~SDMVPSU,
    strat = ~SDMVSTRA,
    weights = ~WTINT2YR,
    nest = TRUE,
    data = nhanes_sample
  )

  # fit negative binomial regression
  fit <- svyglm.nb(total ~ factor(RIAGENDR) * (log(age) + factor(RIDRETH1)), des)

  # print coefficients and standard errors
  fit
}
#> Loading required package: survey
#> Loading required package: grid
#> Loading required package: survival
#> 
#> Attaching package: ‘survey’
#> The following object is masked from ‘package:sjstats’:
#> 
#>     cv
#> The following object is masked from ‘package:graphics’:
#> 
#>     dotchart
#>                                   term    irr std.error conf.low conf.high
#> 2                          (Intercept) 9.8463    0.1556   7.2578   13.3580
#> 3                    factor(RIAGENDR)2 0.4511    0.1805   0.3167    0.6426
#> 4                             log(age) 2.9163    0.2331   1.8467    4.6056
#> 5                    factor(RIDRETH1)2 1.0859    0.1477   0.8130    1.4504
#> 6                    factor(RIDRETH1)3 1.0977    0.1779   0.7746    1.5556
#> 7                    factor(RIDRETH1)4 2.2686    0.2974   1.2665    4.0634
#> 8                    factor(RIDRETH1)5 1.0589    0.3789   0.5039    2.2250
#> 9           factor(RIAGENDR)2:log(age) 0.2947    0.2651   0.1753    0.4955
#> 10 factor(RIAGENDR)2:factor(RIDRETH1)2 0.8314    0.2611   0.4984    1.3870
#> 11 factor(RIAGENDR)2:factor(RIDRETH1)3 1.8285    0.1931   1.2523    2.6698
#> 12 factor(RIAGENDR)2:factor(RIDRETH1)4 1.0668    0.3747   0.5119    2.2232
#> 13 factor(RIAGENDR)2:factor(RIDRETH1)5 1.4564    0.4427   0.6116    3.4680
#>       p.value
#> 2  <0.001 ***
#> 3  <0.001 ***
#> 4  <0.001 ***
#> 5  0.5769    
#> 6  0.6003    
#> 7  0.0059 ** 
#> 8  0.8800    
#> 9  <0.001 ***
#> 10 0.4795    
#> 11 0.0018 ** 
#> 12 0.8630    
#> 13 0.3957    
#> 
#> Dispersion parameter Theta: 0.8062
#>    Standard Error of Theta: 0.0216
#> 
#> Showing robust standard errors on link-scale (untransformed).