A bivariate hierarchical Bayesian model for estimating cropland cash rental rates at the county level
Section 3. Bivariate hierarchical Bayesian model

• Previous
• Next

The correlation between the 2009 and 2010 cash rental rates observed in Section 2.1.1 suggests that using the information in the data from 2009 has the potential to improve the predictions for 2010. A bivariate hierarchical model for a state is specified as a way to incorporate the data for both years. Let $a_{ij\text{}\mathrm{,}t}$ and $y_{ij\text{}\mathrm{,}t}$ be the acres and dollars per acre, respectively, rented by operator $j$ in county $i$ and year $t\left(t\mathrm{<mo>=</mo>\; 09,}10\right),$ and let $x_{ij\text{}\mathrm{,}t}$ be the associated column vector of auxiliary variables with dimension $p_t.$ For covariates that are constant across years and individuals, $x_{ij\text{}\mathrm{,}t}\mathrm{<mo>=</mo>}{x}_{i109}.$ Let $w_{ij\text{}\mathrm{,}t}\mathrm{<mo>=</mo>}{a}_{ij\text{}\mathrm{,}t}{N}_{g\left(ijt\right)}{n}_{g\left(ijt\right)}^{-1},$ where $N_{g\left(ijt\right)}$ and $n_{g\left(ijt\right)}$ are the population size and number of respondents, respectively, in year $t$ for the stratum $g$ that contains unit $\left(ij\right).$

To specify the model, we divide the respondents into three sets:

• Set 1 consists of units $\left(ij\right)$ that report a non-irrigated cash rental rate in both 2009 and 2010.
• Set 2 consists of units $\left(ij\right)$ that only report a non-irrigated cash rental rate in 2009.
• Set 3 consists of units $\left(ij\right)$ that only report a non-irrigated cash rental rate in 2010.

We assume that observations in set 1 satisfy the bivariate model

$$\left(\begin{array}{c}{y}_{ij\text{}\mathrm{,}09}\\ y_{ij\text{}\mathrm{,}10}\end{array}\right)\mathrm{<mo>=</mo>}\left(\begin{array}{c}{x}_{ij\text{}\mathrm{,}09}^{\prime }{\beta }_{09}+{\nu }_{i\text{}\mathrm{,}09}+e_{ij\text{}\mathrm{,}09}\\ x_{ij\text{}\mathrm{,}10}^{\prime }{\beta }_{10}+{\nu }_{i\text{},10}+e_{ij\text{}\mathrm{,}10}\end{array}\right)\mathrm{,}(3.1)$$

where

$$\left(\begin{array}{c}{e}_{ij\text{}\mathrm{,}09}\\ e_{ij\text{}\mathrm{,}10}\end{array}\right)\sim \text{N}\left(0\mathrm{,}{D}_{wij}^{-\text{0}\text{.5}}{\Sigma }_{ee}{D}_{wij}^{-\text{0}\text{.5}}\right)\mathrm{,}(3.2)$$

$D_{wij}=\text{diag}\left(w_{ij\text{}\mathrm{,}09}\mathrm{,}{w}_{ij\text{}\mathrm{,}10}\right),$ and

$$\left(\begin{array}{c}{\nu }_{i\text{}\mathrm{,}09}\\ {\nu }_{i10}\end{array}\right)\sim \text{N}\left(0\mathrm{,}{\Sigma }_{\nu \nu }\right)\mathrm{.}(3.3)$$

We denote the diagonal elements of ${\Sigma }_{ee}$ corresponding to 2009 and 2010 by ${\sigma }_{ee09}$ and ${\sigma }_{ee10},$ respectively. For units $\left(ij\right)$ in set 2 or 3, we assume

$$y_{ij\text{},t}\mathrm{<mo>=</mo>}{x}_{ij\text{}\mathrm{,}t}^{\prime }{\beta }_t+{\nu }_{i\text{}\mathrm{,}t}+e_{ij\text{}\mathrm{,}t}^{*}\mathrm{,}(3.4)$$

where $e_{ij\text{}\mathrm{,}t}^{*}\sim \text{N}\left(\mathrm{0,}{w}_{ij\text{}\mathrm{,}t}^{-1}{\tau }_{e\text{}\mathrm{,}t}^2\right),$ $t\mathrm{<mo>=</mo>\; 09}$ for set 2, and $t\mathrm{<mo>=</mo>\; 10}$ for set 3. The model not only allows the variances for the unit-level errors to differ across time points but also allows the variances of unit-level errors for units that respond in both time points to differ from the variances for units that only respond in one time-point. The quantity to predict for 2010 is

