R Survey Package

http://r-survey.r-forge.r-project.org/survey/
Admin
The R Survey Package provides tools for statistical analysis of complex survey data, including...
Tokens:66,071
Snippets:516
Trust Score:9.9
Update:3 months ago
Show doc for...
Tokens:
Raw
### Multistage Survey Design Examples

Source: http://r-survey.r-forge.r-project.org/survey/html/svyrecvar.html

Examples demonstrating the setup of two-stage cluster samples and variance estimation using the survey package.

```R
data(mu284)
dmu284<-svydesign(id=~id1+id2,fpc=~n1+n2, data=mu284)
svytotal(~y1, dmu284)

data(api)
# two-stage cluster sample
dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
summary(dclus2)
svymean(~api00, dclus2)
svytotal(~enroll, dclus2,na.rm=TRUE)

# two-stage `with replacement'
dclus2wr<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
summary(dclus2wr)
svymean(~api00, dclus2wr)
svytotal(~enroll, dclus2wr,na.rm=TRUE)
```

--------------------------------

### Compute replicate-based variance estimates

Source: http://r-survey.r-forge.r-project.org/survey/html/withReplicates.html

Example demonstrating the setup of a survey design with replicate weights and the calculation of a ratio using both expression and function interfaces.

```R
data(scd)
repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
c(0,1,0,1,1,0))
scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)

a<-svyratio(~alive, ~arrests, design=scdrep)
print(a$ratio)
print(a$var)
withReplicates(scdrep, quote(sum(.weights*alive)/sum(.weights*arrests)))
withReplicates(scdrep, function(w,data) sum(w*data$alive)/sum(w*data$arrests))
```

--------------------------------

### Initialize a database-backed survey design

Source: http://r-survey.r-forge.r-project.org/survey/svy-dbi.html

Setup a survey design object using an SQLite database connection.

```r
> library(survey)
> library(RSQLite)
Loading required package: DBI
> 
> dhanes<-svydesign(id=~SDPPSU6, strat=~SDPSTRA6, weight=~WTPFQX6,  nest=TRUE, data="set1", dbname='~/nhanes/imp.db', dbtype="SQLite")
>
> dhanes
DB-backed Stratified 1 - level Cluster Sampling design (with replacement)
With (98) clusters.
svydesign(ids = ids, probs = probs, strata = strata, data = design.data, 
    fpc = fpc, variables = variables, nest = nest, check.strata = check.strata, 
    weights = weights)
```

--------------------------------

### Factor analysis examples

Source: http://r-survey.r-forge.r-project.org/survey/html/svyfactanal.html

Examples demonstrating factor analysis on survey data and comparison with population data.

```R
data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)

svyfactanal(~api99+api00+hsg+meals+ell+emer, design=dclus1, factors=2)

svyfactanal(~api99+api00+hsg+meals+ell+emer, design=dclus1, factors=2, n="effective")

##Population dat for comparison
factanal(~api99+api00+hsg+meals+ell+emer, data=apipop, factors=2)
```

--------------------------------

### Initialize an ODBC-backed survey design

Source: http://r-survey.r-forge.r-project.org/survey/svy-dbi.html

Setup a survey design object using an ODBC data source name.

```r
odhanes<-svydesign(id=~SDPPSU6, strat=~SDPSTRA6, weight=~WTPFQX6,  nest=TRUE, data="set1", dbname='nhanes3', dbtype="ODBC")
```

--------------------------------

### Examples of barplot and dotchart usage

Source: http://r-survey.r-forge.r-project.org/survey/html/barplot.svystat.html

Demonstrates creating barplots and dotcharts from survey design objects and statistical summaries.

```R
data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)

a<-svymean(~stype, dclus1)
barplot(a)
barplot(a, names.arg=c("Elementary","High","Middle"), col="purple", main="Proportions of school level")

b<-svyby(~enroll+api.stu, ~stype, dclus1, svymean)
barplot(b,beside=TRUE,legend=TRUE)
dotchart(b)
```

--------------------------------

### Comparing pchisqsum approximation methods

Source: http://r-survey.r-forge.r-project.org/survey/html/pchisqsum.html

Example demonstrating the use of different approximation methods (integration, Satterthwaite, and saddlepoint) and plotting the results.

```R
x <- 5*rnorm(1001)^2+rnorm(1001)^2
x.thin<-sort(x)[1+(0:100)*10]
p.invert<-pchisqsum(x.thin,df=c(1,1),a=c(5,1),method="int" ,lower=FALSE)
p.satt<-pchisqsum(x.thin,df=c(1,1),a=c(5,1),method="satt",lower=FALSE)
p.sadd<-pchisqsum(x.thin,df=c(1,1),a=c(5,1),method="sad",lower=FALSE)

