2
$\begingroup$

I have some power-law data sets coming from an Ecology study. Some of them seem to be best modeled using linear regression after log-transformation of the data, but other data sets seem to be best modeled using nonlinear regression. A third group of data sets could be modeled using an average model.

For every data set, I calculate the likelihood that the data was generated from a normal distribution with additive error ($y=Bx^\beta+\epsilon$ with $\epsilon\sim\mathcal{N}(0,\sigma^2)$) and the likelihood that the data was generated from a lognormal distribution with multiplicative error ($y=Bx^\beta e^\epsilon$ with $\epsilon\sim\mathcal{N}(0,\sigma^2)$, so, in the logarithmic scale, $\log(y)=\log(B)+\beta\log(x)+\epsilon$). Next, I calculate AICc for each model, compare them and if neither model is clearly favored, I get the AICc weights of the two models, calculate the model-weighted power-law parameters and, finally, I get their CIs by bootstrapping.

I would like to calculate the coefficient of determination for such average model. Not completely sure about, but I think the best approach in this situation is to apply the generalized $R^2$ that, for continuous models, was defined as [Magee, 1990]:

$$ R^2 = 1 - \left(\frac{\mathcal{L}(0)}{\mathcal{L}(\hat\theta)}\right)^\frac{2}{n} $$

where $\mathcal{L}(\hat\theta)$ and $\mathcal{L}(0)$ denote the likelihoods of the fitted and the 'null' model, respectively, and $n$ is the sample size. In terms of the loglikelihoods, the generalized coefficient of determination would be:

$$ R^2 = 1 - \mathrm{e}^{-\frac{2}{n}\left(\mathcal{l}(\hat\theta)-\mathcal{l}(0)\right)} $$

In this situation, the likelihood related with the error structure with normal distribution becomes:

$$\mathcal{L}_\text{norm}(\hat\theta) = \prod_{i=1}^n\left[\frac{1}{\sqrt{2\pi\sigma^2_\text{norm}}}\;\exp{\left(-\frac{\left(y_i-B_\text{norm}x_i^{\beta_\text{norm}}\right)^2}{2\sigma^2_\text{norm}}\right)}\right]$$

So, the loglikelihood is: \begin{eqnarray*} \mathcal{l}_\text{norm}(\hat\theta) &=& -\frac{n}{2}\log\left|2\pi\sigma^2_\text{norm}\right| - \frac{1}{2\sigma^2_\text{norm}}\underbrace{\sum_{i=1}^n\left(y_i-B_\text{norm}x_i^{\beta_\text{norm}}\right)^2}_{\mathrm{RSS }_\text{norm}}\\ &=& -\frac{n}{2}\log\left|2\pi\sigma^2_\text{norm}\right|-\frac{\mathrm{RSS}_\text{norm}}{2\sigma^2_\text{norm}} \end{eqnarray*}

On the other hand, the likelihood related with the error structure with lognormal distribution would be:

$$\mathcal{L}_\text{logn}(\hat\theta) = \prod_{i=1}^n\left[\frac{1}{y_i\sqrt{2\pi\sigma^2_\text{LR}}}\;\exp{\left(-\frac{\left(\log|y_i|-\log|B_\text{logn}|-\beta_\text{logn}\log|x_i|\right)^2}{2\sigma^2_\text{logn}}\right)}\right]$$

So, the loglikelihood is: \begin{eqnarray*} \mathcal{l}_\text{logn}(\hat\theta) &=& -\frac{n}{2}\log\left|2\pi\sigma^2_\text{logn}\right| - \sum_{i=1}^n\log|y_i| -\\ &&\qquad-\frac{1}{2\sigma^2_\text{logn}}\underbrace{\sum_{i=1}^n\left(\log|y_i|-\log|B_\text{logn}|-\beta_\text{logn}\log|x_i|\right)^2}_{\mathrm{RSS }_\text{logn}}\\ &=& -\frac{n}{2}\log\left|2\pi\sigma^2_\text{logn}\right| - \frac{\mathrm{RSS}_\text{logn}}{2\sigma^2_\text{logn}} - \sum_{i=1}^n\log|y_i| \end{eqnarray*}

But what about the 'null' models? I understand that they are the models with only the intercept. So for the gaussian additive error model:

