计算机代考程序代写 chain algorithm Gibbs sampler for missing data: random rounding – cscodehelp代写

Gibbs sampler for missing data: random rounding

Gibbs sampler for missing data: random rounding

October 2021

Background

Statistical agencies are concerned to protect the confidentiality of respondents’ data. One technique used by Statistics New Zealand is to randomly round cell counts in tables that are either published or specially requested by users. The random rounding algorithm (RR3) is as follows:

If a count is a multiple of 3 it is left unchanged;
If a count is not a multiple of 3 it rounded to the nearest multiple of 3 with probability 2/3 and to the second nearest multiple of 3 with probability 1/3

Thus an observed count of 7 is rounded to 6 with probability 2/3 and to 4 with probability 1/3. Observed counts of (0,3,6,…) are left unchanged.

Given an observed count (Y) the random rounding induces a distribution on the non-negative integers that are multiples of three, but this distribution has positive probability mass on, at most, two of these multiples of 3. An unusual distribution.

Since a user of tabular data produced by Stats NZ will see only randomly rounded counts, inference given these counts can be viewed as a missing data problem where the missing data are the true counts. We can use a Gibbs sampler to solve this inference problem by alternately simulating the true counts, given the randomly rounded counts, and simulating the conditional posterior of underlying model parameters given the true counts.

Setting up the model

Let (mathbf{Y}=(Y_{1},Y_{2},ldots,Y_{k})) denote a set of counts and (mathbf{R}) the corresponding set of randomly rounded counts; e.g we might observe (mathbf{R}=(0,0,3,6,3,0,6,3,0,6,0).)

Suppose that our model for the underlying counts is a simple Poisson model with a common expected value for each cell.

[
egin{aligned}
Y_{i}| heta & stackrel{mathrm{indep}}{sim} mathrm{Poisson}( heta);, , i=1,ldots,k \
heta & sim mathrm{Gamma}(alpha,eta)
end{aligned}
] However, because we observe only the randomly rounded counts, our full model for the observed data is [
egin{aligned}
R_{i}|Y_{i}, heta & stackrel{mathrm{indep}}{sim} RR3(Y_{i}); , i=1,ldots,k\
Y_{i}| heta & stackrel{mathrm{indep}}{sim} mathrm{Poisson}( heta);, , i=1,ldots,k \
heta & sim mathrm{Gamma}(alpha,eta)
end{aligned}
] where (RR3) is the distribution induced by the random rounding algorithm. Note the (RR3) distribution does not depend on ( heta) as indicated by the fact that ( heta) does not appear on the right hand side of the first line in the equation above.

Some useful points to note are:

If (R=0) the only possible (Y) values are ((0,1,2))
If (R=r >0) the only possible (Y) values are ((r-2,r-1,r,r+1,r+2))

Moreover given an observed randomly rounded count we can easily evaluate the probabilities for each possible (Y) value, by Bayes theorem.

[
egin{aligned}
Pr(Y=y|R=r, heta) & =frac{Pr(R=r|Y=y)Pr(Y=y| heta)}
{sum_{a}Pr(R=r|Y=a)Pr(Y=a| heta)} \
& = frac{Pr(R=r|Y=y)mathrm{Poisson}(y| heta)}
{sum_{a}Pr(R=r|Y=a)mathrm{Poisson}(a| heta)}
end{aligned}
] The summation in the denominator above is over the very small set of (Y) values (at most 5) that are possible given the observed (R) value. Given our Poisson model for (Y) we can directly calculate the posterior (to R) probabilities for each possible value of (Y). It is therefore straightforward to simulate (Y) values from the posterior distribution given (R). For example we could just use the sample() function to sample from the possible Y values with probabilities for each possible value computed as above. In fact we really need to compute only the numerator value – unnormalised posterior, as sample() samples with probabilities proportional to the specified values.

Gibbs- sampler

Our main interest is the posterior for ( heta,) (p( heta|mathbf{R}).) However since ( heta) is the parameter of the model for the cell counts (mathbf{Y},) we “extend the conversation to include (mathbf{Y})’’ by noting [p( heta|mathbf{R}) = int p( heta, mathbf{Y}|mathbf{R}) mathrm{d}mathbf{Y} ] In practice we therefore need to obtain (p( heta,mathbf{Y}|mathbf{R}).) A Gibbs sampler for approximating this joint posterior alternates between drawing from

(p( heta|mathbf{Y},mathbf{R})=p( heta|mathbf{Y}))
(p(mathbf{Y}| heta,mathbf{R}))