plot(p.invert, p.satt,type="l",log="xy")
abline(0,1,lty=2,col="purple")
plot(p.invert, p.sadd,type="l",log="xy")
abline(0,1,lty=2,col="purple")
```

--------------------------------

### Load API data

Source: http://r-survey.r-forge.r-project.org/survey/example-design.html

Load the example API dataset for survey design analysis.

```R
  data(api)
```

--------------------------------

### Examples of svyttest

Source: http://r-survey.r-forge.r-project.org/survey/html/svyttest.html

Demonstrates two-sample and one-sample t-tests using a survey design object.

```R
data(api)
dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
svyttest(enroll~comp.imp, dclus2)

svyttest(I(api00-api99)~0, dclus2)
```

--------------------------------

### Raking survey designs with examples

Source: http://r-survey.r-forge.r-project.org/survey/html/rake.html

Comprehensive examples demonstrating raking on replicate weight designs, standard designs, and comparisons with calibration and post-stratification.

```R
data(api)
dclus1 <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
rclus1 <- as.svrepdesign(dclus1)

svymean(~api00, rclus1)
svytotal(~enroll, rclus1)

## population marginal totals for each stratum
pop.types <- data.frame(stype=c("E","H","M"), Freq=c(4421,755,1018))
pop.schwide <- data.frame(sch.wide=c("No","Yes"), Freq=c(1072,5122))

rclus1r <- rake(rclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))

svymean(~api00, rclus1r)
svytotal(~enroll, rclus1r)

## marginal totals correspond to population
xtabs(~stype, apipop)
svytable(~stype, rclus1r, round=TRUE)
xtabs(~sch.wide, apipop)
svytable(~sch.wide, rclus1r, round=TRUE)

## joint totals don't correspond 
xtabs(~stype+sch.wide, apipop)
svytable(~stype+sch.wide, rclus1r, round=TRUE)

## Do it for a design without replicate weights
dclus1r<-rake(dclus1, list(~stype,~sch.wide), list(pop.types, pop.schwide))

svymean(~api00, dclus1r)
svytotal(~enroll, dclus1r)

## compare to raking with calibrate()
dclus1gr<-calibrate(dclus1, ~stype+sch.wide, pop=c(6194, 755,1018,5122),
           calfun="raking")
svymean(~stype+api00, dclus1r)
svymean(~stype+api00, dclus1gr)

## compare to joint post-stratification
## (only possible if joint population table is known)
##
pop.table <- xtabs(~stype+sch.wide,apipop)
rclus1ps <- postStratify(rclus1, ~stype+sch.wide, pop.table)
svytable(~stype+sch.wide, rclus1ps, round=TRUE)

svymean(~api00, rclus1ps)
svytotal(~enroll, rclus1ps)

## Example of raking with partial joint distributions
pop.imp<-data.frame(comp.imp=c("No","Yes"),Freq=c(1712,4482))
dclus1r2<-rake(dclus1, list(~stype+sch.wide, ~comp.imp),
               list(pop.table, pop.imp))
svymean(~api00, dclus1r2)
svytotal(~enroll, dclus1r2)
```

--------------------------------

### Create survey plots

Source: http://r-survey.r-forge.r-project.org/survey/html/svyplot.html

Examples demonstrating different plotting styles including bubble, transparent, hexagonal binning, and subsampling.

```R
data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)

svyplot(api00~api99, design=dstrat, style="bubble")
svyplot(api00~api99, design=dstrat, style="transparent",pch=19)
## Not run: 
## these two require the hexbin package from Bioconductor
svyplot(api00~api99, design=dstrat, style="hex", xlab="1999 API",ylab="2000 API")
svyplot(api00~api99, design=dstrat, style="grayhex",legend=0)
## End(Not run)

dclus2<-svydesign(id=~dnum+snum,  weights=~pw,
                    data=apiclus2, fpc=~fpc1+fpc2)
svyplot(api00~api99, design=dclus2, style="subsample")
svyplot(api00~api99, design=dclus2, style="subsample",
          amount=list(x=25,y=25))

svyplot(api00~api99, design=dstrat,
  basecol=function(df){c("goldenrod","tomato","sienna")[as.numeric(df$stype)]},
  style="transparent",pch=19,alpha=c(0,1))
legend("topleft",col=c("goldenrod","tomato","sienna"), pch=19, legend=c("E","H","M"))
```

--------------------------------

### Example of compressing replicate weights

Source: http://r-survey.r-forge.r-project.org/survey/html/compressWeights.html

Demonstrates creating a survey design object and comparing compressed versus uncompressed replicate weight designs.

```R
data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
rclus1c<-as.svrepdesign(dclus1,compress=TRUE)
rclus1<-as.svrepdesign(dclus1,compress=FALSE)
```

--------------------------------

### Ratio Estimation Examples

Source: http://r-survey.r-forge.r-project.org/survey/html/svyratio.html

Demonstration of ratio estimation using various survey designs, including replicate weights and stratified samples.

```R
data(scd)