$${\theta }_{i\text{}\mathrm{,}10}={\overline{x}}_{N_i\text{}\mathrm{,}10}^{\prime }{\beta }_{10}+{\nu }_{i\text{}\mathrm{,}10}\mathrm{,}(3.5)$$

where ${\overline{x}}_{N_i\text{}\mathrm{,}10}$ is the population mean of the covariates for county $i.$

The variances of the unit-level errors, $e_{ij\text{}\mathrm{,}t}$ and $e_{ij\text{}\mathrm{,}t}^{*},$ are assumed to be inversely proportional to the weight, $w_{ij\text{}\mathrm{,}t},$ for two reasons. First, incorporating the weights in the model aims to reduce bias that could arise if the design is informative for the model. As explained in Section 2, the weights depend on the dollar value of the land rented from the previous year. Therefore, the possibility that the sample design may be informative for a model without the weights is plausible. If ${\Sigma }_{ee}$ and ${\Sigma }_{vv}$ are diagonal, and if ${\tau }_{e\text{}\mathrm{,}t}^2\mathrm{<mo>=</mo>}{\sigma }_{eet},$ then in a frequentist framework, an empirical best linear unbiased predictor for the county $i$ mean in year $t$ is the design-consistent pseudo-eblup of You and Rao (2002). The second reason to incorporate the weights is that the variances of residuals from preliminary analyses decrease as the acres increase.

Diffuse, proper priors are specified for the unknown regression coefficients and variances. Specifically, ${\beta }_t\sim \text{N}\left(0\mathrm{,}{10}^{6}I\right),$ and ${\tau }_{e\text{}\mathrm{,}t}^2\sim \text{inverse}-\text{gamma}\left(\text{0}\text{.001}\mathrm{,}\text{0}\text{.001}\right).$ The covariance matrices, ${\Sigma }_{ee}$ and ${\Sigma }_{\nu \nu }$ have inverse-Wishart prior distributions with shape parameter 0.01 and a diagonal scale matrix with diagonal elements 0.001. The parameterizations for the inverse-gamma and inverse-Wishart distributions are from Gelman, Carlin, Stern and Rubin (2009). We choose priors with conjugate forms for computational simplicity. The choices of the hyperparameters are selected to be un-informative relative to the data for the Cash Rents Survey application.

3.1  Gibbs sampling and posteriors

We use Gibbs sampling to obtain a Monte Carlo approximation to the posterior distribution. An analysis of BGR statistics (Gelman et al., 2009) based on three MCMC chains, each with 20,000 iterations, indicated that 1,000 iterations is sufficient for burn-in. The analyses in Section 4 are based on one chain of length 20,000 for each of the three states, Iowa, Kansas and Texas, where the first 1,000 iterations are discarded for burn-in. By the choices of the likelihood and the priors, the full conditional distributions are known distributions. See Appendix A.

3.2  Prediction and MSE estimation

If ${\overline{x}}_{N_i\text{}\mathrm{,}10}$ is known, the Bayes predictor of ${\theta }_{i\text{}\mathrm{,}10}$ for squared error loss is

$${\tilde{\theta }}_{i\text{}\mathrm{,}10}^B\mathrm{<mo>=</mo>}E\left[{\theta }_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\mathrm{,}{\overline{x}}_{N_i\text{}\mathrm{,}10}\right]\mathrm{<mo>=</mo>}{\overline{x}}_{N_i\text{}\mathrm{,}10}^{\prime }{\widehat{\beta }}_{10}+E\left[v_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\right]\mathrm{,}(3.6)$$

where ${\widehat{\beta }}_{10}\mathrm{<mo>=</mo>}E\left[{\beta }_{10}|\left(y\mathrm{,}x\right)\right],\left(y\mathrm{,}x\right)$ denotes the observed cash rental rates and covariates for the two years, and the second equality in (3.6) follows from (3.5) and linearity of expectation. The posterior mean squared error of ${\tilde{\theta }}_{i\text{}\mathrm{,}10}^B$ is

$$E\left[{\left({\tilde{\theta }}_{i\text{}\mathrm{,}10}^B-{\theta }_{i\text{}\mathrm{,}10}\right)}^2|\left(y\mathrm{,}x\right)\mathrm{,}{\overline{x}}_{N_i\text{}\mathrm{,}10}\right]\mathrm{<mo>=</mo>}V\left\{{\theta }_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\mathrm{,}{\overline{x}}_{N_i\text{}\mathrm{,}10}\right\}\mathrm{.}(3.7)$$

