this worked for me -- and doesn't require removing the PSUs from the
design :)
options( survey.lonely.psu = "adjust" )
svyhist (~dthage,
subset (nhis, xspd2=='No SPD'), breaks=MyBreaks, main= " ",
col="grey80",
xlab="Age at Death Distribution"
)
lines (svysmooth(~dthage, bandwidth=5,subset(nhis, xspd2=='No SPD')), lwd=2)
Dr. Lumley has written quite a bit about single-PSU strata here:
http://faculty.washington.edu/tlumley/survey/exmample-lonely.html
On Fri, Oct 5, 2012 at 7:16 PM, David Winsemius <dwinsem...@comcast.net>wrote:
>
> On Oct 5, 2012, at 3:33 PM, Muhuri, Pradip (SAMHSA/CBHSQ) wrote:
>
> > Hello,
> >
> > I was trying to draw histograms of age at death and got the following
> 2 error messages:
> >
> >
> > 1) Error in tapply(1:NROW(x), list(factor(strata)), function(index) { :
> >
> > arguments must have same length
>
> This is the top of the output of str applied to the data argument you
> offered to svyhist:
>
>
> > str(subset (nhis, xspd2==2) )
> List of 9
> $ cluster :'data.frame': 0 obs. of 1 variable:
> ..$ psu: Factor w/ 47 levels "109.1","115.2",..:
> ..- attr(*, "terms")=Classes 'terms', 'formula' length 2 ~psu
> .. .. ..- attr(*, "variables")= language list(psu)
> .. .. ..- attr(*, "factors")= int [1, 1] 1
> .. .. .. ..- attr(*, "dimnames")=List of 2
> .. .. .. .. ..$ : chr "psu"
> .. .. .. .. ..$ : chr "psu"
>
> At least one problem seems pretty clear. No data. That can be corrected by
> wrapping as.numeric() around the factor on which you are subsetting in two
> places.
>
> Another problem may arise when you restrict to one class only, namely
> there won't any "design" to work with. All the clusters .... there would be
> only one .... no longer have any multiplicity, and svyhist apparently
> isn't built to handle situation, at least with that design argument.
>
> Error in onestrat(x[index, , drop = FALSE], clusters[index],
> nPSU[index][1], :
> Stratum (2) has only one PSU at stage 1
>
> Taking the 'stratum' argument out of the design() spec allows it to
> proceed, but I do not know if that is introducing invalidity in the
> analysis.
> --
> David.
>
> >
> >
> > 2) Error in findInterval(mm[, i], gx) : 'vec' contains NAs
> >
> > In addition: Warning messages:
> >
> > 1: In min(x) : no non-missing arguments to min; returning Inf
> >
> > 2: In max(x) : no non-missing arguments to max; returning -Inf
> >
> >
> >
> > I would appreciate if someone could help me resolve these issues.
> >
> >
> >
> > Below is reproducible example.
> >
> > Thanks,
> >
> > Pradip Muhuri
> >
> >
> >
> > setwd ("E:/RDATA")
> > options(width = 120)
> > library (survey)
> > library (KernSmooth)
> > xd1 <-
> > "dthage ypll_75 xspd2 psu stratum wt8
> > 56 19 2 2 33 1512.7287
> > 86 0 2 2 129 1830.6400
> > 81 0 2 1 67 536.1400
> > 47 28 2 1 17 519.8350
> > 71 4 1 1 225 254.4087
> > 72 3 1 1 238 424.4787
> > 75 0 2 2 115 407.0987
> > 83 0 2 2 46 622.5137
> > 79 -4 2 1 300 509.1212
> > 78 -3 2 1 133 517.3325
> > 71 4 2 2 328 1179.3063
> > 64 11 2 1 2 301.5250
> > 78 -3 2 1 62 253.9025
> > 65 10 2 2 260 932.6575
> > 75 0 2 1 247 145.5900
> > 63 12 2 2 156 247.0650
> > 71 4 2 1 146 829.4787
> > 76 -1 2 2 234 432.5437
> > 76 0 2 1 109 859.6888
> > 68 7 2 1 228 1236.2975
> > 64 11 2 2 167 347.5788
> > 62 13 2 2 312 354.0500
> > 77 0 2 2 275 882.1938
> > 78 -3 2 1 28 481.5975
> > 81 0 2 1 180 1285.5425
> > 79 0 2 2 205 576.0000
> > 70 5 2 1 173 128.3725
> > 75 0 2 2 189 359.3863
> > 78 0 2 1 332 512.8062
> > 74 1 2 2 14 449.0800
> > 77 0 2 1 242 283.0013
> > 92 0 2 1 152 915.3200
> > 69 6 2 2 217 672.7663
> > 53 22 2 1 290 1430.8812
> > 81 0 2 2 90 699.1075
> > 67 8 2 2 316 607.6500
> > 85 0 2 1 171 312.9850
> > 93 0 2 2 119 936.1275
> > 82 0 2 1 118 186.4450
> > 71 4 2 2 329 729.1213
> > 43 32 2 1 215 887.6313
> > 74 1 2 1 180 569.9338
> > 89 0 2 1 324 1054.0887
> > 81 0 2 2 47 532.0987
> > 70 5 2 1 53 450.8750
> > 75 0 1 1 38 557.9750
> > 56 19 2 1 17 512.6363
> > 90 0 2 2 29 569.7888
> > 70 5 2 1 251 554.2138
> > 56 19 2 2 14 1114.1762"
> > tor <- read.table (textConnection(xd1), header=TRUE, sep='', as.is=TRUE)
> >
> >
> > # Grouping variable (xspd) to be factor
> > tor <- within(tor, {
> > xspd2 <- factor(xspd2,levels=c (1,2),
> > labels=c('SPD', 'No SPD'), ordered=TRUE)
> > }
> > )
> > # object with survey design variables and data
> > nhis <- svydesign (id=~psu,strat=~stratum, weights=~wt8, data=tor,
> nest=TRUE)
> >
> > MyBreaks <- c(18,35,45,55,65,75,85,95)
> >
> > png("svyhist_age_at_death.png")
> >
> > svyhist (~dthage,
> > subset (nhis, xspd2==2), breaks=MyBreaks, main= " ",
> > col="grey80",
> > xlab="Age at Death Distribution"
> > )
> > lines (svysmooth(~dthage, bandwidth=5,subset(nhis, xspd2==2)), lwd=2)
> >
> > dev.off ()
> >
> >
> >
> >
> >
> >
> > Pradip K. Muhuri
> > Statistician
> > Substance Abuse & Mental Health Services Administration
> > The Center for Behavioral Health Statistics and Quality
> > Division of Population Surveys
> > 1 Choke Cherry Road, Room 2-1071
> > Rockville, MD 20857
> >
> > Tel: 240-276-1070
> > Fax: 240-276-1260
> > e-mail: pradip.muh...@samhsa.hhs.gov<mailto:pradip.muh...@samhsa.hhs.gov
> >
> >
> > The Center for Behavioral Health Statistics and Quality your feedback.
> Please click on the following link to complete a brief customer survey:
> http://cbhsqsurvey.samhsa.gov<http://cbhsqsurvey.samhsa.gov/>
> >
> >
> > [[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.
>
> David Winsemius, MD
> Alameda, CA, USA
>
> ______________________________________________
> 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.
>
[[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.