## survey design objects
scddes<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
nest=TRUE, fpc=rep(5,6))
scdnofpc<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
nest=TRUE)

# convert to BRR replicate weights
scd2brr <- as.svrepdesign(scdnofpc, type="BRR")

# use BRR replicate weights from Levy and Lemeshow
repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1),
c(0,1,0,1,1,0))
scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)

# ratio estimates
svyratio(~alive, ~arrests, design=scddes)
svyratio(~alive, ~arrests, design=scdnofpc)
svyratio(~alive, ~arrests, design=scd2brr)
svyratio(~alive, ~arrests, design=scdrep)

data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)

## domain means are ratio estimates, but available directly
svyratio(~I(api.stu*(comp.imp=="Yes")), ~as.numeric(comp.imp=="Yes"), dstrat)
svymean(~api.stu, subset(dstrat, comp.imp=="Yes"))

## separate and combined ratio estimates of total
(sep<-svyratio(~api.stu,~enroll, dstrat,separate=TRUE))
(com<-svyratio(~api.stu, ~enroll, dstrat))

stratum.totals<-list(E=1877350, H=1013824, M=920298)

predict(sep, total=stratum.totals)
predict(com, total=sum(unlist(stratum.totals)))

SE(com)
coef(com)
coef(com, drop=FALSE)
confint(com)
```

--------------------------------

### Create conditioned survey plots

Source: http://r-survey.r-forge.r-project.org/survey/html/svycoplot.html

Examples demonstrating the use of hexbin and transparent styles for conditioned survey data visualization.

```R
data(api)
dclus2<-svydesign(id=~dnum+snum,  weights=~pw,
                    data=apiclus2, fpc=~fpc1+fpc2)

svycoplot(api00~api99|sch.wide*comp.imp, design=dclus2, style="hexbin")
svycoplot(api00~api99|sch.wide*comp.imp, design=dclus2, style="hexbin", hexscale="absolute")

svycoplot(api00~api99|sch.wide, design=dclus2, style="trans")

svycoplot(api00~meals|stype,design=dclus2,
  style="transparent",
  basecol=function(d) c("darkred","purple","forestgreen")[as.numeric(d$stype)],
  alpha=c(0,1))
```

--------------------------------

### Visualize Hadamard matrices and size distribution

Source: http://r-survey.r-forge.r-project.org/survey/html/hadamard.html

Examples demonstrating the construction of different types of Hadamard matrices and plotting the resulting matrix sizes against the input n.

```R
par(mfrow=c(2,2))
## Sylvester-type
image(hadamard(63),main=quote("Sylvester: "*64==2^6))
## Paley-type
image(hadamard(59),main=quote("Paley: "*60==59+1))
## from NJ Sloane's library
image(hadamard(27),main=quote("Stored: "*28))
## For n=90 we get 96 rather than the minimum possible size, 92.
image(hadamard(90),main=quote("Constructed: "*96==2^3%*%(11+1)))

par(mfrow=c(1,1))
plot(2:150,sapply(2:150,function(i) ncol(hadamard(i))),type="S",
     ylab="Matrix size",xlab="n",xlim=c(1,150),ylim=c(1,150))
abline(0,1,lty=3)
lines(2:150, 2:150-(2:150 %% 4)+4,col="purple",type="S",lty=2)
legend(c(x=10,y=140),legend=c("Actual size","Minimum possible size"),
     col=c("black","purple"),bty="n",lty=c(1,2))
```

--------------------------------

### Examples of svyciprop with survey designs

Source: http://r-survey.r-forge.r-project.org/survey/html/svyciprop.html

Demonstrates calculating confidence intervals using different methods on both cluster and replicate survey designs.

```R
data(api)
dclus1<-svydesign(id=~dnum, fpc=~fpc, data=apiclus1)

svyciprop(~I(ell==0), dclus1, method="li")
svyciprop(~I(ell==0), dclus1, method="lo")
svyciprop(~I(ell==0), dclus1, method="as")
svyciprop(~I(ell==0), dclus1, method="be")
svyciprop(~I(ell==0), dclus1, method="me")