As discussed in Section 2, the population mean of the covariates, ${\overline{x}}_{N_i\text{}\mathrm{,}10},$ is not available for unit-level covariates in the Cash Rent Survey application. To define a predictor, we add a model for the covariate mean. See Lohr and Prasad (2003) for an approach that begins with a model specification for the unit level covariates. Partition $x_{ij\text{}\mathrm{,}10}$ into two sub-vectors, $x_{ij\text{}\mathrm{,}10}^{\left(1\right)}$ and $x_{ij\text{}\mathrm{,}10}^{\left(2\right)},$ where $x_{ij\text{}\mathrm{,}10}^{\left(1\right)}$ contains county-level covariates, and $x_{ij\text{}\mathrm{,}10}^{\left(2\right)}$ contains unit-level covariates. Assume ${\overline{x}}_{wi10}|{\overline{x}}_{N_i\text{}\mathrm{,}10}\sim \text{N}\left({\overline{x}}_{N_i\text{}\mathrm{,}10}\mathrm{,}{V}_{xxi\text{}\mathrm{,}10}\right),$ where ${\overline{x}}_{wi10}\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{n_{i10}}}{w}_{ij\text{}\mathrm{,}10}\right)}^{-1}\text{}\left({\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{n_{i10}}}{w}_{ij\text{}\mathrm{,}10}{x}_{ij\text{}\mathrm{,}10}\right),$ $n_{i10}$ is the sum of the number of units in set 1 and in set 3, and $V_{xxi\text{}\mathrm{,}10}$ is known. The elements of $V_{xxi\text{}\mathrm{,}10}$ corresponding to $x_{ij\text{}\mathrm{,}10}^{\left(1\right)}$ are 0, and we explain how we obtain the elements of $V_{xxi\text{}\mathrm{,}10}$ corresponding to unit-level covariates in Appendix B. The Central Limit Theorem supports the assumption of normality for ${\overline{x}}_{wi10}$ even if the distribution of the unit-level covariate values is not normal (Kim, Park and Lee, 2017). Assuming ${\overline{x}}_{N_i\text{}\mathrm{,}10}$ has a flat prior, ${\overline{x}}_{N_i\text{}\mathrm{,}10}|{\overline{x}}_{wi10}\sim \text{N}\left({\overline{x}}_{wi10}\mathrm{,}{V}_{xxi\text{}\mathrm{,}10}\right).$ The Bayes predictor of ${\theta }_{i\text{}\mathrm{,}10}$ for squared error loss under the extended model in which the population mean of the covariates is unknown is

$${\widehat{\theta }}_{i\text{}\mathrm{,}10}^B\mathrm{<mo>=</mo>}{\overline{x}}_{wi10}^{\prime }{\widehat{\beta }}_{10}+E\left[v_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\right]\mathrm{.}(3.8)$$

The posterior mean squared error of ${\widehat{\theta }}_{i\text{}\mathrm{,}10}^B$ is

$$\begin{array}{ll}E\left[{\left({\widehat{\theta }}_{i\text{}\mathrm{,}10}^B-{\theta }_{i\text{}\mathrm{,}10}\right)}^2|\left(y\mathrm{,}x\right)\right]\hfill & \mathrm{<mo>=</mo>}E\left\{{\left({\widehat{\theta }}_{i\text{}\mathrm{,}10}^B-{\tilde{\theta }}_{i\text{}\mathrm{,}10}^B+{\tilde{\theta }}_{i\text{}\mathrm{,}10}^B-{\theta }_{i\text{}\mathrm{,}10}\right)}^2|\left(y\mathrm{,}x\right)\right\}\hfill \\ \hfill & \mathrm{<mo>=</mo>}E\left\{{\left({\widehat{\theta }}_{i\text{}\mathrm{,}10}^B-{\tilde{\theta }}_{i\text{}\mathrm{,}10}^B\right)}^2|\left(y\mathrm{,}x\right)\right\}\hfill \\ \hfill & +2E\left\{E\left[\left({\widehat{\theta }}_{i\text{}\mathrm{,}10}^B-{\tilde{\theta }}_{i\text{}\mathrm{,}10}^B\right)\left({\tilde{\theta }}_{i\text{}\mathrm{,}10}^B-{\theta }_{i\text{}\mathrm{,}10}\right)|\left(y\mathrm{,}x\right)\mathrm{,}{\overline{x}}_{N_i\text{}\mathrm{,}10}\right]|\left(y\mathrm{,}x\right)\right\}\hfill \\ \hfill & +V\left\{{\theta }_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\right\}\hfill \\ \hfill & \mathrm{<mo>=</mo>}{\widehat{\beta }}_{10}^{\prime }V\left\{{\overline{x}}_{N_i\text{}\mathrm{,}10}|{\overline{x}}_{wi10}\right\}{\widehat{\beta }}_{10}+V\left\{{\theta }_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\right\}\hfill \\ \hfill & \mathrm{<mo>=</mo>}{\widehat{\beta }}_{10}^{\prime }V\left\{{\overline{x}}_{N_i\text{}\mathrm{,}10}|{\overline{x}}_{wi10}\right\}{\widehat{\beta }}_{10}\hfill \\ \hfill & +V\left\{{\overline{x}}_{wi10}^{\prime }{\beta }_{10}+v_{i\text{}\mathrm{,}10}+{\left({\overline{x}}_{N_i\text{}\mathrm{,}10}-{\overline{x}}_{wi10}\right)}^{\prime }{\beta }_{10}|\left(y\mathrm{,}x\right)\right\}\hfill \\ \hfill & \approx {\widehat{\beta }}_{10}^{\prime }V\left\{{\overline{x}}_{N_i\text{}\mathrm{,}10}|{\overline{x}}_{wi10}\right\}{\widehat{\beta }}_{10}+V\left\{{\overline{x}}_{wi10}^{\prime }{\beta }_{10}+v_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\right\}\mathrm{,}(3.9)\hfill \end{array}$$