Under our model (1) is straightforward since by the usual conjugacy results this amounts to simulating from a (mathrm{Gamma}(alpha + sum_{i}Y_{i},eta+k)) distribution. Step (2) of the Gibbs sampler can be implemented by direct simulation using the probabilities (Pr(Y=y|R=r)) as outlined above.

Updating functions

The Gibbs sampler alternates by drawing from (1) (p( heta|mathbf{Y})) and (2) and drawing imputed (mathbf{Y}) values from (p(mathbf{Y}|mathbf{R}, heta).)

First we define a function for updating ( heta) assuming a gamma prior with prior parameters (alpha) and (eta.) The function is just a call to rgamma

update_theta <-function(Y,alpha,beta) { theta <- rgamma(n=1,shape=alpha+sum(Y),rate=beta+length(Y)) return(theta) } Next we define functions to simulate the (mathbf{Y}) values given (mathbf{R}) and ( heta.) First we define a function to compute (Pr(R=r|Y=y)) for all (y) values compatible with and observed (R) value RRprobs <- function(R) { if (R > 0) {
y <- seq(from=R-2,to=R +2) pRy <- c(1/3,2/3,1,2/3,1/3) } else { #R=0 y <- c(0,1,2) pRy <- c(1,2/3,1/3) } return(cbind(y,pRy)) } test3 <- RRprobs(3) test3 ## y pRy ## [1,] 1 0.3333333 ## [2,] 2 0.6666667 ## [3,] 3 1.0000000 ## [4,] 4 0.6666667 ## [5,] 5 0.3333333 test0 <- RRprobs(0) test0 ## y pRy ## [1,] 0 1.0000000 ## [2,] 1 0.6666667 ## [3,] 2 0.3333333 We can now define the function to simulate the (mathbf{Y}) values - this will call RRprobs(). Yimpute_poisson <- function(R,theta){ ##vector of R randomly rounded counts ##Yimpute is the ##count theta is the assumed common success probability among ##these trials. ##Could try and vectorize RRprobs but we will just loop over the ##elements of R Y <- vector(length=length(R),mode="integer") for (i in 1:length(R)) { ##obtain possible Y values and RR3 probabilities - the likelihood for this case Yprobs <- RRprobs(R[i]) Yposs <- Yprobs[,1] #possible Y values pRy <- Yprobs[,2] # likelihood for R at each possible value of Y ##get Poisson probabilities for the possible Y values priory <- dpois(Yposs,lambda=theta) likbyprior <- pRy*priory ##unnormalized posterior pmf we could normalize but ##the sample function is happy with unnormalized pmfs Y[i] <- sample(Yposs,size=1,prob=likbyprior) } return(Y) } ##test the functions testR <- c(0,3,6) Yimpute_poisson(testR,theta=1) ## [1] 0 3 4 Y <- Yimpute_poisson(R=testR,theta=1) Y ## [1] 0 1 6 testtheta <- update_theta(Y=Y,alpha=1,beta=1) testtheta ## [1] 1.043106 Gibbs sampler ####Set-up ##the observed data R <- c(0,0,3,6,3,0,6,3,0,6,0,3) ##These are randomly rounded ##counts from 12 cells; nchains <- 5 nsim <- 2000 burnin <- 1000 ##Structures for storing simulated values theta_store <- matrix(NA,nrow=nsim,ncol=nchains) Ystore <- array(NA,dim=c(length(R),nsim,nchains)) ##specify parameters of gamma prior for theta ##Question specifies alpha=1, beta=1 alpha <- 0.01 beta <- 0.01 ##Check what this looks like testvalues <- seq(from=0,to=10,by=0.01) gamma_prior <- dgamma(testvalues,shape=alpha,rate=beta) plot(testvalues,gamma_prior,type="l") Run the Gibbs sampler for (j in 1:nchains) { ##draw an initial value for theta. We draw from an over-dispersed #approximation to the ##posterior. We base this approximation on the posterior that would #be obtained if the R values #were the actual Y values but we use only half the actual number of values ##so we use the first 6 R values as though they were Y values newtheta <- rgamma(n=1,shape=alpha+sum(R[1:6]), rate=beta+length(R[1:6])) for (i in 1:nsim) { # update Y newY <- Yimpute_poisson(R=R,theta=newtheta) newtheta <- update_theta(Y=newY,alpha=alpha,beta=beta) Ystore[,i,j] <- newY theta_store[i,j] <- newtheta } ##end loop over simulations } #end loop over chains Check convergence Define a function to do convergence and effective sample size checking using the Gelman-Rubin split chains approach GRchunk <- function(post){ require(coda) possize <- nrow(post) n1 <- round(possize/2) ##size of first chunk chunk1 <- post[1:n1,] chunk2 <- post[(round(possize/2)+1):possize,] post_chunked <- cbind(chunk1,chunk2) ##obtain the chain means chain_mean <- colMeans(post) ##obtain the within chain variance chain_sd <- apply(post,MARGIN=2,FUN=sd) chain_var <- chain_sd^2 W = mean(chain_var) ##get between chain variance possize_chunked <- nrow(post_chunked) B <- possize_chunked*(sd(chain_mean))^2 ##Now build the components of the Gelman-Rubin statistic Vplus <- ((possize_chunked-1)/possize_chunked) * W + (1/possize_chunked)*B Rhat <- sqrt(Vplus/W) Rhat ##compute estimated effective sample size, using Coda post.mcmc <- coda::as.mcmc(post_chunked ) codaNeff_chains <- coda::effectiveSize(post.mcmc) codaNeff <- sum(codaNeff_chains) outlist <- list(Rhat=Rhat,codaNeff=codaNeff,codaNeff_chains=codaNeff_chains) return(outlist) } Check convergence for ( heta) postheta <- theta_store[(burnin+1):nsim,] posY <- Ystore[,(burnin+1):nsim,] possize <- nsim - burnin GRtheta <- GRchunk(postheta) ## Loading required package: coda GRtheta ## $Rhat ## [1] 0.9994041 ## ## $codaNeff ## [1] 2783.266 ## ## $codaNeff_chains ## var1 var2 var3 var4 var5 var6 var7 var8 ## 283.6742 221.2449 300.5721 257.4455 280.8503 287.9361 297.1881 297.9748 ## var9 var10 ## 285.9517 270.4280 Check convergence of Y’s str(posY) ## num [1:12, 1:1000, 1:5] 1 1 1 6 2 1 7 2 1 4 ... #We will loop over the posterior samples for each of 12 cells extracting the $hat{R}$ and Neff values. #Set-up structures for storing the values. RhatY <- vector(length=length(R),mode="numeric") NeffY <- vector(length=length(R),mode="numeric") for (i in 1:length(R)) { postY <- posY[i,,] GRpostY <- GRchunk(postY) RhatY[i] <- GRpostY$ Y[i] <- GRpostY$codaNeff } RhatY ## [1] 0.9993558 0.9992224 1.0001400 0.9991890 0.9999582 1.0001919 0.9998094 ## [8] 0.9991818 0.9991458 0.9996396 0.9992984 0.9998293 NeffY ## [1] 4838.190 4899.516 4931.166 4702.286 4746.767 4996.341 4829.146 4780.166 ## [9] 4672.349 4560.586 4585.306 4726.874 Produce posterior summaries It looks like the sampler has converged so produce some posterior summaries plot(density(postheta),main="Posterior for theta",lwd=2) ###Compare with a density plot for a simulation of the same size for ###a naive analysis when the ###randomly rounded counts are treated as actual counts plot(density(rgamma(n=5000,shape=alpha+sum(R),rate=beta+length(R))), main="naive posterior for theta",lwd=2) summary((postheta)) ##get summaries by chain ## V1 V2 V3 V4 ## Min. :1.072 Min. :1.023 Min. :1.213 Min. :1.126 ## 1st Qu.:2.147 1st Qu.:2.201 1st Qu.:2.207 1st Qu.:2.184 ## Median :2.528 Median :2.575 Median :2.567 Median :2.553 ## Mean :2.564 Mean :2.600 Mean :2.595 Mean :2.575 ## 3rd Qu.:2.917 3rd Qu.:2.954 3rd Qu.:2.922 3rd Qu.:2.901 ## Max. :4.863 Max. :4.468 Max. :4.493 Max. :5.055 ## V5 ## Min. :1.160 ## 1st Qu.:2.196 ## Median :2.558 ## Mean :2.597 ## 3rd Qu.:2.973 ## Max. :4.782 ##Obtain desired quantiles for 90% credible intervals and median qprobs <- c(0.05,0.5,0.95) quant_theta <- quantile(postheta,probs=qprobs) quant_theta ## 5% 50% 95% ## 1.740068 2.557226 3.558667 ##What would the quantiles be for the naive analysis treating the R values as though they were ###the actual counts. There is no justification for such an analysis but it is the way randomly ###rounded counts are usually handled. quant_naive <- qgamma(qprobs,shape=alpha+sum(R),rate=beta+length(R)) quant_naive ## [1] 1.798707 2.471052 3.293296 quant_theta ## 5% 50% 95% ## 1.740068 2.557226 3.558667 Not a huge difference here, but our interval is wider than the naive interval because we correctly account for the extra random-ness induced by the random-rounding. The approach could clearly be extended to bigger tables of counts and more complex models.