rclus1<-as.svrepdesign(dclus1)
svyciprop(~I(emer==0), rclus1, method="li")
svyciprop(~I(emer==0), rclus1, method="lo")
svyciprop(~I(emer==0), rclus1, method="as")
svyciprop(~I(emer==0), rclus1, method="be")
svyciprop(~I(emer==0), rclus1, method="me")
```

--------------------------------

### Survey Analysis with survey Package

Source: http://r-survey.r-forge.r-project.org/survey/html/surveysummary.html

Demonstrates basic survey design setup, statistical estimation, and result extraction using the survey package.

```R

  data(api)

  ## one-stage cluster sample
  dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)

  svymean(~api00, dclus1, deff=TRUE)
  svymean(~factor(stype),dclus1)
  svymean(~interaction(stype, comp.imp), dclus1)
  svyquantile(~api00, dclus1, c(.25,.5,.75))
  svytotal(~enroll, dclus1, deff=TRUE)
  svyratio(~api.stu, ~enroll, dclus1)

  v<-svyvar(~api00+api99, dclus1)
  v
  print(v, cov=TRUE)
  as.matrix(v)

  # replicate weights - jackknife (this is slower)
  dstrat<-svydesign(id=~1,strata=~stype, weights=~pw,
        data=apistrat, fpc=~fpc)
  jkstrat<-as.svrepdesign(dstrat)

  svymean(~api00, jkstrat)
  svymean(~factor(stype),jkstrat)
  svyvar(~api00+api99,jkstrat)

  svyquantile(~api00, jkstrat, c(.25,.5,.75))
  svytotal(~enroll, jkstrat)
  svyratio(~api.stu, ~enroll, jkstrat)

  # coefficients of variation
  cv(svytotal(~enroll,dstrat))
  cv(svyratio(~api.stu, ~enroll, jkstrat))

  # extracting information from the results
  coef(svytotal(~enroll,dstrat))
  vcov(svymean(~api00+api99,jkstrat))
  SE(svymean(~enroll, dstrat))
  confint(svymean(~api00+api00, dclus1))

  # Design effect
  svymean(~api00, dstrat, deff=TRUE)
  svymean(~api00, dstrat, deff="replace")
  svymean(~api00, jkstrat, deff=TRUE)
  svymean(~api00, jkstrat, deff="replace")
 (a<-svytotal(~enroll, dclus1, deff=TRUE))
  deff(a)

 
```

--------------------------------

### Executing survey-sampled MLE

Source: http://r-survey.r-forge.r-project.org/survey/example-programming.html

Run the svymle function with the defined objective, gradient, design, and starting parameters.

```R
> summary(svymle(loglike = lcens, gradient = gcens, design = dpbc, 
+     formulas = list(mean = I(cbind(time, status)) ~ bili + protime + 
+         alb, sd = ~1), start = list(c(10, 0, 0, 0), c(1))))
Survey-sampled mle: 
svymle(loglike = lcens, gradient = gcens, design = dpbc, formulas = list(mean = I(cbind(time, 
    status)) ~ bili + protime + alb, sd = ~1), start = list(c(10, 
    0, 0, 0), c(1)))
                       Coef         SE p.value
mean.(Intercept)  8.9803185 1.27851656  <0.001
mean.bili        -0.0927236 0.01648777  <0.001
mean.protime     -0.3689851 0.09495636  <0.001
mean.alb          0.9535612 0.17204268  <0.001
sd.(Intercept)    0.9888172 0.07561830  <0.001
Stratified Independent Sampling design
svydesign(id = ~1, prob = ~randprob, strata = ~edema, data = subset(pbc, 
    trt > 0))
```

--------------------------------

### Visualize survey data with histograms and boxplots

Source: http://r-survey.r-forge.r-project.org/survey/html/svyhist.html

Example demonstrating weighted histogram and boxplot creation using the apistrat dataset.

```R
data(api)
dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
    fpc = ~fpc)
opar<-par(mfrow=c(1,3))
svyhist(~enroll, dstrat, main="Survey weighted",col="purple",ylim=c(0,1.3e-3))
hist(apistrat$enroll,  main="Sample unweighted",col="purple",prob=TRUE,ylim=c(0,1.3e-3))
hist(apipop$enroll,  main="Population",col="purple",prob=TRUE,ylim=c(0,1.3e-3))

par(mfrow=c(1,1))
svyboxplot(enroll~stype,dstrat)
svyboxplot(enroll~1,dstrat)
par(opar)
```

--------------------------------

### Post-stratification Example

Source: http://r-survey.r-forge.r-project.org/survey/html/postStratify.html

Demonstrates applying post-stratification to both replicate weight designs and standard survey designs using population frequency data.

```R
data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
rclus1<-as.svrepdesign(dclus1)

svymean(~api00, rclus1)
svytotal(~enroll, rclus1)

# post-stratify on school type
pop.types <- data.frame(stype=c("E","H","M"), Freq=c(4421,755,1018))
#or: pop.types <- xtabs(~stype, data=apipop)
#or: pop.types <- table(stype=apipop$stype)

rclus1p<-postStratify(rclus1, ~stype, pop.types)
summary(rclus1p)
svymean(~api00, rclus1p)
svytotal(~enroll, rclus1p)

## and for svydesign objects
dclus1p<-postStratify(dclus1, ~stype, pop.types)
summary(dclus1p)
svymean(~api00, dclus1p)
svytotal(~enroll, dclus1p)
```

--------------------------------

### Subset survey designs and estimate means

Source: http://r-survey.r-forge.r-project.org/survey/html/subset.survey.design.html

Examples demonstrating how to subset survey designs and perform domain estimation using the survey package.

```R
data(fpc)
dfpc<-svydesign(id=~psuid,strat=~stratid,weight=~weight,data=fpc,nest=TRUE)
dsub<-subset(dfpc,x>4)
summary(dsub)
svymean(~x,design=dsub)

