This is the second of a two-course sequence introducing the fundamentals of Bayesian statistics. It builds on the course Bayesian Statistics: From Concept to Data Analysis, which introduces Bayesian methods through use of simple conjugate models. Real-world data often require more sophisticated models to reach realistic conclusions. This course aims to expand our “Bayesian toolbox” with more general models, and computational techniques to fit them. In particular, we will introduce Markov chain Monte Carlo (MCMC) methods, which allow sampling from posterior distributions that have no analytical solution. We will use the open-source, freely available software R (some experience is assumed, e.g., completing the previous course in R) and JAGS (no experience required). We will learn how to construct, fit, assess, and compare Bayesian statistical models to answer scientific questions involving continuous, binary, and count data. This course combines lecture videos, computer demonstrations, readings, exercises, and discussion boards to create an active learning experience. The lectures provide some of the basic mathematical development, explanations of the statistical modeling process, and a few basic modeling techniques commonly used by statisticians. Computer demonstrations provide concrete, practical walkthroughs. Completion of this course will give you access to a wide range of Bayesian analytical tools, customizable to your data.
JAGS model (linear regression)
Loading...
Del curso dictado por University of California, Santa Cruz
Bayesian Statistics: Techniques and Models
141 calificaciones
De la lección
Common statistical models
Linear regression, ANOVA, logistic regression, multiple factor ANOVA
Conoce a los instructores
Matthew HeinerDoctoral Student
Applied Mathematics and Statistics
Let's recreate this model now using JAGS.
First, we'll need to load the rjags library.
And we'll specify the hierarchical form of the model, let's call this mod1.
There will be several different models.
So let's number this one and we're going to start with the model string.
The syntax for this in JAGS is to say model with curly braces and
fill in the different pieces.
Let's start with the likelihood for i in 1:n the length of the data.
Observation y[i] will come from
a normal distribution dnorm.
Where the mean is mu[i], we'll specify what that is in a moment.
And precision, we'll call it P-R-E-C for now.
Now we need to add the linear model part of this.
mu[i] is linear,
where we have coefficient 1 plus
coefficient 2 times the log income,
the ith value of the log income.
So, y[i] comes from this normal distribution with mean mu[i], and
the mean comes from our linear model of log_income.
We now need to give a prior for these coefficients.
We'll start another loop for j in 1:2, because they're two coefficients.
Beta j will come from a normal distribution, dnorm
With let's say mean 0 and a variance 1 million.
We'll make it pretty large so that this prior is close to
the non-informative prior model we fit earlier.
Remember, we have to enter it as a precision for the normal distribution.
So this will be 1 divided by the variance of 1.0 times 10 to the 6th power.
We'll use an inverse gamma distribution on the variants of the normal
distribution for the observations.
An inverse gamma on the variance implies a gamma prior
on the precision which is the reciprocal of the variance.
So let's put the prior directly on the precision.
Precision will come from a gamma distribution.
And remember, it's very important to check the JAGS documentation,
the JAGS manual, under the chapter on distributions.
To remind ourselves how JAGS parameterizes these distributions.
In this case, it is the same as with R.
The first parameter for the gamma distribution is a shape.
And the second parameter is a rate.
Let's use the same trick we did before, where we parameterized
the prior in terms of a prior sample size and a prior guess.
We'll use a prior sample size of 5.
And recall we have to divide it by 2 to turn it into the gamma parameter.
And then we have to do the prior sample size again
times our prior guess for the variance,
let's say it's 10, and divide that also by 2.
A gamma prior on the precision using this shape and
this rate parameter is equivalent.
To an inverse gamma prior on the variance with this shape parameter and
this as the scale parameter.
Instead of precision we will probably be interested in the variance,
which is equal to 1 over the precision.
Notice we didn't assign a distribution here
because we gave a prior to the precision.
That implies a prior for
sigma squared, the inverse gamma we've been talking about.
So instead of a distribution here, we give it an equals.
This is a deterministic relationship between the precision, and the variance.
Usually when we're modeling, the most interpretable quantity is the standard
deviation instead of the variance.
So let's monitor the standard deviation,
which is just the square root of the variance.
This completes the model where this section tells us
the likelihood of the data.
This section gives us the priors for the coefficients, and
down here we have the prior for the left over, or the unaccounted for variants.
JAGS requires that we put this in a string.
So I'll put these between quotation marks.
And now we need to set up the model.
First, let's set the random seed.
And we'll tell it what the data are.
So data for model 1 that'll go into JAGS needs
to be a list where we named our variable y in the model.
So we have to give it y down here.
y=dat, we're using our modified
dataset where we removed missing values,
dat$loginfant as our response variable.
n is the number of rows in the dat dataset.
And, we also created a variable here called log_income.
So we need to give it the exact same name down here, log_income.
We'll get dat$logincome.
That's how we specified it in the data set.
The parameters that we want to monitor for
model one are the two coefficients,
beta one and beta two.
And since they were specified and indexed together as a vector I can just call it b.
We have our choice if we want to monitor the precision,
the variance, or the standard deviation from the likelihood.
Let's just look at the standard deviation.
We can give the model initial values, inits1 and
this goes in as a function that creates a list.
We'll give an initial value for the beta coefficients.
There are two of them.
So, let's draw two random normals.
So drawing two random normals, let's say, 0,
and variance or standard deviation, 100.
For initial value of the precision, let's draw from a gamma distribution.
We only need one draw.
And let's just give it, Shape 1 and rate 1.
That's actually an exponential distribution.
So now that we've specified the data, told it which parameters we want to monitor and
given initial values for the different parameters.
We can specify the model itself,
jags.model where we create
a textConnection with mod1 string.
mod1_string that we just specified.
The data is data1_jags
that we just created,
data1_jags.
We'll use the initial values from our function inits1.
And let's say we want to run three different chains.
That's another option to write into this jags.model function.
We want it to run three separate chains.
In each chain it'll initialize the chain with these random draws.
So we'll have different starting values for each chain.
We'll initialize the model by running all of this code.
Initialization was essentially instant.
And now that the model is initialized,
let's give it a little bit of a burn in period by updating
the model for, let's say, 1,000 iterations.
So it ran the model for 1,000 iterations but it didn't keep the samples.
Once we have the burn in period run, let's create our actual
posterior simulation that we're going to keep for model inference.
Mod1_sim is what we'll name it.
It will come from coda.samples.
We're going to create coda.samples from the model.
Model 1, where the variable names
are the parameters that we created.
And we'll run this for let's say 5,000 iterations.
We'll run this line to run the model.
And since we ran three different chains,
sometimes it's useful to combine them into one chain.
So we'll call that mod1_csim for combined simulation.
And we're going to do this by stacking the matrices that
contain the simulations themselves.
Inside mod1.sim the simulations are stored as matrices.
To do this we need to use the do.call function in R and rbind.
Remember we used cbind earlier to combine matrices along the columns.
Now we're going to stack the matrices vertically by combining them on
the rows using mod1_sim.
And run that line.
Explora nuestro catálogo
Inscríbete de manera gratuita y obtén recomendaciones personalizadas, actualizaciones y ofertas.
Coursera brinda acceso universal a la mejor educación del mundo, al asociarse con las mejores universidades y organizaciones, para ofrecer cursos en línea.
© 2018 Coursera Inc. Todos los derechos reservados.