$$\mathcal{L}_\text{norm}(0) = \prod_{i=1}^n\left[\frac{1}{\sqrt{2\pi\sigma^2_\text{norm0}}}\;\exp{\left(-\frac{\left(y_i-\bar y\right)^2}{2\sigma^2_\text{norm0}}\right)}\right]$$

So: \begin{eqnarray*} \mathcal{l}_\text{norm}(0) &=& -\frac{n}{2}\log\left|2\pi\sigma^2_\text{norm0}\right| - \frac{1}{2\sigma^2_\text{norm0}}\sum_{i=1}^n\left(y_i-\bar y\right)^2\\ &=& -\frac{n}{2}\left(\log\left|2\pi\sigma^2_\text{norm0}\right| + 1 \right) \end{eqnarray*}

since $\sigma^2_\text{norm0}=\frac{1}{n}\sum\left(y_i-\bar y\right)^2=\frac{1}{n}\mathrm{TSS}_\text{norm}$. Now, coming back to the coefficient of determination, we have:

\begin{eqnarray*} R^2_\text{norm} &=& 1 - \mathrm{e}^{\frac{2}{n}\left(\mathcal{l}_\text{norm}(0)-\mathcal{l}_\text{norm}(\hat\theta)\right)} = 1 - \exp{\left(\frac{\log(\mathrm{RSS}_\text{norm})}{\log(\mathrm{TSS}_\text{norm})}\right)} =\\ &=& 1 - \frac{\mathrm{RSS}_\text{norm}}{\mathrm{TSS}_\text{norm}} = 1 - \frac{\sum_{i=1}^n\left(y_i-B_\text{norm}x_i^{\beta_\text{norm}}\right)^2}{\sum_{i=1}^n\left(y_i-\bar y\right)^2} \end{eqnarray*}

recovering the traditional expression for $R^2$. Using the same approach for calculating $R^2_\text{logn}$, then:

$$\mathcal{L}_\text{logn}(0) = \prod_{i=1}^n\left[\frac{1}{y_i\sqrt{2\pi\sigma^2_\text{logn0}}}\;\exp{\left(-\frac{\left(\log|y_i|-\log|B_\text{logn0}|\right)^2}{2\sigma^2_\text{logn0}}\right)}\right]$$

So: \begin{eqnarray*} \mathcal{l}_\text{logn}(0) &=& -\frac{n}{2}\log\left|2\pi\sigma^2_\text{logn0}\right| - \frac{1}{2\sigma^2_\text{logn0}}\sum_{i=1}^n\left(\log|y_i|-\overline{\log|y|}\right)^2 - \sum_{i=1}^n\log|y_i|\\ &=& -\frac{n}{2}\left(\log\left|2\pi\sigma^2_\text{logn0}\right| + 1 \right) - \sum_{i=1}^n\log|y_i| \end{eqnarray*}

since $\sigma^2_\text{logn0}=\frac{1}{n}\sum\left(\log|y_i|-\overline{\log|y|}\right)^2=\frac{1}{n}\mathrm{TSS}_\text{logn}$. Again, recalling the expression for the generalized coefficient of determination, we have:

\begin{eqnarray*} R^2_\text{logn} &=& 1 - \mathrm{e}^{\frac{2}{n}\left(\mathcal{l}_\text{logn}(0)-\mathcal{l}_\text{logn}(\hat\theta)\right)} = 1 - \exp{\left(\frac{\log(\mathrm{RSS}_\text{logn})}{\log(\mathrm{TSS}_\text{logn})}\right)} =\\ &=& 1 - \frac{\mathrm{RSS}_\text{logn}}{\mathrm{TSS}_\text{logn}} = 1 - \frac{\sum_{i=1}^n\left(\log|y_i|-\log|B_\text{logn}|-\beta_\text{logn}\log|x_i|\right)^2}{\sum_{i=1}^n\left(\log|y_i|-\overline{\log|y|}\right)^2} \end{eqnarray*}

So the question is: how to calculate $R^2$ but for the average model? Many thanks in advance for any comment, suggestion, reference and, of course, answer!

$\endgroup$
  • 1
    $\begingroup$ Average of what model? Do you mean 1) an autocorrelation model? 2) a population model, or model of average responses or, 3) a single response modeled by the average of several models, as for example, in mixture models. $\endgroup$ – Carl Aug 15 '16 at 22:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.