doPRC perform a Principal Response Curve analysis (PRC) using the
S3 formula interface of rda or cca
doPRC(
formula,
scale = FALSE,
referencelevel = NULL,
rank = 2,
flip = rep(FALSE, rank),
scaling = "ms",
method = "rda",
data
)Arguments
- formula
formula specifying how the response (species data) must modeled by predictors and covariates. A typical
formulaisY~treatment*week+Condition(week), whereYis the global response data matrix or data frame. Typically these are community (species) data. If count-like, they probably should be log transformed prior to the analysis, because a multiplicative model has more realism for non-negative data than an additive model.- scale
Scales species to unit variance (like correlations) if
methodis"rda"(default:FALSE), else void. Seerda.- referencelevel
numeric or character level(s) to be used as reference level or levels of the focal treatment(s); default: 1 (first level(s)). With
doPRC, the focal treatment is (or treatments are) given in the resulting list atresult$focal_and_conditioning_factors$`focal factor`- rank
number of axes to be computed; default 2.
- flip
logical value or vector; default: FALSE. Should the axes be flipped, i.e. reversed in orientation?
flipcan be numeric with values -1 and 1 for reverse and no reverse, respectively, i.e. multiply the scores by either -1 or 1.- scaling
character
c("ms", "ss"). Scales the axes of an RDA object to mean square species scores of 1 ("ms") or a sum of squares of 1 ("ss"). The only scaling for a CCA object is"ss"(weighted mean square = 1, with weight = species total abundance). For a dbrda object, there are no species scores and the scaling is like "ss". The default scaling ("ms"for RDA and"ss"for CCA) has the advantage that the PRC treatment and sample scores measure the effect sizes of the 'average' species as their (weigthed) mean square is 1 (sd = 1); it is the scaling used in Canoco5 (www.canoco5.com). Scaling"ss"is the default inprcompand used in ter Braak (2023).- method
character rda or cca. Default:
"rda"- data
Data frame containing the variables on the right hand side of the model formula.
Value
an object as a cca.object with additional entries as output of PRC_scores.
The assigned classes are c("PRC","rda","cca".
Details
The function performs PRC using rda or cca with
the specified formula and applies PRC_scores to the result.
References
ter Braak (2023) Redundancy analysis includes analysis of variance-simultaneous component analysis (ASCA) and outperforms its extensions Chemometrics and Intelligent Laboratory Systems https://doi.org/10.1016/j.chemolab.2023.104898
ter Braak, C.J.F. & te Beest, D.E. 2022. Testing environmental effects on taxonomic composition with canonical correspondence analysis: alternative permutation tests are not equal. Environmental and Ecological Statistics. 29 (4), 849-868. https://doi.org/10.1007/s10651-022-00545-4
van den Brink, P.J. & ter Braak, C. (1998) Multivariate analysis of stress in experimental ecosystems by Principal Response Curves and similarity analysis. Aquatic Ecology, 32, 163-178. http://dx.doi.org/10.1023/A:1009944004756
van den Brink, P.J. & ter Braak, C.J.F. (1999) Principal Response Curves: Analysis of time-dependent multivariate responses of a biological community to stress. Environmental Toxicology and Chemistry, 18, 138-148. https://doi.org/10.1002/etc.5620180207
van den Brink, P.J., den Besten, P., bij de Vaate, A. & ter Braak, C.J.F. (2009) Principal response curves technique for the analysis of multivariate biomonitoring time series. Environmental Monitoring and Assessment, 152, 271-281. http://dx.doi.org/10.1007/s10661-008-0314-6
Examples
library(PRC)
data(pyrifos, package = "vegan") #log-transformed species data from package vegan
#the chlorpyrifos experiment from van den Brink & ter Braak 1999
Design <- data.frame(week=gl(11, 12, labels=c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24)),
dose=factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11)),
ditch = gl(12, 1, length=132))
#
mod_prc <- doPRC(pyrifos ~ dose:week + Condition(week), data = Design)
#>
#> Some constraints or conditions were aliased because they were redundant. This
#> can happen if terms are constant or linearly dependent (collinear):
#> 'dose44:week-4', 'dose44:week-1', 'dose44:week0.1', 'dose44:week1',
#> 'dose44:week2', 'dose44:week4', 'dose44:week8', 'dose44:week12',
#> 'dose44:week15', 'dose44:week19', 'dose44:week24'
mod_prc$focal_and_conditioning_factors$`focal factor` # as expected: "dose"
#> [1] "dose"
print(mod_prc)
#>
#> Call: rda(formula = pyrifos ~ dose:week + Condition(week), data = data,
#> scale = scale)
#>
#> Inertia Proportion Rank
#> Total 288.9920 1.0000
#> Conditional 63.3493 0.2192 10
#> Constrained 96.6837 0.3346 44
#> Unconstrained 128.9589 0.4462 77
#>
#> Inertia is variance
#>
#> -- NOTE:
#> Some constraints or conditions were aliased because they were redundant.
#> This can happen if terms are constant or linearly dependent (collinear):
#> 'dose44:week-4', 'dose44:week-1', 'dose44:week0.1', 'dose44:week1',
#> 'dose44:week2', 'dose44:week4', 'dose44:week8', 'dose44:week12',
#> 'dose44:week15', 'dose44:week19', 'dose44:week24'
#>
#> Eigenvalues for constrained axes:
#> RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11
#> 25.282 8.297 6.044 4.766 4.148 3.857 3.587 3.334 3.087 2.551 2.466
#> RDA12 RDA13 RDA14 RDA15 RDA16 RDA17 RDA18 RDA19 RDA20 RDA21 RDA22
#> 2.209 2.129 1.941 1.799 1.622 1.579 1.440 1.398 1.284 1.211 1.133
#> RDA23 RDA24 RDA25 RDA26 RDA27 RDA28 RDA29 RDA30 RDA31 RDA32 RDA33
#> 1.001 0.923 0.862 0.788 0.750 0.712 0.685 0.611 0.584 0.537 0.516
#> RDA34 RDA35 RDA36 RDA37 RDA38 RDA39 RDA40 RDA41 RDA42 RDA43 RDA44
#> 0.442 0.417 0.404 0.368 0.340 0.339 0.306 0.279 0.271 0.205 0.179
#>
#> Eigenvalues for unconstrained axes:
#> PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
#> 17.156 9.189 7.585 6.064 5.730 4.843 4.518 4.105
#> (Showing 8 of 77 unconstrained eigenvalues)
#>
# Proporti Rank
# Total 1.0000
# Conditional 0.2192 10 (between-week variation)
# Constrained 0.3346 44 (between-treatment variation across weeks)
# Unconstrained 0.4462 77 (within-week variation)
# without lines, threshold is 7
plotPRC(mod_prc)
#> The analysis is based on model
#> pyrifos~dose:week + Condition(week)
#> The current score, condition and treatment(s) are:
#> score (y-axis): PRC1
#> condition (x-axis): week
#> treatments (color ): dose
#> Set these arguments yourself if you wish otherwise.
#> 34 species with names out of 178 species
# this threshold is about P = 0.01 for PRC1 in the regression "species_k ~ week+PRC1"
# and retains 34 species to be plotted with names
sum(mod_prc$species[,"Fratio1"] >7) # 34
#> [1] 34
# only the PRC curves:
plotPRC(mod_prc, plot = 0, symbols_on_curves = TRUE)
#> The analysis is based on model
#> pyrifos~dose:week + Condition(week)
#> The current score, condition and treatment(s) are:
#> score (y-axis): PRC1
#> condition (x-axis): week
#> treatments (color ): dose
#> Set these arguments yourself if you wish otherwise.
#> 34 species with names out of 178 species
# with lines and a stronger threshold
plotPRC(mod_prc, plot = "ditch", threshold = 10,width = c(4,1))
#> The analysis is based on model
#> pyrifos~dose:week + Condition(week)
#> The current score, condition and treatment(s) are:
#> score (y-axis): PRC1
#> condition (x-axis): week
#> treatments (color ): dose
#> Set these arguments yourself if you wish otherwise.
#> 29 species with names out of 178 species
# modifying the plot
gg <- plotPRC(mod_prc, plot = "ditch",width = c(4,1), verbose = FALSE)
# PRC plot of samples (c_dt)
p1 <- gg$separateplots$treatments +
ggplot2::ggtitle(paste("new title:", latex2exp::TeX("$c_{dt}$")))
p2 <- gg$separateplots$species +
ggplot2::ylab("new title: loadings")# loadings of species (b_k)
# Assign these plots to symbols and use grid.arrange to produce the plot
#you like, for example:
gridExtra::grid.arrange(p1+ ggplot2::ylab("")+
ggplot2::xlab("week since chorpyrifos application") + ggplot2::ggtitle(""),
p2, ncol =2, widths = c(4,1),
top = "new title", left ="PRC curves of treaments",
right = "PRC species loadings")
# test significance of PRC/RDA axes ---------------------------------------
## Ditches are randomized, we have a time series, and are only
## interested in the first axis
library(vegan)
#> Loading required package: permute
ctrl <- how(plots = Plots(strata = Design$ditch,type = "free"),
within = Within(type = "none"), nperm = 99)
anova(mod_prc, by = "axis",permutations= ctrl, model = "reduced",
cutoff = 0.10)
#> Permutation test for rda under reduced model
#> Forward tests for axes
#> Plots: Design$ditch, plot permutation: free
#> Permutation: none
#> Number of permutations: 99
#>
#> Model: rda(formula = pyrifos ~ dose:week + Condition(week), data = data, scale = scale)
#> Df Variance F Pr(>F)
#> RDA1 1 25.282 15.0958 0.01 **
#> RDA2 1 8.297 5.0183 0.56
#> RDA3 1 6.044 3.6089
#> RDA4 1 4.766 2.8459
#> RDA5 1 4.148 2.4768
#> RDA6 1 3.857 2.3029
#> RDA7 1 3.587 2.1420
#> RDA8 1 3.334 1.9907
#> RDA9 1 3.087 1.8434
#> RDA10 1 2.551 1.5233
#> RDA11 1 2.466 1.4727
#> RDA12 1 2.209 1.3189
#> RDA13 1 2.129 1.2710
#> RDA14 1 1.941 1.1589
#> RDA15 1 1.799 1.0740
#> RDA16 1 1.622 0.9683
#> RDA17 1 1.579 0.9426
#> RDA18 1 1.440 0.8600
#> RDA19 1 1.398 0.8350
#> RDA20 1 1.284 0.7666
#> RDA21 1 1.211 0.7230
#> RDA22 1 1.133 0.6768
#> RDA23 1 1.001 0.5977
#> RDA24 1 0.923 0.5511
#> RDA25 1 0.862 0.5150
#> RDA26 1 0.788 0.4705
#> RDA27 1 0.750 0.4476
#> RDA28 1 0.712 0.4253
#> RDA29 1 0.685 0.4088
#> RDA30 1 0.611 0.3648
#> RDA31 1 0.584 0.3490
#> RDA32 1 0.537 0.3206
#> RDA33 1 0.516 0.3080
#> RDA34 1 0.442 0.2637
#> RDA35 1 0.417 0.2492
#> RDA36 1 0.404 0.2410
#> RDA37 1 0.368 0.2196
#> RDA38 1 0.340 0.2028
#> RDA39 1 0.339 0.2025
#> RDA40 1 0.306 0.1824
#> RDA41 1 0.279 0.1665
#> RDA42 1 0.271 0.1621
#> RDA43 1 0.205 0.1222
#> RDA44 1 0.179 0.1069
#> Residual 77 128.959
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1