I am trying to age standardize using the svystandardize package in R. I
have successfully managed to hit my SUDAAN based targets for estimates by
sex, but not the total. The total is only a little different, but I'd like
some help knowing why it isn't exact. I've included the SUDAAN code that
generates the targets and my R script (and output) that I have so far. I
can't supply the data since that would violate usage agreements.
* Sorry if the output isn't formatted nicely.
SUDAAN results:
Variance Estimation Method: Taylor Series (WR)
Standardized estimates
by: Variable, Sex of respondent.
--------------------------------------------------------------------------------
| | |
Sex of respondent |
| Variable |
|-----------------------------------------|
| | |
Total | Male | Female |
--------------------------------------------------------------------------------
| | |
| |
|
| Current Smoker: | Sample Size | 14732 | 5582
| 9150 |
| At risk | Weighted Size | 18073876.32 |
8917534.53 | 9156341.79 |
| | Total |
3477083.28 | 2105598.41 | 1371484.87 |
| | Percent |
18.84 | 22.78 | 14.87 |
| | SE Percent | 0.60
| 0.95 | 0.73 |
| | Lower 95% Limit |
| | |
| | Percent |
17.69 | 20.97 | 13.50 |
| | Upper 95% Limit |
| | |
| | Percent | 20.05
| 24.69 | 16.34 |
--------------------------------------------------------------------------------
-------------------------
R Output:
Estimate for the total is off.
> svyciprop(~I(rfsmok=="At risk"), stdes, na.rm=TRUE, ci=TRUE,
vartype="ci")
2.5%
97.5%
I(rfsmok == "At risk") 0.187685 0.176309 0.19962
Estimates by sex are a match.
> svyby(~I(rfsmok=="At risk"), ~sex, stdes, svyciprop, na.rm=TRUE, ci=TRUE,
vartype="ci")
sex I(rfsmok == "At risk") ci_l
ci_u
Male Male 0.2277632725 0.2097154285 0.2468790907
Female Female 0.1486524831 0.1349832066 0.1634444481
-------------------------
R Code:
sample.brfss <- svydesign(id=~brfss.2011$SEQNO, strata=~brfss.2011$ststr,
weights=~brfss.2011$LLCPWT, data=brfss.2011, nest=T)
popage <- c(0.128810, 0.182648, 0.219077, 0.299194, 0.170271)
options("survey.lonely.psu"="adjust")
testing_dataset <- subset(sample.brfss, (!is.na(agegr5) & !is.na(sex) & !
is.na(rfsmok)))
stdes<-svystandardize(testing_dataset, by=~agegr5, over=~sex,
population=popage)
-------------------------
SUDAAN code:
PROC DESCRIPT DATA = "h:\\brfss_data\\brfss11\\pudf11\\txbrfss_pudf11"
DESIGN = WR FILETYPE= SUDXPORT NOTSORTED;
NEST STSTR SEQNO/MISSUNIT;
WEIGHT LLCPWT;
SETENV LINESIZE = 106 PAGESIZE = 45 COLWIDTH = 11 TOPMGN = 1;
VAR RFSMOK;
CATLEVEL 2;
STDVAR AGEGR5;
STDWGT .128810 .182648 .219077 .299194 .170271;
SUBGROUP SEX AGEGR5;
LEVELS 2 5 ;
TABLES (SEX);
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