## These should give the same domain estimates and standard errors
svyby(~x,~I(x>4),design=dfpc, svymean)
summary(svyglm(x~I(x>4)+0,design=dfpc))

data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
rclus1<-as.svrepdesign(dclus1)
svymean(~enroll, subset(dclus1, sch.wide=="Yes" & comp.imp=="Yes"))
svymean(~enroll, subset(rclus1, sch.wide=="Yes" & comp.imp=="Yes"))
```

--------------------------------

### Examples of regTermTest usage

Source: http://r-survey.r-forge.r-project.org/survey/html/regTermTest.html

Demonstrates testing terms in both standard glm models and survey-weighted svyglm models.

```R
 data(esoph)
 model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * 
     alcgp, data = esoph, family = binomial())
 anova(model1)

 regTermTest(model1,"tobgp")
 regTermTest(model1,"tobgp:alcgp")
 regTermTest(model1, ~alcgp+tobgp:alcgp)

 data(api)
 dclus2<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
 model2<-svyglm(I(sch.wide=="Yes")~ell+meals+mobility, design=dclus2, family=quasibinomial())
 regTermTest(model2, ~ell)
 regTermTest(model2, ~ell,df=NULL)
 regTermTest(model2, ~ell, method="LRT")
 regTermTest(model2, ~ell+meals, method="LRT")
```

--------------------------------

### Example of HR function usage

Source: http://r-survey.r-forge.r-project.org/survey/html/HR.html

Demonstrates a basic call to the HR function with a psum value.

```R
HR(0.1)
```

--------------------------------

### Fitting and testing loglinear models

Source: http://r-survey.r-forge.r-project.org/survey/html/svyloglin.html

Example workflow for fitting a model, updating it to include interactions, and performing ANOVA and Wald tests.

```R
 data(api)
 dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
 a<-svyloglin(~stype+comp.imp,dclus1)
 b<-update(a,~.^2)
 an<-anova(a,b)
 an
 print(an, pval="saddlepoint")

 ## Wald test
 regTermTest(b, ~stype:comp.imp)

 ## linear-by-linear association
 d<-update(a,~.+as.numeric(stype):as.numeric(comp.imp))
 an1<-anova(a,d)
 an1
```

--------------------------------

### Construct and verify a Hadamard matrix

Source: http://r-survey.r-forge.r-project.org/survey/html/paley.html

Example of generating a matrix with paley and verifying its properties using is.hadamard and matrix multiplication.

```R
M<-paley(11)

is.hadamard(M)
## internals of is.hadamard(M)
H<-2*M-1
## HH^T is diagonal for any Hadamard matrix
H%*%t(H)
```

--------------------------------

### Defining and Using Custom Calibration Metrics

Source: http://r-survey.r-forge.r-project.org/survey/html/make.calfun.html

Examples showing how to inspect existing calibration functions, define a new Hellinger distance calibration function, and apply it to survey design objects.

```R
str(cal.linear)
cal.linear$Fm1
cal.linear$dF

hellinger <- make.calfun(Fm1=function(u, bounds)  ((1-u/2)^-2)-1,
                    dF= function(u, bounds) (1-u/2)^-3 ,
                    name="hellinger distance")

hellinger

data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)

svymean(~api00,calibrate(dclus1, ~api99, pop=c(6194, 3914069),
         calfun=hellinger))

svymean(~api00,calibrate(dclus1, ~api99, pop=c(6194, 3914069),
         calfun=cal.linear))

svymean(~api00,calibrate(dclus1, ~api99, pop=c(6194,3914069),
          calfun=cal.raking))
```

--------------------------------

### Predictive Margins Examples

Source: http://r-survey.r-forge.r-project.org/survey/html/marginpred.html

Demonstrates calculating adjusted marginal predictions for Cox models, Kaplan-Meier estimators, and logistic regression.

```R
## generate data with apparent group effect from confounding
set.seed(42)
df<-data.frame(x=rnorm(100))
df$time<-rexp(100)*exp(df$x-1)
df$status<-1
df$group<-(df$x+rnorm(100))>0
des<-svydesign(id=~1,data=df)
newdf<-data.frame(group=c(FALSE,TRUE), x=c(0,0))

## Cox model
m0<-svycoxph(Surv(time,status)~group,design=des)
m1<-svycoxph(Surv(time,status)~group+x,design=des)
## conditional predictions, unadjusted and adjusted
cpred0<-predict(m0, type="curve", newdata=newdf, se=TRUE)
cpred1<-predict(m1, type="curve", newdata=newdf, se=TRUE)
## adjusted marginal prediction
mpred<-marginpred(m0, adjustfor=~x, predictat=newdf, se=TRUE)

