<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/gelmans-mrp-in-r-what-is-this-all-about/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://hagutierrezro.blogspot.com/2017/01/gelman-mrp-in-r-what-is-this-all-about.html"> Data Literacy - The blog of Andrés Gutiérrez</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Multilevel regression with poststratification (MrP) is a useful technique to predict a parameter of interest within small domains through modeling the mean of the variable of interest conditional on poststratification counts. This method (or methods) was first proposed by Gelman and Little (1997) and is widely used in political science where the voting intention is modeling conditional on the interaction of classification variables.</p>
<p>The aim fo this methodology is to provide reliable estimates on strata based on census counts. For those who have some background on survey sampling, this method should look very similar to the Raking method, where sampling weights are adjusted due to known census cell counts. However, a significant difference with Raking is that MrP is a model-based approach, rather than a design-based method. This way, even in the presence of a (maybe complex) survey design, MrP does not take it into account for inference. In other words, sampling design will be considered as ignorable. So, the probability measure that governs the whole inference is based on modeling the voting intention (variable of interest) to demographic categories (auxiliary variables).</p>
<p>Is this a major technical issue to ignore the complex survey design? Yes, because in any case we are considering a probability measure to draw the sample. However, when it comes to voting intention (a major area where this technique is used), we rarely find a sophisticated complex design. Moreover, this kind of studies is barely based on probabilistic polls. So, if the survey lacks a proper sampling plan, it is always better to model the response variable.</p>
<p>Therefore, the ultimate goal of this technique it to estimate a parameter of interest (totals, means, proportions, etc.) for all of the strata (domains, categories or subgroups) in a finite population. From now on, let’s assume that:</p>
<ol>
<li>a population is divided into $H$ strata of interest (for example states),</li>
<li>the parameters of interest are the means (same rules apply for proportions) in each strata $\theta_h$ ($h=1, \ldots, H$), </li>
<li>every stratum is cross-classified by some demographics $j \in H$ (from now on defined as post-srata), besides every population count $N_j$ is known, and</li>
<li>all of the population means $\mu_j$ can be estimated by using some statistical technique, such as multilevel regression.</li>
</ol>
<p>This way, the mean in stratum $h$ ($\theta_h$), is defined as a function of means in post-strata $j$ ($\mu_j$), and post-strata counts ($N_j$):</p>
<p>$$\theta_h = \frac{\sum_{j \in h} N_j \mu_j }{\sum_{j \in h} N_j}$$</p>
<p>The first part of MrP is defined by a multilevel regression (MR). This kind of models is a particular case of mixed effect models. The core components of this phase are</p>
<ul>
<li>the variable of interest (for example, voting intention), </li>
<li>the auxiliary variables (classification on census demographic cells) and, </li>
<li>the random effects (strata of interest, that usually are states or counties).</li>
</ul>
<p>The second part of MrP is cell poststratification (P), and the predicted response variable is aggregated into states and adjusted by corresponding poststratification weights.</p>
<h2>R workshop</h2>
<p>Let’s consider the Lucy database of the TeachingSampling package. This population contains some economic variables of 2396 industrial companies in a particular year. Assume that we want to estimate the mean income of industries by each of the five existing zones (strata of interest) on that database. This way, our parameters of interest are $\theta_1, \ldots, \theta_5$.</p>
<p>Now, we also know that the population is divided into three levels (small, medium and big industries) and we have access to the total number of industries within each cross group. That is, we know exactly how many small industries are on each of the five zones, and how many medium industries are on each of the five zones, and so on. </p>
<p>The following code shows how to load the database and obtain the cell counts. </p>
<pre><span>></span> rm(<span>list</span> <span>=</span> ls())<br><span>></span> set.seed(<span>123</span>)<br><span>></span><br><span>></span> library(TeachingSampling)<br><span>></span> library(dplyr)<br><span>></span> library(lme4)<br><span>></span><br><span>></span> data(<span>&quot;Lucy&quot;</span>)<br><span>></span> <span># Number of industries per level</span><br><span>></span> table(Lucy<span>$</span>Level)<br><br>   Big Medium  Small<br>    <span>83</span>    <span>737</span>   <span>1576</span><br><span>></span> <span># Number of industries per zone</span><br><span>></span> table(Lucy<span>$</span>Zone)<br><br>  A   B   C   D   E<br><span>307</span> <span>727</span> <span>974</span> <span>223</span> <span>165</span><br><span>></span> <span># Size of post-strata</span><br><span>></span> (Np <span>&lt;-</span> table(Lucy<span>$</span>Level, Lucy<span>$</span>Zone))<br>        <br>           A   B   C   D   E<br>  Big     <span>30</span>  <span>13</span>   <span>1</span>  <span>16</span>  <span>23</span><br>  Medium <span>180</span> <span>121</span> <span>111</span> <span>187</span> <span>138</span><br>  Small   <span>97</span> <span>593</span> <span>862</span>  <span>20</span>   <span>4</span></pre>
<p>Of course, this technique works over a selected sample. That’s why we are going to select a random sample of size $n = 1000$. We can also create some tables showing the counts on the sample.</p>
<pre><span>></span> <span># A sample is selected</span><br><span>></span> SLucy <span>&lt;-</span> sample_n(Lucy, <span>size</span> <span>=</span> <span>1000</span>)<br><span>></span> table(SLucy<span>$</span>Level)<br><br>   Big Medium  Small<br>    <span>33</span>    <span>280</span>    <span>687</span><br><span>></span> table(SLucy<span>$</span>Zone)<br><br>  A   B   C   D   E<br><span>130</span> <span>295</span> <span>426</span>  <span>86</span>  <span>63</span></pre>
<p>The first step of MRP is Multilevel Regression in order to estimate post-strata means. The following code shows how to estimate them by using the <strong>lmer</strong> function. The object <strong>Mupred</strong> contains the corresponding $\mu_j$ ($j$ is defined for each level) regarding each stratum (zones).</p>
<pre><span>></span> # <span>Step</span> <span>1</span><span>:</span> <span>&lt;&lt;</span><span>MR</span><span>>></span> <span>-</span> <span>Multilevel</span> regression<br><span>></span> <span>M1</span> <span>&lt;-</span> lmer(<span>Income</span> <span>~</span> <span>Level</span> <span>+</span> (<span>1</span> <span>|</span> <span>Zone</span>), <span>data</span> <span>=</span> <span>SLucy</span>)<br><span>></span> coef(<span>M1</span>)<br><span>$</span><span>Zone</span><br>  (<span>Intercept</span>) <span>LevelMedium</span> <span>LevelSmall</span><br><span>A</span>    <span>1265.596</span>   <span>-</span><span>579.1851</span>  <span>-</span><span>893.8958</span><br><span>B</span>    <span>1138.337</span>   <span>-</span><span>579.1851</span>  <span>-</span><span>893.8958</span><br><span>C</span>    <span>1189.285</span>   <span>-</span><span>579.1851</span>  <span>-</span><span>893.8958</span><br><span>D</span>    <span>1248.658</span>   <span>-</span><span>579.1851</span>  <span>-</span><span>893.8958</span><br><span>E</span>    <span>1284.322</span>   <span>-</span><span>579.1851</span>  <span>-</span><span>893.8958</span><br><br>attr(,<span>&quot;class&quot;</span>)<br>[<span>1</span>] <span>&quot;coef.mer&quot;</span><br><span>></span> <span>SLucy</span><span>$</span><span>Pred</span> <span>&lt;-</span> predict(<span>M1</span>)<br><span>></span><br><span>></span> # <span>Summary</span><br><span>></span> grouped <span>&lt;-</span> group_by(<span>SLucy</span>, <span>Zone</span>, <span>Level</span>)<br><span>></span> <span>sum</span> <span>&lt;-</span> summarise(grouped, mean2 <span>=</span> mean(<span>Pred</span>))<br><span>></span> (<span>Mupred</span> <span>&lt;-</span> matrix(<span>sum</span><span>$</span>mean2, ncol <span>=</span> <span>5</span>, nrow <span>=</span> <span>3</span>))<br>          [,<span>1</span>]      [,<span>2</span>]      [,<span>3</span>]      [,<span>4</span>]      [,<span>5</span>]<br>[<span>1</span>,] <span>1265.5959</span> <span>1138.3370</span> <span>1189.2852</span> <span>1248.6575</span> <span>1284.3224</span><br>[<span>2</span>,]  <span>686.4107</span>  <span>559.1518</span>  <span>610.1000</span>  <span>669.4724</span>  <span>705.1373</span><br>[<span>3</span>,]  <span>371.7001</span>  <span>244.4412</span>  <span>295.3894</span>  <span>354.7618</span>  <span>390.4267</span><br></pre>
<p>Now we have estimated the post-strata means, it’s time to weight every strata by their corresponding counts in order to obtain an estimate of the mean income by zone. As we know each post-strata size, we simply use the function <strong>aggregate</strong> to obtain the MRP estimator of the parameters of interest.</p>
<pre><span>></span> <span># Step 2: &lt;&lt;P>> - Post-stratification</span><br><span>></span> <span># Mean income estimation per zone</span><br><span>></span> colSums(Np <span>*</span> Mupred) / table(Lucy$Zone)<br><br>       A        B        C        D        E<br>643.5724 312.8052 332.1726 682.8031 778.2428</pre>
<p>How accurate was the estimation? It was good compared to the true parameters on the finite population. </p>
<pre><span>></span> <span># True parameters</span><br><span>></span> aggregate(Lucy<span>$</span>Income, <span>by</span> <span>=</span> <span>list</span>(Lucy<span>$</span>Zone), <span>FUN</span> <span>=</span> mean)<br>  Group.1        x<br><span>1</span>       A <span>652.2834</span><br><span>2</span>       B <span>320.7469</span><br><span>3</span>       C <span>331.0195</span><br><span>4</span>       D <span>684.9821</span><br><span>5</span>       E <span>767.3879</span><br></pre>
<p> </p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://hagutierrezro.blogspot.com/2017/01/gelman-mrp-in-r-what-is-this-all-about.html"> Data Literacy - The blog of Andrés Gutiérrez</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/J1gbajtwKW4"></a></div>