title: Iteratively extracting fixed effects coefficients from a lmer model for use in bootstrap estimation
body: <p>The following bit of code was developed in order to generate a c.i. for the location of a peak in a modeled relationship. A simple polynomial approach is used where the RD term and its interactions determines that location (the 'zeros' of the first derivative are what is being solved here: peak.temp=-(m1.coef[10]*0.16+m1.coef[9]<em>0.19)/(2</em>m1.coef[8]), from line 19 below).</p>
  <p>The code was initially developed for use with nlme/lme and it works fine in that setting. The issue I'm having is related to transitioning for use with lme4/lmer models, which are better suited to dealing with the random effects.</p>
  <p>The first error being thrown is: 'Error in m1.new$coefficients : $ operator not defined for this S4 class', corresponding to line=17 below. Searching on that error I (believe I) was able to determine that the syntax needs to be updated to work with the S4 class, and that appears to involve, at least in part, substituting a @ for the $. However, simply making that change, e.g. m1.new@coefficients@fixed leads to another error, which is what has me stumped. That message is: Error: no slot of name &quot;coefficients&quot; for this object of class &quot;lmerModLmerTest&quot;</p>
  <p>While I was able to get the available slot names for the model (m1), I'm not clear on which, if any of these, are what I should be substituting for the previous 'coefficients' and possibly also 'fixed' terms in that line of code.</p>
  <p>m1 slots are: &quot;vcov_varpar&quot;, &quot;Jac_list&quot;,&quot;vcov_beta&quot;,&quot;sigma&quot;,&quot;resp&quot;,&quot;Gp&quot;,&quot;call&quot;,&quot;frame&quot;,&quot;flist&quot;,&quot;cnms&quot;&quot;lower&quot;,&quot;theta&quot;,&quot;beta&quot;,&quot;u&quot;,&quot;devcomp&quot;,&quot;pp&quot;,&quot;optinfo&quot;</p>
  <p>I suspect the source of this issue will be obvious to someone having experience with both packages nlme/lme4, and hope I have provided enough background to lead to a solution. Thanks much in advance!</p>
  <p>#----------------------------------------------</p>
  <pre><code>### bootstrap for m1 ###
  nreps=1000
  df.rows=dim(data)[1]
  peak.m1=numeric(nreps)
  for (i in 1:nreps){
    # resample the rows of the data frame
    new.rows=sample(1:df.rows,df.rows,replace=T)
    new.df=data[new.rows,]
    # fit the model
    m1.new=lmer(gPAI.on ~ WAM.sc + H + EH.sc + skew.sc + QMD.R.14.sc + RD + I(RD^2) +
                 EH.sc*RD + skew.sc*RD + skew.sc*EH.sc +
                 (1|mu/plot) + (1|site.yr),
               weights=plot.wt,
               control = lmerControl(optimizer=&quot;bobyqa&quot;),
               REML=F, data = new.df)
    # extract the fixed coefficients
  m1.coef=m1.new$coefficients$fixed
    # compute the peak
    peak.temp=-(m1.coef[10]*0.16+m1.coef[9]*0.19)/(2*m1.coef[8])
    #save the peak
    peak.m1[i]=peak.temp
  }
  plot(density(peak.m1))
  # find the mean and se for the m1 peak
  (mean.sd=c(mean(peak.m1),sd(peak.m1)))
  # a 95% confidence interval (assuming normality)
  c(mean.sd[1]-2*mean.sd[2],mean.sd[1]+2*mean.sd[2])
  # an empirical CI
  quantile(peak.m1,c(0.025,0.975))
  # check the normality of the peak estimates
  shapiro.test(peak.m1) 
  </code></pre>