Hi, thank you both for your replies.  However I remain a little confused.  Firstly apologies as I noticed that in my original posting I had left out the sqrt in the formula for the residual standard error given by Sprugel (1983) - RSE = sqrt((sum(residuals^2))/(N - no.of parameters)).  In addition to a random effects term for the intercept my model (see below) also has a random effects term for both explanatory variables (dbh and height) so I guess this changes things and makes it more complicated? 

Although I understand what the code you suggested is doing (summing the variances) I'm afraid my understanding of statistics is perhaps not good enough to know how and if this is equivalent to the formula given by Sprugel (which seems to be the accepted way of working out the residual error and correcting estimates for log-bias).  Can Sprugels formula as it's given not be used with mixed effects models?  I used resid(model) to grab the residuals for use in the above formula, maybe this wasn't correct?  I'm unclear as to why all the variances should be used to correct the final estimates of biomass for each tree and not just the residual variance (as in Sprugel)?

When I work out the CF using the above formula I get a value of ~1.0444 (very reasonable) and using the code you sent me I get a value of 1.529! (not very reasonable and quite worrying).  So any further help in explaining this would be much appreciated as I am now quite concerned about the goodness-of-fit of my model if the value of ~1.5 is correct.

Many thanks for your help,


Rachel


model<-lmer((log.mass)~centre.log.dbh+centre.log.height+ (1+centre.log.dbh+centre.log.height|species_site),data=allometry_2,REML=T) 


________________________________
 From: "Baldwin, Jim -FS" <jbaldwin@fs.fed.us>

r-sig-mixed-models@r-project.org" <r-sig-mixed-models@r-project.org> 
Sent: Friday, 3 August 2012, 2:37
Subject: RE: [R-sig-ME] basic question

If the objective is to obtain an approximately unbiased estimate for the prediction of a future observation that would have the same variance components as found in the fitted model, then multiplying the backtransformed value by exp of "something like sum(sapply(VarCorr(model),diag))+sigma(model)^2 ..."/2 is a reasonable thing to do as Dr. Bolker suggests.

But if the approximately unbiased prediction is for a different situation which might not include either the same or all of the variance components in the fitted model, you'll need to think about which ones to use.

Also, (and this is likely very obvious) if the multiplicative correction is much larger than around 1.5, I'd have to wonder about the goodness-of-fit of the distributional assumptions of the model.  Sometimes this bias correction can give results that are way outside of observed values on the original scale.  (And, at least for me, I'd also check my code in such situations.  It's many times my fault and not the data's fault.)

Jim

Jim Baldwin
Station Statistician
Pacific Southwest Research Station
USDA Forest Service


-----Original Message-----
From: r-sig-mixed-models-bounces@r-project.org [mailto:r-sig-mixed-models-bounces@r-project.org] On Behalf Of Ben Bolker
Sent: Thursday, August 02, 2012 5:56 PM
To: Rachel Cohen; r-sig-mixed-models@r-project.org
Subject: Re: [R-sig-ME] basic question