where the final approximation assumes that the $Cov\left\{{\overline{x}}_{wi10}^{\prime }{\beta }_{10}+v_{i\text{}\mathrm{,}10}\mathrm{,}{\left({\overline{x}}_{N_i\text{}\mathrm{,}10}-{\overline{x}}_{wi10}\right)}^{\prime }{\beta }_{10}|\left(y\mathrm{,}x\right)\right\}$ is negligible. A comparison of (3.7) and (3.9) shows that the term ${\widehat{\beta }}_{10}^{\prime }V\left\{{\overline{x}}_{N_i\text{}\mathrm{,}10}|{\overline{x}}_{wi10}\right\}{\widehat{\beta }}_{10}$ accounts for the increase in posterior MSE due to replacing ${\overline{x}}_{N_i\text{}\mathrm{,}10}$ in (3.6) with ${\overline{x}}_{wi10}$ in (3.8). To quantify the posterior MSE of ${\widehat{\theta }}_{i\text{}\mathrm{,}10}^B,$ we use

$$\text{MSE}\left({\widehat{\theta }}_{i\text{}\mathrm{,}10}^B\right)\mathrm{<mo>=</mo>}{\widehat{\text{MSE}}}_{1i}+{\widehat{\text{MSE}}}_{2i}\mathrm{,}(3.10)$$

where ${\widehat{\text{MSE}}}_{1i}\mathrm{<mo>=</mo>}V\left\{{\overline{x}}_{wi10}^{\prime }{\beta }_{10}+{\nu }_{i\text{}\mathrm{,}10}|\left(y\mathrm{,}x\right)\right\},$ and ${\widehat{\text{MSE}}}_{2i}\mathrm{<mo>=</mo>}{\widehat{\beta }}_{10}^{\prime }{V}_{xxi\text{}\mathrm{,}10}{\widehat{\beta }}_{10}$ . In the application of Section 4, we evaluate the effect of including the term ${\widehat{\text{MSE}}}_{2i},$ which accounts for the increase in posterior MSE due to use of the sample mean of the covariate instead of the population mean, on the posterior MSE of the predictor.

3.3  Two-stage benchmarking

NASS obtains estimates of cash rental rates at the state level using data from a national survey conducted in June (the June Area Survey) in addition to the Cash Rent Survey. The state estimates are published before the county-level data from the Cash Rent Survey are fully processed. NASS also establishes estimates of cash rental rates for agricultural statistics districts. To retain internal consistency, appropriately weighted sums of county estimates must equal the district estimates and appropriately weighted sums of district estimates must equal the previously published state estimate. Letting ${\widehat{\theta }}_{i10}$ be the benchmarked predictor for 2010, the benchmarking restrictions for a single time-point are defined by

$${\displaystyle \sum _{i\in d_k}}{w}_{i10}{\widehat{\theta }}_{i10}\mathrm{<mo>=</mo>}{\widehat{\lambda }}_{k10}\mathrm{,}(3.11)$$

and

$${\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}{\eta }_{k10}{\widehat{\lambda }}_{k10}\mathrm{<mo>=</mo>}{\theta }_{\text{pub}10}\mathrm{,}(3.12)$$