plot(cpred0)
lines(cpred1[[1]],col="red")
lines(cpred1[[2]],col="red")
lines(mpred[[1]],col="blue")
lines(mpred[[2]],col="blue")

## Kaplan--Meier
s2<-svykm(Surv(time,status>0)~group, design=des)
p2<-marginpred(s2, adjustfor=~x, predictat=newdf,se=TRUE)
plot(s2)
lines(p2[[1]],col="green")
lines(p2[[2]],col="green")

## logistic regression
logisticm <- svyglm(group~time, family=quasibinomial, design=des)
newdf$time<-c(0.1,0.8)
logisticpred <- marginpred(logisticm, adjustfor=~x, predictat=newdf)
```

--------------------------------

### Estimate and visualize survival functions

Source: http://r-survey.r-forge.r-project.org/survey/html/svykm.html

Example workflow for creating a survey design, estimating survival curves, and performing diagnostic plots and quantile calculations.

```R
data(pbc, package="survival")
pbc$randomized <- with(pbc, !is.na(trt) & trt>0)
biasmodel<-glm(randomized~age*edema,data=pbc)
pbc$randprob<-fitted(biasmodel)

dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))

s1<-svykm(Surv(time,status>0)~1, design=dpbc)
s2<-svykm(Surv(time,status>0)~I(bili>6), design=dpbc)

plot(s1)
plot(s2)
plot(s2, lwd=2, pars=list(lty=c(1,2),col=c("purple","forestgreen")))

quantile(s1, probs=c(0.9,0.75,0.5,0.25,0.1))

s3<-svykm(Surv(time,status>0)~I(bili>6), design=dpbc,se=TRUE)
plot(s3[[2]],col="purple")

confint(s3[[2]], parm=365*(1:5))
quantile(s3[[1]], ci=TRUE)
```

--------------------------------

### Convert survey design to BRR replicate weights

Source: http://r-survey.r-forge.r-project.org/survey/html/brrweights.html

Examples demonstrating conversion of a survey design to BRR replicate weights and subsequent estimation.

```R
data(scd)
scdnofpc<-svydesign(data=scd, prob=~1, id=~ambulance, strata=~ESA,
nest=TRUE)

## convert to BRR replicate weights
scd2brr <- as.svrepdesign(scdnofpc, type="BRR")
svymean(~alive, scd2brr)
svyratio(~alive, ~arrests, scd2brr)

## with user-supplied hadamard matrix
scd2brr1 <- as.svrepdesign(scdnofpc, type="BRR", hadamard.matrix=paley(11))
svymean(~alive, scd2brr1)
svyratio(~alive, ~arrests, scd2brr1)
```

--------------------------------

### Constructing and adjusting non-response weights

Source: http://r-survey.r-forge.r-project.org/survey/html/nonresponse.html

Example workflow for creating a non-response object, identifying sparse cells, joining cells, and extracting final weights.

```R
data(api)
## pretend the sampling was stratified on three variables
poptable<-xtabs(~sch.wide+comp.imp+stype,data=apipop)
sample.count<-xtabs(~sch.wide+comp.imp+stype,data=apiclus1)
sample.weight<-xtabs(pw~sch.wide+comp.imp+stype, data=apiclus1)

## create a nonresponse object
nr<-nonresponse(sample.weight,sample.count, poptable)

## sparse cells
sparseCells(nr)

## Look at neighbours
neighbours(3,nr)
neighbours(11,nr)

## Collapse some contiguous cells
nr1<-joinCells(nr,3,5,7)

## sparse cells now
sparseCells(nr1)
nr2<-joinCells(nr1,3,11,8)

nr2

## one relatively sparse cell
sparseCells(nr2)
## but nothing suitable to join it to
neighbours(3,nr2)

## extract the weights
weights(nr2)
```

--------------------------------

### bootweights(strata, psu, replicates = 50, fpc = NULL, fpctype = c("population", "fraction", "correction"), compress = TRUE)

Source: http://r-survey.r-forge.r-project.org/survey/html/bootweights.html

Creates bootstrap weights for infinite populations by sampling n PSUs with replacement.

```APIDOC
## bootweights

### Description
Creates bootstrap weights by sampling n PSUs with replacement. This is most useful when strata are large and may have large sampling fractions.

### Arguments
- **strata** - Identifier for sampling strata (top level)
- **psu** - Identifier for primary sampling units
- **replicates** - Number of bootstrap replicates
- **fpc** - Finite population correction
- **fpctype** - Is fpc the population size, sampling fraction, or 1-sampling fraction?
- **compress** - Should the replicate weights be compressed?

### Value
A set of replicate weights
```

--------------------------------

### subbootweights(strata, psu, replicates = 50, compress = TRUE)

Source: http://r-survey.r-forge.r-project.org/survey/html/bootweights.html

Creates bootstrap weights by sampling n-1 PSUs from the n available, useful when strata have a small number of PSUs.

```APIDOC
## subbootweights

