title: Random Effects Esimation on cluter sample using plm on a unbalanced panel
body: <p>I normally work within a generalized least squares framework estimating, what Wooldridge's <em>Introductory</em> (2013) calls, Random Effects and Fixed Effects models on longitudinal data indexed by an individual and a time dimension.</p>
  
  <p>I've been using the <em>Feasible GLS</em> estimation in <code>plm()</code>, from the <a href="/questions/tagged/plm" class="post-tag" title="show questions tagged &#39;plm&#39;" rel="tag">plm</a> package, to estimate the Random Effects Model – what some stats literature term the <em>Mixed Model</em>.  The <code>plm()</code> function takes an <code>index</code> argument where I indicate the individual and time indexes. However, I’m now faced with some data where each individual has several measures at each time-point, i.e. what a group-wise structure.</p>
  
  <p>I’ve found out that it’s possible to set up such a model using <code>lmer()</code> from the <a href="/questions/tagged/lme4" class="post-tag" title="show questions tagged &#39;lme4&#39;" rel="tag">lme4</a> package, however I am a bit confused by the differences in jargon, and also the likelihood framework, and I wanted to know if specified the model correctly. I fear I could overlook at more substantial as I am not familiar with the framework and this terminology.</p>
  
  <p>I can replicate my usual  <code>plm()</code> model using <code>lmer()</code>, but I am unsure as to how I could add the grouping. I’ve tried to illustrate my question in the following.</p>
  
  <p>I found some data that looks somewhat like my data to illustrate my situation. First some packages that are needed,</p>
  
  <pre><code>install.packages(c("mlmRev", "plm", "lme4", "stargazer"), dependencies = TRUE)
  </code></pre>
  
  <p>and then the data</p>
  
  <pre><code>data(egsingle, package = "mlmRev")
  </code></pre>
  
  <p><code>egsingle</code> is a unbalanced panel consisting of 1721 school children, grouped in 60 schools, across five time points. These data are originally distributed with the HLM software package (Bryk, Raudenbush and Congdon, 1996), but can be found the <a href="/questions/tagged/mlmrev" class="post-tag" title="show questions tagged &#39;mlmrev&#39;" rel="tag">mlmrev</a> package, for details see <code>? mlmRev::egsingle</code></p>
  
  <p>Some light data management</p>
  
  <pre><code>dta &lt;- egsingle
  dta$Female &lt;- with(dta, ifelse(female == 'Female', 1, 0))
  </code></pre>
  
  <p>Here’s a snippet for the data</p>
  
  <pre><code>dta[118:127,c('schoolid','childid','math','year','size','Female')]
  #&gt;     schoolid   childid   math year size Female
  #&gt; 118     2040 289970511 -1.830 -1.5  502      1
  #&gt; 119     2040 289970511 -1.185 -0.5  502      1
  #&gt; 120     2040 289970511  0.852  0.5  502      1
  #&gt; 121     2040 289970511  0.573  1.5  502      1
  #&gt; 122     2040 289970511  1.736  2.5  502      1
  #&gt; 123     2040 292772811 -3.144 -1.5  502      0
  #&gt; 124     2040 292772811 -2.097 -0.5  502      0
  #&gt; 125     2040 292772811 -0.316  0.5  502      0
  #&gt; 126     2040 293550291 -2.097 -1.5  502      0
  #&gt; 127     2040 293550291 -1.314 -0.5  502      0
  </code></pre>
  
  <p>Here’s how I would set a <em>random effects model</em> without the  <code>schoolid</code> using <code>plm()</code>,</p>
  
  <pre><code>library(plm)
  reg.re.plm &lt;- plm(math~Female+size+year, dta, index = c("childid", "year"), model="random")
  # summary(reg.re.plm)
  </code></pre>
  
  <p>I can reproduce these results <code>lme4</code> like this</p>
  
  <pre><code>require(lme4)
  dta$year &lt;- as.factor(dta$year)
  reg.re.lmer &lt;- lmer(math~Female+size+year+(1|childid), dta)
  # summary(reg.re.lmer)
  </code></pre>
  
  <p>Now, from reading <a href="http://lme4.r-forge.r-project.org/book/" rel="nofollow noreferrer">chapter 2 in Bates (2010) “lme4: Mixed-effects modeling
  with R”</a> I believe I’ve this is how I would specific the model including the cluster level,  <code>schoolid</code>,</p>
  
  <pre><code>reg.re.lmer.in.school &lt;- lmer(math~Female+size+year+(1|childid)+(1|schoolid), dta)
  # summary(reg.re.lmer.in.school)
  </code></pre>
  
  <p>However, when I look at the results I am not too convinced I’ve actually specified it correctly (see below).</p>
  
  <p>In my actual data the repeated measures are within individuals, but I take that I can use this data as example. I would appreciate any advice on how to proceed. Maybe a reference to a worked example with notation/terminology not too far from what is used in Wooldridge (2013). And, how do I work backwards and write up the specification for the <code>reg.re.lmer.in.school</code> model?</p>
  
  <pre><code># library(stargazer)
  stargazer::stargazer(reg.re.plm, reg.re.lmer, reg.re.lmer.in.school, type="text")
  #&gt; =====================================================================
  #&gt;                                    Dependent variable:               
  #&gt;                     -------------------------------------------------
  #&gt;                                           math                       
  #&gt;                                panel                   linear        
  #&gt;                               linear                mixed-effects    
  #&gt;                                 (1)                (2)        (3)    
  #&gt; ---------------------------------------------------------------------
  #&gt; Female                        -0.025              -0.025     0.008   
  #&gt;                               (0.046)            (0.047)    (0.042)  
  #&gt;                                                                      
  #&gt; size                        -0.0004***          -0.0004***  -0.0003  
  #&gt;                              (0.0001)            (0.0001)   (0.0002) 
  #&gt;                                                                      
  #&gt; year-1.5                     0.878***            0.876***   0.866*** 
  #&gt;                               (0.059)            (0.059)    (0.059)  
  #&gt;                                                                      
  #&gt; year-0.5                     1.882***            1.880***   1.870*** 
  #&gt;                               (0.059)            (0.058)    (0.058)  
  #&gt;                                                                      
  #&gt; year0.5                      2.575***            2.574***   2.562*** 
  #&gt;                               (0.059)            (0.059)    (0.059)  
  #&gt;                                                                      
  #&gt; year1.5                      3.149***            3.147***   3.133*** 
  #&gt;                               (0.060)            (0.059)    (0.059)  
  #&gt;                                                                      
  #&gt; year2.5                      3.956***            3.954***   3.939*** 
  #&gt;                               (0.060)            (0.060)    (0.060)  
  #&gt;                                                                      
  #&gt; Constant                     -2.671***          -2.669***  -2.693*** 
  #&gt;                               (0.085)            (0.086)    (0.152)  
  #&gt;                                                                      
  #&gt; ---------------------------------------------------------------------
  #&gt; Observations                   7,230              7,230      7,230   
  #&gt; R2                             0.735                                 
  #&gt; Adjusted R2                    0.735                                 
  #&gt; Log Likelihood                                  -8,417.815 -8,284.357
  #&gt; Akaike Inf. Crit.                               16,855.630 16,590.720
  #&gt; Bayesian Inf. Crit.                             16,924.490 16,666.460
  #&gt; F Statistic         2,865.391*** (df = 7; 7222)                      
  #&gt; =====================================================================
  #&gt; Note:                                     *p&lt;0.1; **p&lt;0.05; ***p&lt;0.01
  </code></pre>