[cc'ing back to r-sig-mixed-models]

On 12-08-02 07:14 AM, Rachel Cohen wrote:
> Hi, thanks for your reply.  I was reading the Pinheiro and Bates book
> (pg.21) which says that for a basic model the parameters are: β, σ2b ,
> and σ2 which is why I thought the residual variance was counted as a
> parameter?  Also, on the same page: "Although the random effects, bi,
> i = 1, . . .,M may behave like parameters, formally they are just
> another level of random variation in the model so we do not “estimate”
> them as such."  That lead me to believe that the random effects are
> not counted as parameters?

   First, you need to make the distinction between the b_i value ('BLUPs', or 'conditional modes', or 'random effect estimates' -- the estimates of what's happening in each grouping level) -- and the theta values (the variances of the random effects).
>
> I am intending on using the RSE in the formula: CF = exp (RSE^2
> /2) which is supposed to correct for log-bias in estimates (in this
> case estimated tree biomass).  Each unlogged estimate is multiplied by
> the correction factor (CF).  I haven't come across any mention in the
> literature of whether it is appropriate to use the RSE from a
> mixed-effects model in this way.  I hope so?

  If you want to fit on the log scale but estimate the *mean* effects on the unlogged scale, you should (I think) add the exp(TOTAL_variance/2); in this case TOTAL_variance will be the sum of *all of the variance
components* in the model, including the residual variance (I'm going to assume you have only intercept random effects in the model, otherwise things would get a bit more complicated ...) I *think* it would be something like sum(sapply(VarCorr(model),diag))+sigma(model)^2 ...  I think there has been a reasonably recent post by Jarrod Hadfield answering someone's question related to this topic ...

  If I'm right, then there's no need to go messing around with numbers of parameters ...
  Ben Bolker





------------------------
> *From:* Ben Bolker <bbolker@gmail.com>
> *To:* r-sig-mixed-models@r-project.org
> *Sent:* Thursday, 2 August 2012, 2:51
> *Subject:* Re: [R-sig-ME] basic question
>
> Rachel Cohen <stat.list@...> writes:
>
>>  Hi, I have a very basic query (I think).  If I fit a model of the
>> specification given below how many model parameters are there?  Am I
>> correct in thinking that the model parameters are the fixed effects
>> intercept and coefficients (in this case two of them) and also a
>> parameter for within and for between group variance?  So in total 5
>> fitted parameters? or are there also covariance parameters?.  I am
>> trying to calculate the standard error of the residuals as:
>
>   Unfortunately, this is an easy question to state but not (in my
> opinion) an easy question to answer unambiguously.
>
>
>> RSE = ((sum(residuals^2))/(N - no.of parameters))
>>
>> and therefore need to know how many parameters my model actually has.
>> I've gotten a little confused by the
>> literature definitions of what constitutes a parameter.
>
> It's confusing!
>
> [snip]
>
> A *reasonable* definition (although not the only one) would be to
> count the number of fixed-effect parameters ('beta' in much of the
> literature) and the number of random-effect parameters (generally
> referred to as 'theta' in the lme4 documentation, but varying a great
> deal among references):
>
> length(fixef(model))+length(getME(model,"theta"))
>
>   in your case that's 3 parameters for the fixed effects and 6 RE
> parameters (you have a 3x3 variance-covariance matrix of the random
> effects, the matrix is symmetric, so counting the diagonal plus one
> triangle gives 3*(3+1)/2 = 6).  Generally one doesn't count the
> residual variance since that is estimated from the (penalized)
> residual sum of squares.
>
>   It really depends what you want to use the RSE for.  It may very
> well not have the properties you're expecting (i.e. the properties
> that it has in a simple (non-mixed) linear model ...)
>
>   There's a bit more about parameter-counting issues at
> http://glmm.wikidot.com/faq, I think ...
>
>> model<-lmer((log.mass)~centre.log.dbh+centre.log.height+
>> (1+centre.log.dbh+centre.log.height|species_site),
>> data=allometry_2,REML=T)
>>
>>  model summary output table:
>>
>> Linear mixed model fit by REML
>> Formula: (log.mass) ~ centre.log.dbh + centre.log.height +
>> (1 + centre.log.dbh + centre.log.height | species_site)
>
> [snip]
>> Random effects:
>>  Groups       Name              Variance  Std.Dev. Corr
>>  species_site (Intercept)       0.061636  0.24827
>>               centre.log.dbh    0.279362  0.52855   0.575
>>               centre.log.height 0.285753  0.53456  -0.441 -0.755
>>  Residual                       0.086963  0.29489
>> Number of obs: 337, groups: species_site, 18
>>
>
> _______________________________________________
> R-sig-mixed-models@r-project.org
> <mailto:R-sig-mixed-models@r-project.org> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models 
>
>

_______________________________________________
R-sig-mixed-models@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models 




This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.
	[[alternative HTML version deleted]]

