<p>The <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=survey">survey-package</a> from Thomas Lumley is a great toolkit when analyzing complex samples. It provides <code>svyglm()</code>, to fit generalised linear models to data from a complex survey design. <code>svyglm()</code> covers all families that are also provided by R’s <code>glm()</code> – however, the survey-package has no function to fit negative binomial models, which might be useful for overdispersed count models. Yet, the package provides a generic <code>svymle()</code> to fit user-specified likelihood estimations. In <a target="_blank" rel="nofollow" href="http://r-survey.r-forge.r-project.org/svybook/">his book</a>, Appendix E, Thomas Lumley describes how to write your own likelihood-function, passed to <code>svymle()</code>, to fit <em>negative binomial models</em> for complex samples. So I wrote a small „wrapper“ and implemented a function <code>svyglm.nb()</code> in my <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=sjstats">sjstats-package</a>.</p>
<pre># ------------------------------------------
# This example reproduces the results from
# Lumley 2010, figure E.7 (Appendix E, p256)
# ------------------------------------------
library(sjstats)
library(survey)
data(nhanes_sample)

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

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

# print coefficients and standard errors
round(cbind(coef(fit), survey::SE(fit)), 2)

#>                                          [,1] [,2]
#> theta.(Intercept)                        0.81 0.05
#> eta.(Intercept)                          2.29 0.16
#> eta.factor(RIAGENDR)2                   -0.80 0.18
#> eta.log(age)                             1.07 0.23
#> eta.factor(RIDRETH1)2                    0.08 0.15
#> eta.factor(RIDRETH1)3                    0.09 0.18
#> eta.factor(RIDRETH1)4                    0.82 0.30
#> eta.factor(RIDRETH1)5                    0.06 0.38
#> eta.factor(RIAGENDR)2:log(age)          -1.22 0.27
#> eta.factor(RIAGENDR)2:factor(RIDRETH1)2 -0.18 0.26
#> eta.factor(RIAGENDR)2:factor(RIDRETH1)3  0.60 0.19
#> eta.factor(RIAGENDR)2:factor(RIDRETH1)4  0.06 0.37
#> eta.factor(RIAGENDR)2:factor(RIDRETH1)5  0.38 0.44</pre>
<p>The functions returns an object of class <code>svymle</code>, so all methods provided by the survey-package for this class work – it’s just that there are only a few, and common methods like <code>predict()</code> are currently not implemented. Maybe, hopefully, future updates of the survey-package will include such features.</p>
<p> Tagged: <a target="_blank" rel="nofollow" href="https://strengejacke.wordpress.com/tag/r/">R</a>, <a target="_blank" rel="nofollow" href="https://strengejacke.wordpress.com/tag/rstats/">rstats</a>, <a target="_blank" rel="nofollow" href="https://strengejacke.wordpress.com/tag/statistik/">Statistik</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/strengejacke.wordpress.com/1257/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/strengejacke.wordpress.com/1257/" border="0"></a> <img alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" border="0" data-recalc-dims="1"></p>
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/g0eg5uGbvZw">