### Description
Creates bootstrap weights by sampling n-1 PSUs from the n available (Rao and Wu). This is useful when strata have a small number of PSUs and sampling fractions are small.

### Arguments
- **strata** - Identifier for sampling strata (top level)
- **psu** - Identifier for primary sampling units
- **replicates** - Number of bootstrap replicates
- **compress** - Should the replicate weights be compressed?

### Value
A set of replicate weights
```

--------------------------------

### Database-backed survey design with replication weights

Source: http://r-survey.r-forge.r-project.org/survey/html/svrepdesign.html

Shows how to initialize a survey design object using data from a SQLite database and manage the connection.

```R
## Not run:
## Needs RSQLite
library(RSQLite)
db_rclus1<-svrepdesign(weights=~pw, repweights="wt[1-9]+", type="JK1", scale=(1-15/757)*14/15,
data="apiclus1rep",dbtype="SQLite", dbname=system.file("api.db",package="survey"))
svymean(~api00+api99,db_rclus1)

summary(db_rclus1)

## closing and re-opening a connection
close(db_rclus1)
db_rclus1
try(svymean(~api00+api99,db_rclus1))
db_rclus1<-open(db_rclus1)
svymean(~api00+api99,db_rclus1)


## End(Not run)
```

--------------------------------

### Fitting a Proportional Odds Model

Source: http://r-survey.r-forge.r-project.org/survey/html/svyolr.html

Example of fitting a proportional odds model using the api dataset.

```R
data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
dclus1<-update(dclus1, mealcat=cut(meals,c(0,25,50,75,100)))

svyolr(mealcat~avg.ed+mobility+stype, design=dclus1)
```

--------------------------------

### Configure survey options and load data

Source: http://r-survey.r-forge.r-project.org/survey/example-twostage.html

Sets the global option for handling lonely PSUs and loads the survey dataset.

```R
> library(survey)
> options(survey.lonely.psu="remove")
> nhis<-read.dta("nhis.dta")
```

--------------------------------

### Calibrate replicate weight designs

Source: http://r-survey.r-forge.r-project.org/survey/html/calibrate.html

Shows how to apply calibration to replicate weight designs and visualize the results.

```R
## Same syntax with replicate weights
rclus1<-as.svrepdesign(dclus1)

(rclus1g3 <- calibrate(rclus1, ~stype+api99, c(pop.totals, api99=3914069)))

svymean(~api00, rclus1g3)
svytotal(~enroll, rclus1g3)
svytotal(~stype, rclus1g3)

(rclus1agg3 <- calibrate(rclus1, ~stype+api99, c(pop.totals,api99=3914069), aggregate.index=~dnum))

svymean(~api00, rclus1agg3)
svytotal(~enroll, rclus1agg3)
svytotal(~stype, rclus1agg3)

image(rclus1agg3)
```

--------------------------------

### Updating a Survey Design Object

Source: http://r-survey.r-forge.r-project.org/survey/html/update.survey.design.html

Example demonstrating how to add a new variable to a stratified survey design and calculate the mean.

```R
data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat,
fpc=~fpc)
dstrat<-update(dstrat, apidiff=api00-api99)
svymean(~api99+api00+apidiff, dstrat)
```

--------------------------------

### Usage of make.calfun

Source: http://r-survey.r-forge.r-project.org/survey/html/make.calfun.html

The basic syntax for defining a new calibration function.

```R
make.calfun(Fm1, dF, name)
```

--------------------------------

### Fit and predict with survey-weighted Cox models

Source: http://r-survey.r-forge.r-project.org/survey/html/svycoxph.html

Example demonstrating model fitting with survey designs and generating survival curve predictions.

```R
## Somewhat unrealistic example of nonresponse bias.
data(pbc, package="survival")

pbc$randomized<-with(pbc, !is.na(trt) & trt>0)
biasmodel<-glm(randomized~age*edema,data=pbc,family=binomial)
pbc$randprob<-fitted(biasmodel)
if (is.null(pbc$albumin)) pbc$albumin<-pbc$alb ##pre2.9.0

dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
rpbc<-as.svrepdesign(dpbc)

(model<-svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=dpbc))

svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=rpbc)

s<-predict(model,se=TRUE, type="curve",
     newdata=data.frame(bili=c(3,9), protime=c(10,10), albumin=c(3.5,3.5)))
plot(s[[1]],ci=TRUE,col="sienna")
lines(s[[2]], ci=TRUE,col="royalblue")
quantile(s[[1]], ci=TRUE)
confint(s[[2]], parm=365*(1:5))
```

--------------------------------

### Configure Lonely PSU Options

Source: http://r-survey.r-forge.r-project.org/survey/exmample-lonely.html

Global options for handling lonely PSUs in variance estimation.

```R
 options(survey.lonely.psu="fail")
```

```R
 options(survey.lonely.psu="remove")
