2
$\begingroup$

Stata's svy command allows to get standard error estimates from two way contingency tables and the tab(varname) option makes it easy to indicate which variable is the total when obtaining proportions, as the documentation here puts it:

tab(varname) specifies that counts be cell totals of this variable and that proportions (or percentages) be relative to (that is, weighted by) this variable. For example, if this variable denotes income, the cell “counts” are instead totals of income for each cell, and the cell proportions are proportions of income for each cell

Following the previous command, a coworker made this calculation, from which he then obtains the SE by making some calculations on the e(V) (variance matrix) object from this table.

enter image description here

In R, using the survey package I can get the mean by groups and afterwards calculate the proportions of those means by their size and get to the same stata output for the proportions, but I'm having a lot of trouble obtaining the SE from said proportions. Note that the SE from the following image is from the mean and not the proportions, which is what i need.

enter image description here

So, how do I get the SE of the proportion estimation of an "aggregated" variable in R?

$\endgroup$
4
$\begingroup$

To expand on the statistical issue: you can't just divide the standard errors by the same denominator that divides the estimates, because the denominator is (in general) random and correlated with the numerators. You need a proper multivariate delta-method calculation.

In R, transformations of the output of a survey function are done with svycontrast. If I've understood your description correctly, you want svytotal rather than svymean

> library(survey)
> data(api)
> dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
> 
> (totals<-svyby(~enroll,~stype,design=dstrat,svytotal,covmat=TRUE))
  stype    enroll       se
E     E 1842584.3 72581.33
H     H  997128.5 69239.40
M     M  847464.7 55502.96
> 
> svycontrast(totals, list(quote(E/(E+H+M)), quote(H/(E+H+M)), quote(M/(E+H+M))))
       nlcon     SE
[1,] 0.49973 0.0155
[2,] 0.27043 0.0153
[3,] 0.22984 0.0132

matching Stata's

. svy: tabulate stype, tab(enroll) se
(running tabulate on estimation sample)

Number of strata   =         3                  Number of obs     =        200
Number of PSUs     =       200                  Population size   =      6,194
                                                Design df         =        197

----------------------------------
    stype | proportion          se
----------+-----------------------
        E |      .4997       .0155
        H |      .2704       .0153
        M |      .2298       .0132
          | 
    Total |          1            
----------------------------------
  Tabulated variable:  enroll

  Key:  proportion  =  cell proportion
        se          =  linearized standard error of cell proportion

Update: the original post didn't match the Stata output because I forgot to declare the finite population correction in Stata

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you so much for your informative response, I got to the same results as my coworker's stata output. $\endgroup$ – David Jorquera Jul 21 at 3:21

Not the answer you're looking for? Browse other questions tagged or ask your own question.