Skip to content

Commit

Permalink
Merge branch 'restore-partial-R2-adj'
Browse files Browse the repository at this point in the history
This merge re-enables calculation of adj-R2 for partial RDA that was
removed in vegan release 2.5-1. It also implements adj-R2 for partial
CCA in a way that is consistent with RDA and varpart. Before 2.5-1
release CCA was inconsistent.

See github issue #295
  • Loading branch information
Jari Oksanen committed Jan 29, 2019
2 parents 6617fca + 5bfa96f commit 40584a1
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 19 deletions.
38 changes: 23 additions & 15 deletions R/RsquareAdj.R
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,7 @@
r2
}

## rda. Function varpart() estimates adjusted R2 for partial model
## rda(Y, X1, X2) as a difference of adjusted R2 values of models
## rda(Y, cbind(X1,X2)) and rda(Y, X2), and we did the same here from
## vegan release 2.0-3 to 2.4. However, I am not quite convinced that
## this is correct, and therefore I disable adjusted R2 for partial
## models while I study the issue (Jari Oksanen, 23/2/2018). If you
## change handling of partial models, remember to update ordiR2step!
## Use this with rda() results
`RsquareAdj.rda` <-
function(x, ...)
{
Expand All @@ -30,24 +24,38 @@
if (is.null(x$pCCA)) {
radj <- RsquareAdj(R2, n, m)
} else {
radj <- NA
## Partial model: same adjusted R2 as for component [a] in two
## source varpart model
R2p <- x$pCCA$tot.chi/x$tot.chi
p <- x$pCCA$QR$rank
radj <- RsquareAdj(R2 + R2p, n, m + p) - RsquareAdj(R2p, n, p)
}
list(r.squared = R2, adj.r.squared = radj)
}

## cca: for pCCA, see comment for rda. Earlier we used the same
## equation both for pCCA and CCA which was inconsistent with rda. git
## commit 9c25e6fdb4 has a version that was consistent with rda in
## vegan 2.0-3 to 2.4-6.
## cca result: RsquareAdj is calculated similarly as in
## varpart(). This is similar "semipartial" model as for rda() and
## found as a difference of R2-adj values of combined model with
## constraints + conditions and only conditions.
`RsquareAdj.cca` <-
function (x, permutations = 1000, ...)
{
r2 <- x$CCA$tot.chi / x$tot.chi
if (is.null(x$pCCA) || x$pCCA$rank == 0) {
if (is.null(x$pCCA)) {
p <- permutest(x, permutations, ...)
radj <- 1 - ((1 - r2) / (1 - mean(p$num) / x$tot.chi))
radj <- 1 - ((1 - r2) / (1 - mean(p$num / x$tot.chi)))
} else {
radj <- NA
p <- getPermuteMatrix(permutations, nobs(x))
Y <- ordiYbar(x, "initial")
r2tot <- (x$pCCA$tot.chi + x$CCA$tot.chi) / x$tot.chi
r2null <- mean(sapply(seq_len(nrow(p)), function(i)
sum(qr.fitted(x$CCA$QR, Y[p[i,],])^2)))
r2tot <- 1 - ((1-r2tot)/(1-r2null/x$tot.chi))
r2p <- x$pCCA$tot.chi / x$tot.chi
r2null <- mean(sapply(seq_len(nrow(p)), function(i)
sum(qr.fitted(x$pCCA$QR, Y[p[i,],])^2)))
r2p <- 1 - ((1-r2p)/(1-r2null/x$tot.chi))
radj <- r2tot - r2p
}
list(r.squared = r2, adj.r.squared = radj)
}
Expand Down
8 changes: 4 additions & 4 deletions man/RsquareAdj.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -32,9 +32,8 @@ Adjusted R-square
R-squared for a cca. The permutations can be calculated in parallel by
specifying the number of cores which is passed to \code{\link{permutest}}}

\item{\dots}{ Other arguments (ignored) except in the \code{cca}
method, where these are passed to \code{\link{permutest.cca}}.
}
\item{\dots}{ Other arguments (ignored) except in the case of cca in
which these arguments are passed to \code{\link{permutest}}.}
}

\details{ The default method finds the adjusted \eqn{R^2}{R-squared}
Expand All @@ -50,7 +49,8 @@ Adjusted R-square
\code{\link{glm}}.

The adjusted, \eqn{R^2}{R-squared} of \code{cca} is computed using a
permutation approach developed by Peres-Neto et al. (2006).
permutation approach developed by Peres-Neto et al. (2006). By
default 1000 permutations are used.

}

Expand Down

3 comments on commit 40584a1

@CamillaLienart
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hello !
I am trying to get the Adjusted R2 for a RDA analysis but when using the following script it always gives me 1 as R2 result and NA as adjusted r2 ... for any RDA based on any dataset... The vegan package I am using is the 2.5-6. I wonder if there is another way to calculate the adjusted R2 (or at least some r2) ... the varpart doesn't work in my case because we only have 1 explanatory matrix so it gives an error message that it needs more data tables...

Here is the script and results for the adjusted R2

RsquareAdj(rda_Pontaillac_Full)
$r.squared
[1] 1

$adj.r.squared
[1] NA

Thank you for your help.
Camilla

@jarioksa
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should happen (it is so coded) if you have a saturated model which has no residual variation after constraints. Saturated model means that for data of n observations the number of explanatory terms is n–1 or greater. I don't have your data, but if the unadjusted R2 = 1, you have no residual variation and adjusted R2 cannot be meaningfully calculated.

@CamillaLienart
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, Indeed, I see... so it comes from my data which has 12 rows for 26 explanatory variables... I have tried to use vif.cca() to select some of the explanatory variables based on Brocard, Gillet, Legendre 2011. But still gives R2=1 and NA for the adjusted. I am not so sure what to do next...

Please sign in to comment.