R language uses WinBUGS software to build hierarchical (hierarchical) Bayesian model for sat

Time:2021-5-4

Link to the original text:http://tecdat.cn/?p=11974


R2winbugs package provides a convenient function to call WinBUGS from R. It automatically writes data and scripts in WinBUGS readable format for batch processing (since version 1.4). After the WinBUGS process is completed, the result data can be read into R through the package itself (which provides a compact graphical summary of inference and convergence diagnosis), and the output can be further analyzed by using the function of coda package.

  WinBUGS software can be downloaded fromhttp://www.mrc-bsu.cam.ac.uk/…。  

R is a “language for data analysis and graphics processing”, which is an open source and free statistical software package to realize the language. Please refer tohttp://www.R-project.org/.   R and r2winbugs are available from cran, i.ehttp://CRAN.R-Project.org Or its. If you can use an Internet connection, you can install r2winbugs by typing install. Packages (“r2winbugs”) at the R command prompt. Don’t forget to use library (r2winbugs)

example

School data

The academic aptitude test (SAT) measures the ability of high school students to help universities make entrance decisions. Our data came from an experiment conducted in the late 1970s with sat-v (language proficiency test) from eight high schools. Sat-v is a standard multiple choice test administered by the educational testing service. The service is interested in the effectiveness of the coaching program for each school selected.

realization

The implementation of r2winbugs is very simple.   The “main” function bugs() is called by the user. In principle, it is right   The packages of other functions that are called step by step are as follows:

  1. Bugs. Data. Inits() writes the data files’ data. TXT ‘and’ inits1. TXT ‘,’inits2. TXT’… Into the working directory.
  2. Bugs. Script() writes to the file “script. TXT” used by WinBUGS for batch processing.
  3. Bugs. Run() updates the WinBUGS registry, calls WinBUGS, and runs it in batch mode using ‘script. TXT’.
  4. Bugs. Sims() is called only if the parameter codapkg is set to false (the default).
    Otherwise, bugs () returns the file name where the data is stored. For example, these can be imported through packaged coda, which provides convergence diagnosis, Monte Carlo estimation calculation, trace diagram and other functions.
    The bugs. Sims() function reads the simulation in WinBUGS into R, formats it, monitors the convergence, performs the convergence check, and calculates the median and quantile. It also prepares output for bugs () itself.
    These functions are not called directly by the user. Arguments are passed from bugs () to other functions.

example

We apply the functions provided by r2winbugs to the sample data and analyze the output.

School data

Sample data:

> schools 

R language uses WinBUGS software to build hierarchical (hierarchical) Bayesian model for sat


In order to model these data, we use the hierarchical model proposed by Gelman et al. We assume that the estimated values of each school have a normal distribution, and the mean theta and variance tau. Y, and the inverse variance is 1= σ. Y2, whose prior distribution is uniform on (01000). For the mean theta, we use another normal distribution, the mean is mu. Theta and the inverse variance is tau. Theta. For its prior distribution, see the following WinBUGS Code:

model {
for (j in 1:J)
{
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}

The model must be stored in a separate file, such as’ schools. Bug ‘2, in an appropriate directory, such as C: / schools /. In R, the user must prepare the data input required by the bugs() function. This can be a list containing the name of each data vector, for example

> J <- nrow(schools) 

Using these data and model files, we can run MCMC simulation to get the estimates of theta, mu. Theta and sigma. Theta. Before running, the user must determine how many chains (N. chain = 3) and iterations (N. ITER = 1000) to run. In addition, the user must specify the initial value of the chain

> inits <- function(){
+ list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),
+ sigma.theta = runif(1, 0, 100))
+ }

You can start MCMC simulation. The parameter bugs. Directory in R must point to the installation directory of WinBUGS. You can easily output the results in the school. SIM object through print (schools. SIM).  
For this example, you will get similar results

Inference for Bugs model at "c:/schools/schools.bug"
3 chains, each with 1000 iterations (first 500 discarded)
n.sims = 1500 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
theta[1] 11.1 9.1 -3.0 5.0 10.0 16.0 31.8 1.1 39
theta[2] 7.6 6.6 -4.7 3.3 7.8 11.6 21.1 1.1 42
theta[3] 5.7 8.4 -12.5 0.6 6.1 10.8 21.8 1.0 150
theta[4] 7.1 7.0 -6.6 2.7 7.2 11.5 21.0 1.1 42
theta[5] 5.1 6.8 -9.5 0.7 5.2 9.7 18.1 1.0 83
theta[6] 5.7 7.3 -9.7 1.0 6.2 10.2 20.0 1.0 56
theta[7] 10.4 7.3 -2.1 5.3 9.8 15.3 25.5 1.1 27
theta[8] 8.3 8.4 -6.6 2.8 8.1 12.7 26.2 1.0 64
mu.theta 7.6 5.9 -3.0 3.7 8.0 11.0 19.5 1.1 35
sigma.theta 6.7 5.6 0.3 2.8 5.1 9.2 21.2 1.1 46
deviance 60.8 2.5 57.0 59.1 60.2 62.1 66.6 1.0 170
pD = 3 and DIC = 63.8 (using the rule, pD = var(deviance)/2)
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC is an estimate of expected predictive error (lower deviance is better).

In addition, the user can generate the result graph by inputting plot (schools. SIM). The results are shown in the figure. In the figure, the left column shows a quick summary of the following:
Inference and convergence (RB of all parameters are close to 1.0, which indicates that the three chains are well mixed, so they converge approximately); The right column shows the inference for each group of parameters. As you can see from the right column, r2winbugs uses the parameter names in WinBUGS to construct the output as scalar, vector and parameter array.
R language uses WinBUGS software to build hierarchical (hierarchical) Bayesian model for sat