Nothing
## ----echo=FALSE-----------------------------------------------------------------------------------
options(width=100) # width of output
## ----message=FALSE--------------------------------------------------------------------------------
library(survey)
data(api)
# define the regression model
model <- api00 ~ ell + meals + stype + hsg + col.grad + grad.sch
# compute corresponding population totals
XpopT <- colSums(model.matrix(model, apipop))
N <- XpopT[["(Intercept)"]] # population size
# create the survey design object
des <- svydesign(ids=~1, data=apisrs, weights=~pw, fpc=~fpc)
# compute the calibration or GREG estimator
cal <- calibrate(des, formula=model, population=XpopT)
svymean(~ api00, des) # equally weighted estimate
svymean(~ api00, cal) # GREG estimate
## -------------------------------------------------------------------------------------------------
mean(apipop$api00)
## ----message=FALSE--------------------------------------------------------------------------------
library(mcmcsae)
set.seed(1)
sampler <- create_sampler(model, data=apisrs)
sim <- MCMCsim(sampler, verbose=FALSE)
(summ <- summary(sim))
compute_DIC(sim)
## -------------------------------------------------------------------------------------------------
m <- match(apisrs$cds, apipop$cds) # population units in the sample
# use only a sample of 250 draws from each chain
predictions <- predict(sim, newdata=apipop[-m, ], iters=sample(1:1000, 250), show.progress=FALSE)
str(predictions)
samplesum <- sum(apisrs$api00)
summary(transform_dc(predictions, fun = function(x) (samplesum + sum(x))/N))
## -------------------------------------------------------------------------------------------------
summary(predict(sim, newdata=apipop[-m, ], fun=function(x, p) (samplesum + sum(x))/N,
show.progress=FALSE)
)
## -------------------------------------------------------------------------------------------------
n <- nrow(apisrs)
XsamT <- colSums(model.matrix(model, apisrs))
XpopR <- matrix(XpopT - XsamT, nrow=1)
predictions <- predict(sim, X=list(reg1=XpopR), var=N-n, fun=function(x, p) (samplesum + x)/N,
show.progress=FALSE)
summary(predictions)
## ----fig.width=4, fig.height=4, fig.align="center"------------------------------------------------
sampler <- create_sampler(model, data=apisrs,
linpred=list(reg1=matrix(XpopT/N, nrow=1)),
compute.weights=TRUE)
sim <- MCMCsim(sampler, verbose=FALSE)
plot(weights(cal)/N, weights(sim)); abline(0, 1)
sum(weights(sim) * apisrs$api00)
print(summary(sim, "linpred_"), digits=6)
## ----fig.width=4, fig.height=4, fig.align="center", message=FALSE---------------------------------
sampler <- create_sampler(model, formula.V=~vfac(prior=pr_invchisq(df="modeled")),
linpred=list(reg1=matrix(XpopR/N, nrow=1)),
data=apisrs, compute.weights=TRUE)
sim <- MCMCsim(sampler, burnin=1000, n.iter=5000, thin=2, verbose=FALSE)
(summ <- summary(sim))
plot(sim, "vfac1_df")
acceptance_rates(sim)
compute_DIC(sim)
predictions <- predict(sim, newdata=apipop[-m, ], show.progress=FALSE,
fun=function(x, p) (samplesum + sum(x))/N)
summary(predictions)
plot(weights(cal)/N, weights(sim)); abline(0, 1)
summary(get_means(sim, "Q_")[["Q_"]])
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.