```

```R
 options(survey.lonely.psu="certainty")
```

```R
 options(survey.lonely.psu="adjust")
```

--------------------------------

### Apply smoothing to survey data

Source: http://r-survey.r-forge.r-project.org/survey/example-graphics.html

Uses svysmooth to visualize trends in survey data, with examples of subsetting and adding categorical means.

```R
load("ex2.RData")
shs.des <- svydesign(id=~PSU, weights=~GROSSWT,strata=~STRATUM,data=shs)

plot(c(15,83),c(0,65),type='n',xlab='age',ylab='% internet users')
legend("topright",lty=c(1,1),col=c("blue","red"),legend=c("Male","Female"),bty="n")
lf<-svysmooth(I(100*INTUSE)~AGE,subset(shs.des,SEX=="female" & !is.na(AGE)),bandwidth=10)
lines(lf,col="red")
lm<-svysmooth(I(100*INTUSE)~AGE,subset(shs.des,SEX=="male" & !is.na(AGE)),bandwidth=10)
lines(lm, col="blue")

shs.des <- update(shs.des, agegp=cut(AGE, c(0,25,35,45,55,66,75,Inf)))
means<-svyby(~INTUSE,~agegp+SEX,svymean,design=shs.des, na.rm=TRUE)
points((2:8)*10, means[1:7,3]*100,col="blue")
points((2:8)*10, means[7+1:7,3]*100,col="red")
```

--------------------------------

### Fit MLE models with svymle

Source: http://r-survey.r-forge.r-project.org/survey/html/svymle.html

Demonstrates fitting models using svymle with and without user-defined gradients, compared against standard svyglm.

```R
 data(api)

 dstrat<-svydesign(id=~1, strata=~stype, weight=~pw, fpc=~fpc, data=apistrat)

 ## fit with glm
 m0 <- svyglm(api00~api99+ell,family="gaussian",design=dstrat)
 ## fit as mle (without gradient)
 m1 <- svymle(loglike=dnorm,gradient=NULL, design=dstrat, formulas=list(mean=api00~api99+ell, sd=~1),start=list(c(80,1,0),c(20)), log=TRUE)
 ## with gradient
 gr<- function(x,mean,sd,log){
         dm<-2*(x - mean)/(2*sd^2)
         ds<-(x-mean)^2*(2*(2 * sd))/(2*sd^2)^2 - sqrt(2*pi)/(sd*sqrt(2*pi))
         cbind(dm,ds)
      }
 m2 <- svymle(loglike=dnorm,gradient=gr, design=dstrat, formulas=list(mean=api00~api99+ell, sd=~1),
               start=list(c(80,1,0),c(20)), log=TRUE, method="BFGS")

 summary(m0)
 summary(m1,stderr="model")
 summary(m2)

 ## More complicated censored data example
 ## showing that the response variable can be multivariate

 data(pbc, package="survival")
 pbc$randomized <- with(pbc, !is.na(trt) & trt>0)
 biasmodel<-glm(randomized~age*edema,data=pbc)
 pbc$randprob<-fitted(biasmodel)
 dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema,
    data=subset(pbc,randomized))

 lcens<-function(x,mean,sd){
    ifelse(x[,2]==1,
           dnorm(log(x[,1]),mean,sd,log=TRUE),
           pnorm(log(x[,1]),mean,sd,log=TRUE,lower.tail=FALSE)
           )
 }

 gcens<- function(x,mean,sd){

        dz<- -dnorm(log(x[,1]),mean,sd)/pnorm(log(x[,1]),mean,sd,lower.tail=FALSE)

        dm<-ifelse(x[,2]==1,
                   2*(log(x[,1]) - mean)/(2*sd^2),
                   dz*-1/sd)
        ds<-ifelse(x[,2]==1,
                   (log(x[,1])-mean)^2*(2*(2 * sd))/(2*sd^2)^2 - sqrt(2*pi)/(sd*sqrt(2*pi)),
                   ds<- dz*-(log(x[,1])-mean)/(sd*sd))
        cbind(dm,ds)      
 }

if(!is.null(pbc$albumin)){
  svymle(loglike=lcens, gradient=gcens, design=dpbc,
      formulas=list(mean=I(cbind(time,status>0))~bili+protime+albumin,
                    sd=~1),
         start=list(c(10,0,0,0),c(1)))
} else {
  svymle(loglike=lcens, gradient=gcens, design=dpbc,
      formulas=list(mean=I(cbind(time,status>0))~bili+protime+alb,
                    sd=~1),
         start=list(c(10,0,0,0),c(1)))
}
```

--------------------------------

### Initialize Survey Design

Source: http://r-survey.r-forge.r-project.org/survey/example-poststrat.html

Create a survey design object for post-stratification.

```R
 > dpets <- svydesign(id = ~1, weight = ~weight, data = pets, fpc = ~n)
```
Tokens:
Raw