Solved – Survey regression in R with singleton PSUs

rregressionsurvey

I am completely new to R, just downloaded and installed it today. I am familiar with SAS and Stata; I am using R because I have found out that in survey regression analysis, R is capable of using data that have stratum with one PSU. However, I cannot figure out how to write the code at all.

Here is what I have done so far: read a Stata dataset and save the .RData file. I have also put in the MASS, pscl, and survey (for svyglm) packages.

Here's what I need to do:
1) I am using survey data, so I have a "weight" variable, a "strata" variable, and a "PSU" variable. I need to incorporate those; I know how to use svyset in Stata, but no idea in R.
2) I have stratum with singleton PSUs. I need to use an option called survey.lonely.psu I believe, and I have no idea where to even begin with that. This is the reason why I am using R instead of Stata as I do not want to collapse stratum or delete observations.
3) The types of regression models I have to run: survey negative binomial, survey zero-inflated negative binomial (need to also determine the predictors of zeros), survey logistic, and survey OLS regression.
4)I also really can't make much sense in R of how to write the model in R code. In Stata, I can simply write the model as:

svy: nbreg dependent_var independent_var1 independent_var2 independent_var3

I can't figure out how to do that at all in R.

Any and all help will be greatly appreciated.

Best Answer

You need to install the survey package. Here is an example of how to define the survey design you have specified and how to run a linear regression on these data. I assume that the dataset has already been loaded.

require(survey)
options(survey.lonely.psu = "adjust")
design1 <- svydesign(id = ~psuid, strata = ~stratvar, weights = ~weightvar, data = mydata)
model1 <- svyglm(y ~ x1 + x2, design = design1)
summary(model1)

IMHO, Thomas Lumley's homepage is an excellent starting point for this kind of things.

Rather than only installing the survey package, you can install the Official Statistics task view:

install.packages("ctv")
install.views("OfficialStatistics")

This task view gives you a rather nice and complete toolbox to work with survey data.

Note that with Stata's svyset command you have basically the same possibilities than you have in R to handle singleton sampling units.

Related Solutions

Solved – Dealing with singleton strata in survey analysis

Also looking for an answer to this. The closest I've found so far is from the R survey package manual:

Using "adjust" [centered in Stata] is conservative, and it would often be better to combine strata in some intelligent way. The properties of "average" [scaled] have not been investigated thoroughly, but it may be useful when the lonely PSUs are due to a few strata having PSUs missing completely at random.

The "remove" and "certainty" options [certainty] give the same result, but "certainty" is intended for situations where there is only one PSU in the population stratum, which is sampled with certainty (also called ‘self-representing’ PSUs or strata). With "certainty" no warning is generated for strata with only one PSU. Ordinarily, svydesign will detect certainty PSUs, making this option unnecessary.

Other sources (e.g. here) seem to suggest that selecting certainty is only appropriate when there is only a single "large" PSU in this stratum in the population, but also note that using approximate methods can give values that are too high.

Solved – Cluster-robust SE in Stata when using a survey design

A standard error for regress with a cluster() option will always be larger than that from svy: regress, with the ratio of squared standard errors equal to $\dfrac{n-1}{n-k}$, where $n$ is the sample size and $k$ is the number of predictors, including the intercept.

Here's a modification of your example to demonstrate this. The change in definition of cluster should make no difference.

sysuse auto, clear
set seed 92122
gen mkr = substr(make,1,2)
codebook mkr // 23 clusters
*run the "manual" regression
reg price mpg weight [pw = gear], cluster(mkr)
scalar v1 = _se[mpg]^2
scalar k = e(rank) // number of predictors + 1
scalar n = e(N)

*using svy design
svyset mkr [pw = gear]
svy : reg price mpg weight
scalar v2 = _se[mpg]^2
di " v1 = " v1   "    v2 = " v2
di " n =  "n   "    k = " k
di " v1/v2 = "v1/v2
di " (n-1)/(n-k) = "(n-1)/(n-k)

The results of the display statements:

v1 = 5364.5846 v2 = 5217.6097
n = 74 k = 3
v1/v2 = 1.028169
(n-1)/(n-k) = 1.028169

The values in the last two lines are identical. Page 469 of the Stata 14 Manual entry for _robust, (http://www.stata.com/manuals14/p_robust.pdf) refers to a multiplier $\dfrac{n}{n-k}$ that should make the two calculations equivalent. I verified the ratio $\dfrac{n-1}{n-k}$ in other examples. I can't account for the difference.

Best Answer

Related Solutions

Solved – Dealing with singleton strata in survey analysis

Solved – Cluster-robust SE in Stata when using a survey design

Related Question