Description
1. Load the file munichrents.RData and look at head(rents). The variables are
• RentPerM2: Monthly rent per square meter in Euros
• Year: Year the apartment was built
• Location: District index
• NoHotWater: Indicator variable for no hot water
• NoCentralHeat: Indicator variable for no central heat
• NoBathTiles: Indicator variable for no tiles in the bath
• SpecialBathroom: Indicator variable for a special bathroom
• SpecialKitchen: Indicator variable for a special kitchen
• Room1-Room6: Indicator variables for corresponding number of rooms
Fit a linear model relating rent per square meter to the covariates using least squares,
and extract the coefficient estimates. You can ignore the Location variable for now
since we will later treat this as a spatial random effect. Also note that the room
indicator variables include one that is redundant, so treat a single room as the baseline
(i.e., leave Room1 out of the model, so the intercept corresponds to a single room and
coefficients for the others represent adjustments to the intercept for a different number
of rooms).
2. There are two SpatialPolygons objects associated with this dataset, districts.sp
and parks.sp. The first corresponds to city districts in which apartments may be
located. The second corresponds to districts with no possible apartments, such as
parks or fields.
Create an nb object with neighbors for the districts, defining neighbors as districts
that share a common boundary. Make a plot showing the districts, then add the parks
shaded a different color.
1
3. There are 380 districts in districts.sp, and the corresponding district numbers are
indicated by the Location variable in rents.
I’ve included a matrix H that provides a mapping between the districts as they’re
ordered in districts.sp and as they appear in the rents dataframe. Use H to create
a new vector containing the number of observations in each district, and make a color
or grayscale plot to illustrate this.
4. We will now create a Gibbs sampler to sample from the posterior distribution under the
following Bayesian model. Let X be the matrix of covariates, including the intercept
term. Let n be the number of data points in Y and m be the number of spatial locations
in η.
Data model:
Y |β, η, σ2 ∼ MV N(Xβ + Hη, σ2
I)
Process model:
p(η|τ
2
) ∝ (τ
2
)
−(m−1)/2
exp n
−
1
2τ
2
η
0
(Dw − W)η
o
where W is the matrix of 0 and 1 indicating the neighborhood structure from problem
2, and Dw is a diagonal matrix with diagonal entries P
j W1j
, · · · ,
P
j Wnj . That is η
follows an (improper) intrinsic autoregressive model.
Prior model: Specify independent priors for β, σ2
, and τ
2 with
p(β) ∝ 1, σ2
, τ 2 ∼ InverseGamma(0.001, 0.001)
The full conditional distributions for β, η, σ2
, and τ
2 are given at the end of this
assignment. Construct a Gibbs sampler that cycles through each of the full conditionals
and stores the results for B = 10, 000 iterations. The full conditionals are given below.
A few notes to keep in mind when constructing the sampler:
• The matrix W can be computed from your nb object in problem 2; see help(nb2mat).
I also included objects X and y with the data file.
• The function rinvgamma is in the library MCMCpack.
• IMPORTANT: The intrinsic autoregressive model is an example of pairwise difference prior. It defines proper distributions for the differences ηi − ηj
, but it
implicitly contains a distribution for 1
m
Pm
i=1 ηi that has infinite variance. In
practice, since there is also an intercept term in Xβ, we impose the constraint
Pm
i=1 ηi = 0 when we sample from the full conditional for η. Do this numerically
by subtracting the mean 1
m
Pm
i=1 η
(j)
i
from η
(j)
in each iteration j.
Turn in the following:
• Your map with the neighbors from problem 2
• Your map of the apartment counts for each district
2
• Trace plots and ACF plots for σ
2 and τ
2
• A table with posterior means of the β and 95% credible intervals constructed
using the 0.025 and 0.975 quantiles of the posterior samples
• A color or grayscale map of the posterior means for the vector η
• A color or grayscale map of the posterior standard deviations for the vector η
Full conditionals:
β|Rest ∼ MV N((X
0X)
−1X
0
(Y − Hη), σ2
(X
0X)
−1
)
η|Rest ∼ MV N([H
0H/σ2+(Dw−W)/τ 2
]
−1H
0
(Y −Xβ)/σ2
, [H
0H/σ2+(Dw−W)/τ 2
]
−1
)
σ
2
|Rest ∼ InverseGamma(0.001 + n/2, 0.001 + (Y − Xβ − Hη)
0
(Y − Xβ − Hη)/2)
τ
2
|Rest ∼ InverseGamma(0.001 + (m − 1)/2, 0.001 + η
0
(Dw − W)η/2)
3