where $k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K$ index the districts, $w_{i10}\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{i\in d_k}}{z}_{i10}\right)}^{-1}\text{}{z}_{i10},$

$${\eta }_{k10}\mathrm{<mo>=</mo>}{\left({\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}{\displaystyle \sum _{i\in d_k}}{z}_{i10}\right)}^{-1}\text{}{\displaystyle \sum _{i\in d_k}}{z}_{i10}\mathrm{,}$$

$z_{i10}$ is the direct estimator of the acres rented in county $i$ in year 2010, $d_k$ is the index set for the counties in district $k,$ ${\widehat{\lambda }}_{k10}$ is the final estimate of the average cash rental rate for district $k,$ and ${\theta }_{\text{pub}10}$ is the published estimate of the state-level cash rent per acre. We consider estimates for the year 2010 in (3.11) and (3.12) because we focus on estimation for 2010 in the analysis of Section 4.

We use the two-stage benchmarking procedure proposed by Ghosh and Steorts (2013) to define benchmarked estimates. The benchmarked estimates minimize the quadratic form

$$g\left(c\mathrm{,}b\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}{\displaystyle \sum _{i\in d_k}}{\xi }_i{\left({\widehat{\theta }}_{i10}^B-c_i\right)}^2+{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}{\rho }_k{\left({\widehat{\theta }}_{k\mathrm{10,}w}^B-b_k\right)}^2(3.13)$$

subject to the constraints in (3.11) and (3.12), where $c\mathrm{<mo>=</mo>}\left(c_1\mathrm{,}\dots \mathrm{,}{c}_D\right),$ $D$ denotes the total number of counties, $b\mathrm{<mo>=</mo>}\left(b_1\mathrm{,}\dots \mathrm{,}{b}_K\right),$ ${\widehat{\theta }}_{k\mathrm{10,}w}^B\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in d_k}}{w}_{i10}{\widehat{\theta }}_{i10}^B,$ and $\left({\rho }_k\mathrm{,}{\xi }_i\right)$ are constants selected by the analyst. We set ${\xi }_i\mathrm{<mo>=</mo>}{w}_{i10}$ and ${\rho }_k\mathrm{<mo>=</mo>}{\eta }_{k10},$ which gives the benchmarked estimates

$${\widehat{\theta }}_{i10}\mathrm{<mo>=</mo>}{\widehat{\theta }}_{i10}^B+{\widehat{\lambda }}_{k\left(i\right)10}-{\widehat{\theta }}_{k\left(i\right)\mathrm{10,}w}^B\mathrm{,}(3.14)$$

with

$${\widehat{\lambda }}_{k\left(i\right)10}\mathrm{<mo>=</mo>}{\widehat{\theta }}_{k\left(i\right)\mathrm{10,}w}^B+\frac{\left({\theta }_{\text{pub}10}-{\widehat{\theta }}_{w10}^B\right){\eta }_{k\left(i\right)10}{\left(1+{\eta }_{k\left(i\right)10}\right)}^{-1}}{{\displaystyle {\sum }_{i\in d_{k\left(i\right)10}}{\eta }_{k\left(i\right)10}^2}{\left(1+{\eta }_{k\left(i\right)10}\right)}^{-1}}\mathrm{,}(3.15)$$

for county $i$ and district $k\left(i\right),$ respectively, where $k\left(i\right)$ is the district containing county $i.$ In (3.15), ${\widehat{\theta }}_{w10}^B\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^K}{\eta }_{k10}{\widehat{\theta }}_{k\mathrm{10,}w}^B.$ Each of the benchmarked estimates in (3.14) and (3.15) is a sum of the hierarchical Bayes predictor and an adjustment term. If the hierarchical Bayes predictor for the state is larger (smaller) than the previously published state total, then the adjustment is negative (positive), and the benchmarked county and district estimates are smaller (larger) than the hierarchical Bayes predictors. The posterior mean squared error of the benchmarked predictor for year $t$ is

$${\text{MSE}}_{i10}^{B\text{Bench}}\mathrm{<mo>=</mo>}\text{MSE}\left({\widehat{\theta }}_{i10}^B\right)+{\left({\widehat{\theta }}_{i10}^B-{\widehat{\theta }}_{i10}\right)}^2\mathrm{,}(3.16)$$

where $\text{MSE}\left({\widehat{\theta }}_{i10}^B\right)$ is defined in (3.10). See (You, Rao and Dick, 2004) for a derivation of the posterior MSE of a benchmarked predictor.

• Previous
• Article main page
• Next