| Title: | Principal Response Curves (PRC) For Displaying Multivariate Main Effects And Interactions |
|---|---|
| Description: | PRC is a small \code{vegan}- and \code{LMMsolver}-based package to create and display Principal Response Curves (PRC) and smooth Principal Response Curves (PRC). PRC is a special case of Redundancy Analysis (rda) for multivariate responses in designed experiments and repeated measurement designs. PRC expresses ordination scores as deviation from a reference, usually to visualize multivariate A-dependent effects (e.g time-dependent effects) of a treatment B in (un)balanced factorial experiments. The main user functions are \code{\link{doPRC}}, \code{\link{plotPRC}}, \code{\link{plotPRC2d}} and \code{\link{smoothPRC}}. The easiest example is the code demo \code{PRC_pyrifos}. This version includes smoothPRC (\code{\link{smoothPRC}} using \code{\link{set_smoothPRC_model}}). |
| Authors: | Cajo J.F ter Braak [aut, cre] (ORCID: <https://orcid.org/0000-0002-0414-8745>) |
| Maintainer: | Cajo J.F. ter Braak <[email protected]> |
| License: | GPL-3 | file LICENSE |
| Version: | 2.0.4 |
| Built: | 2026-10-14 18:52:07 UTC |
| Source: | https://github.com/CajoterBraak/PRC |
Permutation Test for Smooth PRC
## S3 method for class 'smoothPRC' anova(object, ..., permutations = 19, verbose = TRUE)## S3 method for class 'smoothPRC' anova(object, ..., permutations = 19, verbose = TRUE)
object |
an object from |
... |
unused. |
permutations |
a list of control values for the permutations as
returned by the function |
verbose |
logical for printing progress of the permutations |
The test takes account of additional covariates in fixed of
the LLMsolver model. The test is approximate without permutation of whole
time series.
The test is based on the first axis.
The data frame biomonitoring contains the design and species data from van den Brink et al.(2009). The data is
on biomonitoring of macroinvertebrate species in the Rhine and Meuse rivers, near the
towns of Kampen and Grave, respectively, from 1992 - 2000, sampled six times per year,
except in 1998 (4x), 1999 and 2000 (7x). The first example in the demo monitoring_EMA2009
uses reference coding of an reduncancy analysis without covariates
so as to emphasize the compositional difference from a reference. The second example is a more regular PRC, contrasting
the compositional changes at one site without those in another.
sample character identifier for the sample.
year numeric year of sampling
location location of sampling: factor with levels c("Grave","Kampen")
counts of 139 macroinvertebrate species
van den Brink, PJ, den Besten, PJ, bij de Vaate, A and ter Braak, CJF, 2009. Principal response curves technique for the analysis of multivariate biomonitoring time series. Environmental Monitoring and Assessment. 152, 271-281. doi:10.1007/s10661-008-0314-6
The dataframe CPTNUT_phytoplankton contains the design
and phytoplankton data from
the microcosm example of PRC from van den Brink & ter Braak (1998).
In this experiment from van Donk et al. (1995), there are three
treatments: control (con), nutrients (nut) and nutrients with
chlorpyrifos application (cpf-nut), each with
four cosms (replications) sampled 13 times (twice before nutrients where
added and six times before the additional
application of chorpyrifos).
The phytoplankton data are counts and
there are 26 species in total.
A data frame with 156 rows and 31 columns:
samples character identifier for the sample.
Time factor coding for the week of sampling with respect to first
application of nutrients.
Treatment treatment factor as in
van den Brink & ter Braak (1998).
The levels are con, nut and cpf-nut
for the four control cosms,
the four cosms that receive nutrients after the first
two sampling points, and the four cosms to which both
nutrients and chorpyrifos are applied, respectively.
Note that the cosms of the nutrient-chlorpyrifos treatment are coded
herein as nut in the weeks
before the first chlorpyrifos application, as they are treated
experiment-wise similarly in these weeks. But, note that
the samples from con and nut are coded differently on
pre-treatment times (although they
are similar experiment-wise),
so as to see the pre-treatment variability.
See also the factor con_nut_cpf.
Plot numeric identifier for the microcosm.
con_nut_cpf treatment factor. Before application of nutrients or
chlorpyrifos, the samples are all treated similarly,
but are coded differently (as often in PRC).
cell counts of 25 different phytoplankton taxa in 1L, and density for the 26th species, Volvox. For details, see van den Donk et al. (1995).
van Donk et al. (1995), van den Brink & ter Braak (1998)
van den Brink, PJ & ter Braak, CJF (1998). Multivariate analysis of stress on experimental ecosystems by principal response curves and similarity analysis. Aquatic ecology, 32, 163-178. doi:10.1023/A:1009944004756
van Donk E, Prins H, Voogd HM, Crum SJH and Brock TCM (1995) Effects of nutrient loading and insecticide application on the ecology of Elodea – dominated freshwater microcosms. I. Responses of plankton and zooplanktivorous insects. Arch Hydrobiol 133: 417–439. doi:10.1127/archiv-hydrobiol/133/1995/417
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 )doPRC( formula, scale = FALSE, referencelevel = NULL, rank = 2, flip = rep(FALSE, rank), scaling = "ms", method = "rda", data )
formula |
formula specifying how the response (species data) must modeled by predictors and covariates.
A typical |
scale |
Scales species to unit variance (like correlations) if |
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 |
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? |
scaling |
character |
method |
character rda or cca. Default: |
data |
Data frame containing the variables on the right hand side of the model formula. |
The function performs PRC using rda or cca with
the specified formula and applies PRC_scores to the result.
an object as a cca.object with additional entries as output of PRC_scores.
The assigned classes are c("PRC","rda","cca".
ter Braak (2023) Redundancy analysis includes analysis of variance-simultaneous component analysis (ASCA) and outperforms its extensions Chemometrics and Intelligent Laboratory Systems. doi:10.1016/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. doi: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. doi: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. doi: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. doi:10.1007/s10661-008-0314-6
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) mod_prc$focal_and_conditioning_factors$`focal factor` # as expected: "dose" print(mod_prc) # 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) # 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 # only the PRC curves: plotPRC(mod_prc, plot = 0, symbols_on_curves = TRUE) # with lines and a stronger threshold plotPRC(mod_prc, plot = "ditch", threshold = 10,width = c(4,1)) # 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) 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)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) mod_prc$focal_and_conditioning_factors$`focal factor` # as expected: "dose" print(mod_prc) # 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) # 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 # only the PRC curves: plotPRC(mod_prc, plot = 0, symbols_on_curves = TRUE) # with lines and a stronger threshold plotPRC(mod_prc, plot = "ditch", threshold = 10,width = c(4,1)) # 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) 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)
Returns an object with a fitted() method.
fit_full_model_columns(model, ...)fit_full_model_columns(model, ...)
model |
Model specification. |
... |
Additional arguments. |
Returns an object with a fitted() method.
fit_full_model_rows(model, ...)fit_full_model_rows(model, ...)
model |
Model specification. |
... |
Additional arguments. |
Returns an object with a fitted() method.
fit_null_model_columns(model, ...)fit_null_model_columns(model, ...)
model |
Model specification. |
... |
Additional arguments. |
Returns an object with a fitted() method.
fit_null_model_rows(model, ...)fit_null_model_rows(model, ...)
model |
Model specification. |
... |
Additional arguments. |
fvalues4levels changes levels of a factor to numeric variable
in set_smoothPRC_model and
for plotting of the PRC diagram.
fvalues4levels(object, fctr)fvalues4levels(object, fctr)
object |
a data frame or tibble. |
fctr |
name of a factor in object (character or symbol) |
a numeric named vector (names are levels)
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))) fvalues4levels(Design, "week") # fvalues4levels strips the starting k characters so as to find numbers,e.g. Design$week2 <- factor(paste("Week", Design$week)) levels(Design$week2) fvalues4levels(Design, "week2") # they may not be in increasing order in general, # but will be after sorting: sort(fvalues4levels(Design, "week2")) # the resulting numeric variable: numvar <- fvalues4levels(Design, "week2")[Design[,"week2"]] summary(numvar) length(numvar)== nrow(Design)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))) fvalues4levels(Design, "week") # fvalues4levels strips the starting k characters so as to find numbers,e.g. Design$week2 <- factor(paste("Week", Design$week)) levels(Design$week2) fvalues4levels(Design, "week2") # they may not be in increasing order in general, # but will be after sorting: sort(fvalues4levels(Design, "week2")) # the resulting numeric variable: numvar <- fvalues4levels(Design, "week2")[Design[,"week2"]] summary(numvar) length(numvar)== nrow(Design)
get_focal_and_conditioning_factors derives the focal constraining variables(s) and (if present) the
conditioning variables(s) from a rda
or cca ordination specified via a formula for use in a PRC
(with doPRC or PRC_scores)
get_focal_and_conditioning_factors(object, factors_only = TRUE)get_focal_and_conditioning_factors(object, factors_only = TRUE)
object |
a result of |
factors_only |
logical to exclude numeric variables from the focal factor(s). Default |
A list with element `focal factor` and condition
ter Braak C.J.F. & Šmilauer P. (2018): Canoco reference manual and user's guide: software for ordination, version 5.1x. Microcomputer Power, Ithaca, USA, 536 pp. (http::www.canoco5.com)
data(pyrifos, package = "vegan") #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))) library(vegan) mod_rda <- rda( pyrifos ~ dose:week + Condition(week), data = Design) species_scores= scores(mod_rda,choices =1:2, display="sp",scaling = 1) library(PRC) head(rdacca_taxon_stats(mod_rda)) # stats per species get_focal_and_conditioning_factors(mod_rda) # factors from the model species_scores= vegan::scores(mod_rda,choices =1:2, display="sp",scaling = 2) # 1-d vertical plot of species scores (a cut of RDA1-scores at 0.5) plot_species_scores_bk(species_scores) # a cut of RDA1-scores at 1.0 (as scoresname is unset) plot_species_scores_bk(species_scores, threshold = 14) plot_species_scores_2d(species_scores) plot_species_scores_2d(species_scores, label.repel = TRUE) plot_species_scores_2d(species_scores, max.overlaps = 30, label.repel = TRUE)data(pyrifos, package = "vegan") #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))) library(vegan) mod_rda <- rda( pyrifos ~ dose:week + Condition(week), data = Design) species_scores= scores(mod_rda,choices =1:2, display="sp",scaling = 1) library(PRC) head(rdacca_taxon_stats(mod_rda)) # stats per species get_focal_and_conditioning_factors(mod_rda) # factors from the model species_scores= vegan::scores(mod_rda,choices =1:2, display="sp",scaling = 2) # 1-d vertical plot of species scores (a cut of RDA1-scores at 0.5) plot_species_scores_bk(species_scores) # a cut of RDA1-scores at 1.0 (as scoresname is unset) plot_species_scores_bk(species_scores, threshold = 14) plot_species_scores_2d(species_scores) plot_species_scores_2d(species_scores, label.repel = TRUE) plot_species_scores_2d(species_scores, max.overlaps = 30, label.repel = TRUE)
initializes the power algorthm (internal)
initialize_pwr_alg(model, Y, weights, axes)initialize_pwr_alg(model, Y, weights, axes)
model |
Model specification. |
Y |
response |
weights |
optional weights |
axes |
prior axes and covariable data |
The dataframe Ossenkampen contains the data from
a long-term fertilizer experiment (1958-2007) in grassland with 98 response variables
(abundances of the plant species that grew in the Ossenkampen plots, measured as counts in 100 subsamples),
as analysed in ter Braak (2023).
The experiment started in 1958 with 12 plots arranged in two randomized complete blocks with six treatments: four types of fertilizer (K, P, PK or NPK), liming and no fertilizer. Eight years later (in 1966) each block was extended with two plots, one with N and another with NPK fertilization. The two limed plots behaved rather differently from the rest and are omitted here. In this data, the K and P plots are combined with the no fertilizer plots so that the control treatment (Cntrl) consisted of 6 plots, with the aim to make the variability among replicates more interesting. Three plots of the control and one PK plot were not sampled in 1967 and two NPK samples were not sampled in 1984; these missing samples were not imputed in this data set.
The data analyzed in Berendse et al. (2021) also contained the limed plots and contained, for the multivariate analyses, imputations of the few missing samples after 1966.
Sample numeric identifier for a sample
Plot identifier for a sample (the identifier details the fertilizer treatment and Block)
Block numeric identifier of the block (1 or 2)
A year of sampling
B fertilizer treatment, condensed to four treatments ("Cntrl" is the referencelevel)
AgrCap count out of 100 subsamples of species AgrCap, and so on for the remaining 97 species
Berendse, F., Geerts, R.H.E.M., Elberse, W.T., Bezemer, T.M., Goedhart, P.W., Xue, W., Noordijk, E., ter Braak, C.J.F. & Korevaar, H. (2021) A matter of time: Recovery of plant species diversity in wild plant communities at declining nitrogen deposition. Diversity and Distributions, 27, 1180-1193. doi:10.1111/ddi.13266
Geerts, R.H.E.M., Bufe, C., Schils, P.C., Korevaar, H., Berendse, F. & Struik, P.C. (2021) The Ossekampen Grassland Experiment: Data underlying the publication: A matter of time: recovery of plant species diversity in wild plant communities at declining nitrogen deposition. doi:10.4121/13200398
ter Braak, C.J.F. (2023) Redundancy analysis includes analysis of variance-simultaneous component analysis (ASCA) and outperforms its extensions. Chemometrics and Intelligent Laboratory Systems, 104898. doi:10.1016/j.chemolab.2023.104898
plot_sample_scores_cdt creates a PRC diagram of the treatments without species loadings.
plot_sample_scores_cdt( object, treatment = NULL, condition = NULL, plot = 1, xvals = NULL, score_name = "PRC1", size = c(2, 2), symbols_on_curves = FALSE, verbose = TRUE )plot_sample_scores_cdt( object, treatment = NULL, condition = NULL, plot = 1, xvals = NULL, score_name = "PRC1", size = c(2, 2), symbols_on_curves = FALSE, verbose = TRUE )
object |
a result of |
treatment |
character vector of names of factors in |
condition |
character vector of names of factors in |
plot |
0 or 1 (no/yes individual plot points) or
character name of a factor in |
xvals |
numeric values for the levels of the condition factor,
or name of numeric variable in |
score_name |
name of scores to use as y-axis; default: "PRC1". |
size |
size of symbols and symbols in the legend. default: c(2, 2). If the length is not 2 than the first entry is replicated. |
symbols_on_curves |
logical, draw symbols on the PRC curves at the measurement times (default FALSE).
If |
verbose |
logical for printing the current model, |
plot_sample_scores_cdt can be used to
plot any ordination score or variable against the condition with treatments indicated by colors and line type
by setting score_name to the name of the variable of interest and argument plot to 0.
If object is not a result of doPRC or PRC_scores), set treatment and condition
(or object$'focal_and_conditioning_factors') and set object to a list with element PRCplus;
object$PRCplus should contain the names indicated by the arguments treatment, condition and, optionally, plot
treatment: a factor for the treatments,
condition: a factor or variable condition
If plot!=0, the scores (values) for individual sample points or lines must be in the variable "PRC1_E" or,
in general, paste(score_name,"_E", sep = "").
a ggplot object
PRC_scores, get_focal_and_conditioning_factors, plotPRC and fvalues4levels
data(pyrifos, package = "vegan") #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_rda <- vegan::rda( pyrifos ~ dose:week + Condition(week), data = Design) Design_w_PRCs <- PRC_scores(mod_rda, focal_factor_name= "dose", referencelevel = "0", rank = 2, data= Design, scaling ="ms", flip = FALSE) plot_sample_scores_cdt(Design_w_PRCs, plot = "ditch") plot_sample_scores_cdt(Design_w_PRCs) plot_sample_scores_cdt(Design_w_PRCs, plot = "Ditch")data(pyrifos, package = "vegan") #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_rda <- vegan::rda( pyrifos ~ dose:week + Condition(week), data = Design) Design_w_PRCs <- PRC_scores(mod_rda, focal_factor_name= "dose", referencelevel = "0", rank = 2, data= Design, scaling ="ms", flip = FALSE) plot_sample_scores_cdt(Design_w_PRCs, plot = "ditch") plot_sample_scores_cdt(Design_w_PRCs) plot_sample_scores_cdt(Design_w_PRCs, plot = "Ditch")
plot_smoothPRC_cdt creates a PRC diagram of the treatments without species loadings.
plot_smoothPRC_cdt( object, axis = 1, mod_prc = NULL, with_classical_lines = TRUE, with_unconstrained_scores = TRUE, flip = FALSE, sample_times_only = FALSE )plot_smoothPRC_cdt( object, axis = 1, mod_prc = NULL, with_classical_lines = TRUE, with_unconstrained_scores = TRUE, flip = FALSE, sample_times_only = FALSE )
object |
a result of |
axis |
the number of the axis to plot |
mod_prc |
a results of |
with_classical_lines |
logical to add or suppress classical PRC lines.
Default |
with_unconstrained_scores |
logical to add or suppress the unconstrained
scores in the smooth PRC showing the scatter of the
individual samples (cosms).
Default |
flip |
logical. Should the axis be reversed?
Default |
sample_times_only |
logical to PRClines based on the sampled times
only, i.e. without interpolation. Default |
# Analyze example data ---------------------------------------------------- data(SimData) # extract design Design0 <- SimData[,c("A","B")] Design0$A <- factor(Design0$A); Design0$B <- factor(Design0$B) #put levels in natural order Design0$A<- factor(Design0$A, levels=c(levels(Design0$A)[-2],levels(Design0$A)[2])) Y0 <- as.matrix(SimData[,-(1:3)]) ids <- which(Design0$A %in% "A10" & Design0$B %in% "B5") #design <- "Complete" # design <- "Empty cell" if (design == "Complete"){ Design <- Design0;Y <- Y0} else {Design <-Design0[-ids,]; Y=Y0[-ids,]} print(with(Design,table(A,B))) names(Design) # Design and Y for Fig 1 of ter Braak 2023 ready ------------------------------------------- dim(Y) str(Design) names(Design) <- c("Time","Treatment") mod_prc <- doPRC(Y ~ Treatment:Time + Condition(Time), data = Design) #Mod_prc <- vegan::prc(Y, Design$Treatment,Design$Time) #Mod_prc$CCA$eig# identical to doPRC b_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "sp") x_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "lc") summary(mod_prc) smoothPRC_model <- set_smoothPRC_model(data = Design) smoothPRC_model$spline smooth_PRC <- smoothPRC(Y, model = smoothPRC_model) summary(smooth_PRC$obj[[1]]) (cc <- cor(smooth_PRC$B[,1], b_mod_prc)) plot_species_scores_bk(-smooth_PRC$B, scoresname = "B1") plot_species_scores_bk(b_mod_prc) plot_smoothPRC_cdt(smooth_PRC, flip= TRUE) plotsmoothPRC(smooth_PRC, flip= TRUE,threshold = 0) # Demonstration of the fitted model x1 <- scale(as.numeric(Design$Time))[,1] x2 <- scale(as.numeric(Design$Treatment))[,1] # with the scaling to mean square of 1 of the species scores b summary(lm(smooth_PRC$X[,1]~ x1+ x2+ x1:x2 )) # R^2 ~1! plotPRC(mod_prc)# Analyze example data ---------------------------------------------------- data(SimData) # extract design Design0 <- SimData[,c("A","B")] Design0$A <- factor(Design0$A); Design0$B <- factor(Design0$B) #put levels in natural order Design0$A<- factor(Design0$A, levels=c(levels(Design0$A)[-2],levels(Design0$A)[2])) Y0 <- as.matrix(SimData[,-(1:3)]) ids <- which(Design0$A %in% "A10" & Design0$B %in% "B5") #design <- "Complete" # design <- "Empty cell" if (design == "Complete"){ Design <- Design0;Y <- Y0} else {Design <-Design0[-ids,]; Y=Y0[-ids,]} print(with(Design,table(A,B))) names(Design) # Design and Y for Fig 1 of ter Braak 2023 ready ------------------------------------------- dim(Y) str(Design) names(Design) <- c("Time","Treatment") mod_prc <- doPRC(Y ~ Treatment:Time + Condition(Time), data = Design) #Mod_prc <- vegan::prc(Y, Design$Treatment,Design$Time) #Mod_prc$CCA$eig# identical to doPRC b_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "sp") x_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "lc") summary(mod_prc) smoothPRC_model <- set_smoothPRC_model(data = Design) smoothPRC_model$spline smooth_PRC <- smoothPRC(Y, model = smoothPRC_model) summary(smooth_PRC$obj[[1]]) (cc <- cor(smooth_PRC$B[,1], b_mod_prc)) plot_species_scores_bk(-smooth_PRC$B, scoresname = "B1") plot_species_scores_bk(b_mod_prc) plot_smoothPRC_cdt(smooth_PRC, flip= TRUE) plotsmoothPRC(smooth_PRC, flip= TRUE,threshold = 0) # Demonstration of the fitted model x1 <- scale(as.numeric(Design$Time))[,1] x2 <- scale(as.numeric(Design$Treatment))[,1] # with the scaling to mean square of 1 of the species scores b summary(lm(smooth_PRC$X[,1]~ x1+ x2+ x1:x2 )) # R^2 ~1! plotPRC(mod_prc)
plot_species_scores_2d creates an ordination diagram of species with selection criterion
of which species to plot with names.
plot_species_scores_2d( species_scores, threshold = 7, speciesname = NULL, scorenames = c("RDA1", "RDA2"), selectname = "Fratio2", label.repel = FALSE, withnames_only = FALSE, max.overlaps = 10, mult_expand = 0.1, verbose = TRUE )plot_species_scores_2d( species_scores, threshold = 7, speciesname = NULL, scorenames = c("RDA1", "RDA2"), selectname = "Fratio2", label.repel = FALSE, withnames_only = FALSE, max.overlaps = 10, mult_expand = 0.1, verbose = TRUE )
species_scores |
a species-by-scores matrix, a data frame with rownames (species names) or a tibble
with variable with name |
threshold |
species with criterion (specified by |
speciesname |
name of the variable containing the species names (default |
scorenames |
names of the two columns or variables containing the species scores to be plotted
(default |
selectname |
name of the column or variable containing the criterion for the selection of species to be displayed
Default: |
label.repel |
logical (default: |
withnames_only |
logical to plot the species that exceed the |
max.overlaps |
exclude text labels that overlap too many things. Default: 10 (see |
mult_expand |
fraction of range expansion of the diagram (default = 0.1). See |
verbose |
logical for printing the number of species
that |
The names of species are added using geom_text_repel or geom_label_repel.
These functions may not be able to plot the names of all selected species.
This issue can be diminished by increasing the value of max.overlaps.
a ggplot object
data(pyrifos, package = "vegan") #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))) library(vegan) mod_rda <- rda( pyrifos ~ dose:week + Condition(week), data = Design) species_scores= scores(mod_rda,choices =1:2, display="sp",scaling = 1) library(PRC) head(rdacca_taxon_stats(mod_rda)) # stats per species get_focal_and_conditioning_factors(mod_rda) # factors from the model species_scores= vegan::scores(mod_rda,choices =1:2, display="sp",scaling = 2) # 1-d vertical plot of species scores (a cut of RDA1-scores at 0.5) plot_species_scores_bk(species_scores) # a cut of RDA1-scores at 1.0 (as scoresname is unset) plot_species_scores_bk(species_scores, threshold = 14) plot_species_scores_2d(species_scores) plot_species_scores_2d(species_scores, label.repel = TRUE) plot_species_scores_2d(species_scores, max.overlaps = 30, label.repel = TRUE)data(pyrifos, package = "vegan") #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))) library(vegan) mod_rda <- rda( pyrifos ~ dose:week + Condition(week), data = Design) species_scores= scores(mod_rda,choices =1:2, display="sp",scaling = 1) library(PRC) head(rdacca_taxon_stats(mod_rda)) # stats per species get_focal_and_conditioning_factors(mod_rda) # factors from the model species_scores= vegan::scores(mod_rda,choices =1:2, display="sp",scaling = 2) # 1-d vertical plot of species scores (a cut of RDA1-scores at 0.5) plot_species_scores_bk(species_scores) # a cut of RDA1-scores at 1.0 (as scoresname is unset) plot_species_scores_bk(species_scores, threshold = 14) plot_species_scores_2d(species_scores) plot_species_scores_2d(species_scores, label.repel = TRUE) plot_species_scores_2d(species_scores, max.overlaps = 30, label.repel = TRUE)
plot_species_scores_bk creates a vertical line plot of ordination
scores with selection criterion for which scores to plot with names.
plot_species_scores_bk( species_scores, ylab = latex2exp::TeX("$b_k$"), threshold = 7, y_lab_interval = 0.5, speciesname = NULL, scoresname = "RDA1", selectname = "Fratio1", speciesgroup = NULL, expand = 0.2, verbose = TRUE )plot_species_scores_bk( species_scores, ylab = latex2exp::TeX("$b_k$"), threshold = 7, y_lab_interval = 0.5, speciesname = NULL, scoresname = "RDA1", selectname = "Fratio1", speciesgroup = NULL, expand = 0.2, verbose = TRUE )
species_scores |
a species-by-scores matrix, a data frame with
row names (species names) or a tibble with variable with name
|
ylab |
y-axis label. Default: $b_k$. |
threshold |
species with criterion (specified by |
y_lab_interval |
interval of the y-axis ticks. A tick at no effect (0) is always included; default: 0.5. |
speciesname |
name of the variable containing the species names
(default |
scoresname |
name of the column or variable containing the species
scores to be plotted (default |
selectname |
name of the column or variable containing the criterion
for the selection of species to be displayed. Default: |
speciesgroup |
name of the factor, the levels of which receive different colors in the vertical plot. |
expand |
amount of extension of the line plot (default 0.2). |
verbose |
logical for printing the number of species with names out of
the total number (default: |
The absolute value of the criterion values is taken before selection, so
that also the species scores of the ordination can serve as a criterion
(e.g. selectname = "RDA1"). The function has been copied from
the douconca package at
https://cran.r-project.org/web/packages/douconca.
The function is used in plotPRC
and plotsmoothPRC.
a ggplot object
# Analyze example data ---------------------------------------------------- data(SimData) # extract design Design0 <- SimData[,c("A","B")] Design0$A <- factor(Design0$A); Design0$B <- factor(Design0$B) #put levels in natural order Design0$A<- factor(Design0$A, levels=c(levels(Design0$A)[-2],levels(Design0$A)[2])) Y0 <- as.matrix(SimData[,-(1:3)]) ids <- which(Design0$A %in% "A10" & Design0$B %in% "B5") #design <- "Complete" # design <- "Empty cell" if (design == "Complete"){ Design <- Design0;Y <- Y0} else {Design <-Design0[-ids,]; Y=Y0[-ids,]} print(with(Design,table(A,B))) names(Design) # Design and Y for Fig 1 of ter Braak 2023 ready ------------------------------------------- dim(Y) str(Design) names(Design) <- c("Time","Treatment") mod_prc <- doPRC(Y ~ Treatment:Time + Condition(Time), data = Design) #Mod_prc <- vegan::prc(Y, Design$Treatment,Design$Time) #Mod_prc$CCA$eig# identical to doPRC b_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "sp") x_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "lc") summary(mod_prc) smoothPRC_model <- set_smoothPRC_model(data = Design) smoothPRC_model$spline smooth_PRC <- smoothPRC(Y, model = smoothPRC_model) summary(smooth_PRC$obj[[1]]) (cc <- cor(smooth_PRC$B[,1], b_mod_prc)) plot_species_scores_bk(-smooth_PRC$B, scoresname = "B1") plot_species_scores_bk(b_mod_prc) plot_smoothPRC_cdt(smooth_PRC, flip= TRUE) plotsmoothPRC(smooth_PRC, flip= TRUE,threshold = 0) # Demonstration of the fitted model x1 <- scale(as.numeric(Design$Time))[,1] x2 <- scale(as.numeric(Design$Treatment))[,1] # with the scaling to mean square of 1 of the species scores b summary(lm(smooth_PRC$X[,1]~ x1+ x2+ x1:x2 )) # R^2 ~1! plotPRC(mod_prc)# Analyze example data ---------------------------------------------------- data(SimData) # extract design Design0 <- SimData[,c("A","B")] Design0$A <- factor(Design0$A); Design0$B <- factor(Design0$B) #put levels in natural order Design0$A<- factor(Design0$A, levels=c(levels(Design0$A)[-2],levels(Design0$A)[2])) Y0 <- as.matrix(SimData[,-(1:3)]) ids <- which(Design0$A %in% "A10" & Design0$B %in% "B5") #design <- "Complete" # design <- "Empty cell" if (design == "Complete"){ Design <- Design0;Y <- Y0} else {Design <-Design0[-ids,]; Y=Y0[-ids,]} print(with(Design,table(A,B))) names(Design) # Design and Y for Fig 1 of ter Braak 2023 ready ------------------------------------------- dim(Y) str(Design) names(Design) <- c("Time","Treatment") mod_prc <- doPRC(Y ~ Treatment:Time + Condition(Time), data = Design) #Mod_prc <- vegan::prc(Y, Design$Treatment,Design$Time) #Mod_prc$CCA$eig# identical to doPRC b_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "sp") x_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "lc") summary(mod_prc) smoothPRC_model <- set_smoothPRC_model(data = Design) smoothPRC_model$spline smooth_PRC <- smoothPRC(Y, model = smoothPRC_model) summary(smooth_PRC$obj[[1]]) (cc <- cor(smooth_PRC$B[,1], b_mod_prc)) plot_species_scores_bk(-smooth_PRC$B, scoresname = "B1") plot_species_scores_bk(b_mod_prc) plot_smoothPRC_cdt(smooth_PRC, flip= TRUE) plotsmoothPRC(smooth_PRC, flip= TRUE,threshold = 0) # Demonstration of the fitted model x1 <- scale(as.numeric(Design$Time))[,1] x2 <- scale(as.numeric(Design$Treatment))[,1] # with the scaling to mean square of 1 of the species scores b summary(lm(smooth_PRC$X[,1]~ x1+ x2+ x1:x2 )) # R^2 ~1! plotPRC(mod_prc)
plotPRC creates a one-dimensional PRC diagram with
vertical line plot of species loadings from a result of doPRC or PRC_scores.
plotPRC( object, treatment = NULL, condition = NULL, plot = 1, xvals = NULL, axis = 1, size = c(2, 2), symbols_on_curves = FALSE, widths = c(4, 1), title = NULL, left = "Treatment and plot scores", right = "taxon scores", threshold = 7, y_lab_interval = 0.5, speciesname = NULL, selectname = "Fratio", verbose = TRUE )plotPRC( object, treatment = NULL, condition = NULL, plot = 1, xvals = NULL, axis = 1, size = c(2, 2), symbols_on_curves = FALSE, widths = c(4, 1), title = NULL, left = "Treatment and plot scores", right = "taxon scores", threshold = 7, y_lab_interval = 0.5, speciesname = NULL, selectname = "Fratio", verbose = TRUE )
object |
a result of |
treatment |
character vector of names of factors in |
condition |
character vector of names of factors in |
plot |
0 or 1 (no/yes individual plot points) or
character name of a factor in |
xvals |
numeric values for the levels of the condition factor,
or name of numeric variable in |
axis |
axis shown; default 1, showing the first PRC axis. |
size |
size of symbols and symbols in the legend. default: c(2, 2). If the length is not 2 than the first entry is replicated. |
symbols_on_curves |
logical, draw symbols on the PRC curves at the measurement times (default FALSE).
If |
widths |
relative width of treatment PRC plot and the loading plot (see |
title |
overall title. |
left |
left axis text (see |
right |
right axis text (see |
threshold |
species with criterion (specified by |
y_lab_interval |
interval of the y-axis ticks. A tick at no effect (0) is always included; default: 0.5. |
speciesname |
name of the variable containing the species names (default |
selectname |
name of the column or variable containing the criterion for the selection of species to be displayed
without |
verbose |
logical for plotting and printing the current model, |
The function creates the treatment PRC diagram using plot_sample_scores_cdt and the
vertical line plot of loadings using plot_species_scores_bk. These ggplot2-objects are
combined in a single object using arrangeGrob.
a list of a gridExtra object (plot using arrangeGrob)
and a list of two ggplot2 objects (the treatment PRC-diagram and the vertical line plot).
doPRC, plotPRC2d, plot_sample_scores_cdt,
plot_species_scores_bk and grid.arrange
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) mod_prc$focal_and_conditioning_factors$`focal factor` # as expected: "dose" print(mod_prc) # 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) # 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 # only the PRC curves: plotPRC(mod_prc, plot = 0, symbols_on_curves = TRUE) # with lines and a stronger threshold plotPRC(mod_prc, plot = "ditch", threshold = 10,width = c(4,1)) # 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) 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)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) mod_prc$focal_and_conditioning_factors$`focal factor` # as expected: "dose" print(mod_prc) # 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) # 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 # only the PRC curves: plotPRC(mod_prc, plot = 0, symbols_on_curves = TRUE) # with lines and a stronger threshold plotPRC(mod_prc, plot = "ditch", threshold = 10,width = c(4,1)) # 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) 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)
plotPRC2d creates a two-dimensional PRC diagram with species loadings from a result of doPRC or PRC_scores
plotPRC2d( object, treatment = NULL, condition = NULL, plot = 1, xvals = NULL, axes = c(1, 2), threshold = 7, title = NULL, left = NULL, right = NULL, speciesname = NULL, selectname = "Fratio", xylab = NULL, label.repel = FALSE, withnames_only = FALSE, max.overlaps = 10, mult_expand = 0.1, widths = c(3, 2), verbose = TRUE )plotPRC2d( object, treatment = NULL, condition = NULL, plot = 1, xvals = NULL, axes = c(1, 2), threshold = 7, title = NULL, left = NULL, right = NULL, speciesname = NULL, selectname = "Fratio", xylab = NULL, label.repel = FALSE, withnames_only = FALSE, max.overlaps = 10, mult_expand = 0.1, widths = c(3, 2), verbose = TRUE )
object |
a result of |
treatment |
character vector of names of factors in |
condition |
character vector of names of factors in |
plot |
0 or 1 (no/yes individual plot points) or
character name of a factor in |
xvals |
numeric values for the levels of the condition factor,
or name of numeric variable in |
axes |
axes shown; default |
threshold |
species with criterion (specified by |
title |
overall title. |
left |
left axis text (see |
right |
right axis text (see |
speciesname |
name of the variable containing the species names (default |
selectname |
name of the column or variable containing the criterion for the selection of species to be displayed
Default: |
xylab |
character vector of length 2 defining x- and y-axis title of the species ordination so as to overwrite the defauls axis labels |
label.repel |
logical (default: |
withnames_only |
logical to plot the species that exceed the |
max.overlaps |
exclude text labels that overlap too many things. Default: 10 (see |
mult_expand |
fraction of range expansion of the diagram (default = 0.1). See |
widths |
relative width of treatment PRC plot and the loading plot (see |
verbose |
logical for printing the number of species
that |
The function creates the treatment PRC diagrams for the axes using plot_sample_scores_cdt, with diagrams for the
x=y and x=-y directions, and the
2d species loading plot using plot_species_scores_2d. These ggplot2-objects are
combined in a single object using arrangeGrob.
a list of a gridExtra object (plot using arrangeGrob)
and a list of five ggplot2 objects, arranged in a list of two: treatments = the PRC1, PRC2-diagrams with line x=y and x= -y diagrams
and and species= the 2d species ordination plot).
doPRC, plot_sample_scores_cdt,
plot_species_scores_2d and grid.arrange
library(PRC) data("CPFNUTphytoplankton") # Read example data ---------------------------------------------------- names(CPFNUTphytoplankton)[1:6] Design <- CPFNUTphytoplankton[,c("Time","Treatment","Plot","con_nut_cpf")] summary(Design) print(with(Design,table(Time,Treatment))) # extract response variables Counts <- as.matrix(CPFNUTphytoplankton[,-(1:5)]); rownames(Counts)<-CPFNUTphytoplankton$sample # transform as a multiplicative model is more logical than an additive one. # But need to add pseudocount to avoid log(0) = infinite Y <- log(Counts+1) # Remark: this transformation reproduces the results in # van den Brink & ter Braak 1998 # apart from scaling (but the product b_k*c_dt is identical) and # apart from the PRC1 score at week 4, which is erroneous in the paper. # PRC via vegan ----------------------------------------------------------- mod_prc <- doPRC( Y~ Time * Treatment + Condition(Time), rank = 3, data= Design, flip= c(FALSE, TRUE, FALSE)) round(mod_prc$percExp, 0) plotPRC(mod_prc, plot = "Plot", width = c(2,1), symbols_on_curves = TRUE) plotPRC(mod_prc, plot = "Plot", width = c(2,1), axis = 2) # note that the default selection criterion "Fratio2" also selects species #strongly reacting to axis 1 plotPRC2d(mod_prc, plot = 0, threshold = 1, title = "PRC Phytoplankton van den Brink & ter Braak (1998)") # modifying the plot # Assign the plots to symbols and use grid.arrange to produce the plot you like, #for example: gg <- plotPRC2d(mod_prc, plot = 0, verbose = FALSE) pl_species <- gg$separateplots$species # species ordination pl <- gg$separateplots$treatments names(pl) layout<- rbind(c(1,2), c(3,4),c(3,5)) gridExtra::grid.arrange(pl[["PRC1"]]+ggplot2::ylab("PRC1")+ggplot2::ggtitle(""), pl[["PRC2"]]+ggplot2::ylab("PRC2")+ggplot2::xlab("")+ ggplot2::theme(legend.position="none")+ggplot2::ggtitle(""), pl_species, pl[["line11"]]+ggplot2::ylab("scores x = y")+ ggplot2::xlab("")+ggplot2::theme(legend.position="none")+ ggplot2::ggtitle(""), pl[["line1-1"]]+ggplot2::ylab("scores x = -y")+ ggplot2::theme(legend.position="none")+ ggplot2::ggtitle(""), layout_matrix = layout, widths = c(3,2), top = "mytitle", left ="", right = "")library(PRC) data("CPFNUTphytoplankton") # Read example data ---------------------------------------------------- names(CPFNUTphytoplankton)[1:6] Design <- CPFNUTphytoplankton[,c("Time","Treatment","Plot","con_nut_cpf")] summary(Design) print(with(Design,table(Time,Treatment))) # extract response variables Counts <- as.matrix(CPFNUTphytoplankton[,-(1:5)]); rownames(Counts)<-CPFNUTphytoplankton$sample # transform as a multiplicative model is more logical than an additive one. # But need to add pseudocount to avoid log(0) = infinite Y <- log(Counts+1) # Remark: this transformation reproduces the results in # van den Brink & ter Braak 1998 # apart from scaling (but the product b_k*c_dt is identical) and # apart from the PRC1 score at week 4, which is erroneous in the paper. # PRC via vegan ----------------------------------------------------------- mod_prc <- doPRC( Y~ Time * Treatment + Condition(Time), rank = 3, data= Design, flip= c(FALSE, TRUE, FALSE)) round(mod_prc$percExp, 0) plotPRC(mod_prc, plot = "Plot", width = c(2,1), symbols_on_curves = TRUE) plotPRC(mod_prc, plot = "Plot", width = c(2,1), axis = 2) # note that the default selection criterion "Fratio2" also selects species #strongly reacting to axis 1 plotPRC2d(mod_prc, plot = 0, threshold = 1, title = "PRC Phytoplankton van den Brink & ter Braak (1998)") # modifying the plot # Assign the plots to symbols and use grid.arrange to produce the plot you like, #for example: gg <- plotPRC2d(mod_prc, plot = 0, verbose = FALSE) pl_species <- gg$separateplots$species # species ordination pl <- gg$separateplots$treatments names(pl) layout<- rbind(c(1,2), c(3,4),c(3,5)) gridExtra::grid.arrange(pl[["PRC1"]]+ggplot2::ylab("PRC1")+ggplot2::ggtitle(""), pl[["PRC2"]]+ggplot2::ylab("PRC2")+ggplot2::xlab("")+ ggplot2::theme(legend.position="none")+ggplot2::ggtitle(""), pl_species, pl[["line11"]]+ggplot2::ylab("scores x = y")+ ggplot2::xlab("")+ggplot2::theme(legend.position="none")+ ggplot2::ggtitle(""), pl[["line1-1"]]+ggplot2::ylab("scores x = -y")+ ggplot2::theme(legend.position="none")+ ggplot2::ggtitle(""), layout_matrix = layout, widths = c(3,2), top = "mytitle", left ="", right = "")
Creates a smooth PRC diagram with treatment curves and a vertical species-loading plot.
plotsmoothPRC( object, axis = 1, mod_prc = NULL, with_classical_lines = TRUE, with_unconstrained_scores = TRUE, flip = FALSE, widths = c(4, 1), title = NULL, left = "Treatment curves", right = "taxon scores", threshold = 7, y_lab_interval = 0.5, speciesname = NULL, selectname = "Fratio", sample_times_only = FALSE, verbose = TRUE )plotsmoothPRC( object, axis = 1, mod_prc = NULL, with_classical_lines = TRUE, with_unconstrained_scores = TRUE, flip = FALSE, widths = c(4, 1), title = NULL, left = "Treatment curves", right = "taxon scores", threshold = 7, y_lab_interval = 0.5, speciesname = NULL, selectname = "Fratio", sample_times_only = FALSE, verbose = TRUE )
object |
a result of |
axis |
the number of the axis to plot |
mod_prc |
a results of |
with_classical_lines |
logical to add or suppress classical PRC lines.
Default |
with_unconstrained_scores |
logical to add or suppress the unconstrained
scores in the smooth PRC showing the scatter of the
individual samples (cosms).
Default |
flip |
logical. Should the axis be reversed?
Default |
widths |
relative width of treatment PRC plot and the loading plot (see |
title |
overall title. |
left |
left axis text (see |
right |
right axis text (see |
threshold |
species with criterion (specified by |
y_lab_interval |
interval of the y-axis ticks. A tick at no effect (0) is always included; default: 0.5. |
speciesname |
name of the variable containing the species names
(default |
selectname |
name of the column or variable containing the criterion
for the selection of species to be displayed. Default: |
sample_times_only |
logical to PRClines based on the sampled times
only, i.e. without interpolation. Default |
verbose |
logical for printing the number of species with names out of
the total number (default: |
sample_times_only = TRUE can be useful for comparison with
classical PRC. It avoids interpolation and may work when
the default fails to produce a plot.
PRC_scores calculates scores to display principal response curves with
scores for individual plots from a vegan constrained ordination object (cca.object).
PRC_scores( object, focal_factor_name = NULL, referencelevel = NULL, rank = 2, flip = rep(FALSE, rank), scaling = "ms", data )PRC_scores( object, focal_factor_name = NULL, referencelevel = NULL, rank = 2, flip = rep(FALSE, rank), scaling = "ms", data )
object |
a |
focal_factor_name |
name of one or two focal factors (the treatment(s) of interest); if unspecified,
derived from the formula in |
referencelevel |
numeric or character level(s) to be used as reference of the |
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? |
scaling |
character |
data |
data frame of the experimental design used in the call
to |
PRC_scores uses scores to
calculate the loadings (species scores),if obtainable from the object, and the constrained and unconstrained sample scores
of the units in data. The PRC treatment and sample scores are obtained by subtracting the reference scores
from the constrained and unconstrained scores, respectively.
The reference scores are calculated by applying predict
with type = "lc" on newdata, which
equals data, except that the level of the focal factor of each unit is changed to the reference level.
By this simple method, prc scores can be applied with and without terms in Condition()
in the formula of the call to rda, cca,
dbrda or capscale.
There can thus be additional (perhaps quantitative) covariates,
A:B+ Condition(A+other_covariates) where the focal factor is B,
and also with three (instead of two factors) e.g.
with formula A*B*C + Condition(A*C) so as to display the A*C-dependent effects of focal factor B
on the multivariate response in the call to rda, cca,
dbrda or capscale.
The element object$'focal_and_conditioning_factors' in value is set
by applying the function get_focal_and_conditioning_factors
to the ordination model specified in object.
CCA-based PRC is useful in data with many zeroes and which are strictly compositional, i.e. have row sums that have no meaning for the question of interest (ter Braak and te Beest, 2022). A layman's example is that the weight of a cake is not important for its taste; taste depends on its composition only. CCA can be applied without log-transforming the data and thereby avoids the issue of replacement of zeroes (ter Braak and te Beest, 2022).
A named list:
PRCplus: a data frame that combines the design in data with
the PRC and RDA/CCA scores
with (_E) and without error (the unconstrained and constrained
scores (sites and LC), respectively),
and, finally, the reference scores. Note that
PRC1_E = RDA1_E - RDA1_Ref, etc for other axes
species: a matrix with the loadings with respect to the PRC scores;
NULL if the object lacks loadings.
reference_scores: a matrix with the scores of each unit in
data with the level of the focal factor
replaced by the reference level
focal_factor_name: name of the focal factor in the call to
PRC_scores
referencelevel: reference level in the call to PRC_scores
coefficients: list of PRC-coefficients for each axis (matrix if rank is equal to 1)
focal_and_conditioning_factors: list of focal and
conditioning factors; obtained from
get_focal_and_conditioning_factors
if argument focal_factor_name in the call is unset.
terms: object$terms. See cca.object.
terminfo: object$terminfo. See cca.object.
method: object$method. See cca.object.
percExp: percentage variance of each axis out of the variance that can be explained by the predictors and covariates.
data(pyrifos, package = "vegan") #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_rda <- vegan::rda( pyrifos ~ dose:week + Condition(week), data = Design) Design_w_PRCs <- PRC_scores(mod_rda, focal_factor_name= "dose", referencelevel = "0", rank = 2, data= Design, scaling ="ms", flip = FALSE) plot_sample_scores_cdt(Design_w_PRCs, plot = "ditch") plot_sample_scores_cdt(Design_w_PRCs) plot_sample_scores_cdt(Design_w_PRCs, plot = "Ditch")data(pyrifos, package = "vegan") #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_rda <- vegan::rda( pyrifos ~ dose:week + Condition(week), data = Design) Design_w_PRCs <- PRC_scores(mod_rda, focal_factor_name= "dose", referencelevel = "0", rank = 2, data= Design, scaling ="ms", flip = FALSE) plot_sample_scores_cdt(Design_w_PRCs, plot = "ditch") plot_sample_scores_cdt(Design_w_PRCs) plot_sample_scores_cdt(Design_w_PRCs, plot = "Ditch")
rdacca_taxon_stats calculates fit statistics for taxa in an rda
or cca ordination.
rdacca_taxon_stats(object, rank = 1)rdacca_taxon_stats(object, rank = 1)
object |
|
rank |
number of axes for which fits are calculated (1:rank). Default: 1. |
The code uses ordiYbar and scores and also works for unconstrained ordination
using rdaand cca. The software Canoco 5 (ter Braak & Smilauer, 2018) contains Cfit and a rescaled
version of the column Var but not the columns Fratio.k.
A matrix with taxa in rows and taxon specific statistics in columns:
sample variance of the taxon.
R squared: fraction of the sample variance explained by the predictors and covariates.
fraction of the sample variance explained by the predictors after adjustment for the covariates.
fraction of the sample variance explained by the first k axes (the axes are constrained in a constrained ordination and unconstrained otherwise).
F-ratio of the regression model with k axes (as predictors) and the model without axes (i.e. with covariates or the intercept only).
ter Braak C.J.F. & Šmilauer P. (2018): Canoco reference manual and user's guide: software for ordination, version 5.1x. Microcomputer Power, Ithaca, USA, 536 pp. (http::www.canoco5.com)
data(pyrifos, package = "vegan") #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))) library(vegan) mod_rda <- rda( pyrifos ~ dose:week + Condition(week), data = Design) species_scores= scores(mod_rda,choices =1:2, display="sp",scaling = 1) library(PRC) head(rdacca_taxon_stats(mod_rda)) # stats per species get_focal_and_conditioning_factors(mod_rda) # factors from the model species_scores= vegan::scores(mod_rda,choices =1:2, display="sp",scaling = 2) # 1-d vertical plot of species scores (a cut of RDA1-scores at 0.5) plot_species_scores_bk(species_scores) # a cut of RDA1-scores at 1.0 (as scoresname is unset) plot_species_scores_bk(species_scores, threshold = 14) plot_species_scores_2d(species_scores) plot_species_scores_2d(species_scores, label.repel = TRUE) plot_species_scores_2d(species_scores, max.overlaps = 30, label.repel = TRUE)data(pyrifos, package = "vegan") #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))) library(vegan) mod_rda <- rda( pyrifos ~ dose:week + Condition(week), data = Design) species_scores= scores(mod_rda,choices =1:2, display="sp",scaling = 1) library(PRC) head(rdacca_taxon_stats(mod_rda)) # stats per species get_focal_and_conditioning_factors(mod_rda) # factors from the model species_scores= vegan::scores(mod_rda,choices =1:2, display="sp",scaling = 2) # 1-d vertical plot of species scores (a cut of RDA1-scores at 0.5) plot_species_scores_bk(species_scores) # a cut of RDA1-scores at 1.0 (as scoresname is unset) plot_species_scores_bk(species_scores, threshold = 14) plot_species_scores_2d(species_scores) plot_species_scores_2d(species_scores, label.repel = TRUE) plot_species_scores_2d(species_scores, max.overlaps = 30, label.repel = TRUE)
Define the model for smooth PRC using LMMsolver (Boer 2023) and smoothPRC
set_smoothPRC_model( spline = NULL, splineH0 = ~spl1D(time, nseg = 6, degree = 3, pord = 2), fixed = y ~ 1, fixedH0 = fixed, random = NULL, residual = NULL, weights = NULL, randomH0 = random, treatment.coding = TRUE, pord = c(time = 2, dose = 2, Treatment = 2), nseg = c(time = 6, dose = 6), degree = 3, start_time = 0, referencelevel = NULL, method = "rda", scaling = "ms", data )set_smoothPRC_model( spline = NULL, splineH0 = ~spl1D(time, nseg = 6, degree = 3, pord = 2), fixed = y ~ 1, fixedH0 = fixed, random = NULL, residual = NULL, weights = NULL, randomH0 = random, treatment.coding = TRUE, pord = c(time = 2, dose = 2, Treatment = 2), nseg = c(time = 6, dose = 6), degree = 3, start_time = 0, referencelevel = NULL, method = "rda", scaling = "ms", data )
spline |
A formula for the spline part of the model.
Should be of the form "~ spl1D()", ~ spl2D()" or "~spl3D()".
Generalized Additive Models (GAMs) can also be used,
for example "~ spl1D() + spl2D()".Default: spl1D(time, nseg = 6, degree = 3, pord = 2) + spl1D(D1, nseg = 6, degree = 3, pord = 2) + spl1D(D2, nseg = 6, degree = 3, pord = 2) + spl1D(D3, nseg = 6, degree = 3, pord = 2) with |
splineH0 |
As spline but for the null model. Default: a one dimensional
smooth in time.
|
fixed |
A formula for the fixed part of the model. Should be of the form "response ~ pred" |
fixedH0 |
Formula for the null model.
The default is equal to |
random |
Formula of the form "~ pred" for the random term
in the full model fit. The default is |
residual |
A formula for the residual part of the model. Should be of the form "~ pred". |
weights |
A character string identifying the column of data to use as relative weights in the fit. Default value NULL, weights are all equal to one. |
randomH0 |
The default is equal to |
treatment.coding |
logical. Default |
pord |
named integer vector of the difference order of the
spline penalty.
Default: |
nseg |
named integer vector of number of segments in P-splines,
can be large.
Default |
degree |
degree of the B-splines. Default |
start_time |
numeric time of the start of the treatment application.
Default: |
referencelevel |
name of the reference treatment. Default: first level of factor(data$Treatment). |
method |
character rda or cca. Default: |
scaling |
The default gives mean square species scores of 1, otherwise the sum of squares of the species scores is 1. See doPRC. Default: "ms". |
data |
A data.frame containing the modeling data. |
An atypical aspect of the function is that the time and treatment factors
need to have the names Time and Treatment in data, i.e.
set_smoothPRC_model assumes that the data contain two factor or
character variables
Time and Treatment that can be converted to the numeric
variables time and dose, respectively, using
fvalues4levels. The values of dose are set to 0 for
negative values of time and time ==0,
after subtraction of start_time.
In this case, the default spline model is:
~ spl2D(time, dose, seg = c(6, 6), degree = 3, pord = 2).
To enable Treatment to be treated as a factor (instead of as a
quantitative variable), the function generates variables with
names D1, ..., Dk, with k the number of
treatments beyond the reference treatment,
if treatment.coding is FALSE.
These variables generated in the function, contain the time
since the start of the experiment for each treatment.
D1, D2, ... should therefore not be variable names in data,
nor should time or dose be variable names
in data; these are reserved names in smoothPRC.
An example of the formula for spline
when treatment.coding is FALSE
and the number of levels of Treatment is 4, is
spl1D(time, nseg = 6, degree = 3, pord = 2) + spl1D(D1, nseg = 6, degree = 3, pord = 2) + spl1D(D2, nseg = 6, degree = 3, pord = 2) + spl1D(D3, nseg = 6, degree = 3, pord = 2)
This is the default model if treatment.coding is FALSE.
The default spline model assumes that the Treatment factor has a quantitative
basis, so that the response to Time and Treatment can be analyzed using a
2D spline model.
Boer, Martin P. 2023. Tensor Product P-Splines Using a Sparse Mixed Model Formulation. Statistical Modelling 23 (5-6): 465–79. doi:10.1177/1471082X231178591
data("Ossenkampen") names(Ossenkampen)[4:5] names(Ossenkampen)[4:5] <- c("Time", "Treatment") test <- FALSE if (test){ years <- sort(unique(Ossenkampen$Time)) years<- years[c(-(1:9))] years <- years[-which(years%in% c(1984))] years ids <- which(Ossenkampen$Time %in% years) }else ids <- 1:nrow(Ossenkampen) Y <- Ossenkampen[ids,-(1:5)] Design <- Ossenkampen[ids,1:5] Design$Block <- factor(Design$Block) # in doPRC and smoothPRC Time must be a factor, not a quantitative variable Design$Time <- factor(Design$Time) with(Design,table(Time,Treatment)) # unbalanced!!! if test <- FALSE mod_prc <- doPRC(Y ~ Time:Treatment + Condition(Time+Block), data = Design) b_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "sp") smoothPRC_model <- set_smoothPRC_model(data= Design, treatment.coding =FALSE, start_time = 1958) smoothPRC_model$fixed smoothPRC_model$spline smooth_PRC <- smoothPRC(Y, model = smoothPRC_model, formula_covariates = ~Block) smooth_PRC$iter if (1==0){ # start the iterations from previous opt<- smooth_PRC$options_iter opt$init <- smooth_PRC$X smooth_PRC2 <- smoothPRC(Y, model = smoothPRC_model, formula_covariates = ~Block, options_iter = opt) smooth_PRC2$iter round(diag(cor(smooth_PRC$X,smooth_PRC2$X )),6) round(diag(cor(smooth_PRC$B[,1:2],smooth_PRC2$B[,1:2])),6) } summary(smooth_PRC$obj[[1]]) cor(smooth_PRC$B[,1], b_mod_prc) plotsmoothPRC(smooth_PRC, flip = TRUE)data("Ossenkampen") names(Ossenkampen)[4:5] names(Ossenkampen)[4:5] <- c("Time", "Treatment") test <- FALSE if (test){ years <- sort(unique(Ossenkampen$Time)) years<- years[c(-(1:9))] years <- years[-which(years%in% c(1984))] years ids <- which(Ossenkampen$Time %in% years) }else ids <- 1:nrow(Ossenkampen) Y <- Ossenkampen[ids,-(1:5)] Design <- Ossenkampen[ids,1:5] Design$Block <- factor(Design$Block) # in doPRC and smoothPRC Time must be a factor, not a quantitative variable Design$Time <- factor(Design$Time) with(Design,table(Time,Treatment)) # unbalanced!!! if test <- FALSE mod_prc <- doPRC(Y ~ Time:Treatment + Condition(Time+Block), data = Design) b_mod_prc <- vegan::scores(mod_prc,choices= 1, display = "sp") smoothPRC_model <- set_smoothPRC_model(data= Design, treatment.coding =FALSE, start_time = 1958) smoothPRC_model$fixed smoothPRC_model$spline smooth_PRC <- smoothPRC(Y, model = smoothPRC_model, formula_covariates = ~Block) smooth_PRC$iter if (1==0){ # start the iterations from previous opt<- smooth_PRC$options_iter opt$init <- smooth_PRC$X smooth_PRC2 <- smoothPRC(Y, model = smoothPRC_model, formula_covariates = ~Block, options_iter = opt) smooth_PRC2$iter round(diag(cor(smooth_PRC$X,smooth_PRC2$X )),6) round(diag(cor(smooth_PRC$B[,1:2],smooth_PRC2$B[,1:2])),6) } summary(smooth_PRC$obj[[1]]) cor(smooth_PRC$B[,1], b_mod_prc) plotsmoothPRC(smooth_PRC, flip = TRUE)
Simulated data with a 10x5 factorial design with two replications. The true PRC curves are diverging straight lines as given in FigS2 of ter Braak (2023) with an interaction effect that is half the main effect of the second factor. The data and design as analysed in ter Braak (2023) has one empty cell (cell A10, B5).
The data has 100 unit and 10 colums as follows:
case identifier for a sample e.g. A2B3.1 the replicate of cell
for the levels 2 and 3 of the factors A and B, respectively.
A factor A, taken as Time variable in the classical PRC of ter Braak (2023)
B factor B, taken as Treatment variable in the classical PRC of ter Braak (2023)
V1 value of response variable 1 and so on for the remaining 6
response variables
ter Braak, C.J.F. (2023) Redundancy analysis includes analysis of variance-simultaneous component analysis (ASCA) and outperforms its extensions. Chemometrics and Intelligent Laboratory Systems, 104898. doi:10.1016/10.1016/j.chemolab.2023.104898
Smooth PRC using LMMsolver (Boer, 2023)
smoothPRC( response, model, options_iter = list(init = NULL, k = 2, mA = 4, tol = 1e-08, maxiter = 50), n_axes = 2, formula_covariates = NULL, weights = NULL, n_species = 1 )smoothPRC( response, model, options_iter = list(init = NULL, k = 2, mA = 4, tol = 1e-08, maxiter = 50), n_axes = 2, formula_covariates = NULL, weights = NULL, n_species = 1 )
response |
matrix or data frame of the abundance data or a distance
matrix or object of class |
model |
a named list.
Usually the result of
|
options_iter |
list of five elements
|
n_axes |
the number of axes to extract. Default |
formula_covariates |
one-sided formula with row covariates
(sample-specific predictors) that the model needs to adjust for.
Default |
weights |
a list of user-defined weights for columns (species) first
and, second, rows (sites) Default |
n_species |
number of species in |
An atypical aspect of the function is that the time and treatment factors
need to have the names Time and Treatment in data.
Classical PRC is a reduced rank model of rank 1 (or 2) with model
~Treatment*Time with the Time factor is a concomitant variable.
In smoothPRC the factor Time is taken quantitatively and
the main effect
of time is not an un-ordered set of values but a smooth P-spline model.
Also, the time-dependent effects of the treatments are smooth P-splines.
As the main and interaction effects are splines, the main effect of time
cannot be adjusted for in the usual way. It is taken out via a separate
spline fit using the formulas in splineH0 and fixedH0 of
set_smoothPRC_model.
The algorithm solves a penalized reduced rank model for the
dominant eigen vector (i.e. rank 1).
The model is a P-spline model and
the penalty parameters are obtained by REML using
LMMsolve (Boer 2023).
The model is linear in, possibly tensor, B-splines,
but through the B-splines non-linear in the original quantitative variables,
here time or dose
and thus non-linear from the perspective of the user.
With quadratic ridge-type difference penalties on the coefficients,
the B-splines become P-splines (penalized B-splines).
For an initial vector b of species scores,
the unconstrained sites
are calculated, which are then used as response vector in
LMMsolve to give a constrained site scores
x by the P-spline model specified in model. This
is part of a power algorithm with block and Anderson acceleration to solve
the dominant eigen vector of the
underlying eigen problem while optimizing simultaneously the penalties
using LMMsolve.
The smoothPRC algorithm uses one and two-sided power iterations.
It thus iterates through the transition formulas of the method with
(P-spline) regression fits.
It uses two acceleration methods:
blocking and Anderson iteration.
The response is taken as a distance object if
it inherits class dist or if it a matrix-like object with
equal dimensions using dim and all zero diagonal using
sum(diag(response))==0.
Initialization of the power algorithm expands the
dist object in response
to a square symmetrix matrix.
This expansion can be avoided by explicit initialization
using options_iter$init.
If response is a distance object and the maximum distance exceeds 4
then the distance is divided by n-1 with n the number of samples.
This choice is taken from vegan::dbrda and usually makes the eigenvalues
equal to those of vegan::rda if the distance is 'euclidean'.
A list.
eigenvalues.
matrix with response variable (species) scores B's
with Fratio's with attached axis number, e.g. B1.
matrix with the eigen vectors (i.e. the constrained sample scores) in columns.
matrix with the unconstrained sample scores.
matrix with the PRC scores, X scores compared to the X-scores of
the reference level in Treatment in data
matrix with the unconstrained PRC scores.
names list of columns weights and row weights. Names are
K, R, sqrtK, sqrtR. The first two are
column (species) and row (site)totals in with method != "rda".
They not normalized; none sums to 1. All 1 in method =="rda"
.
a named list of LMMsolve objects,
one for each axis.
a named list of LMMsolve objects
of the null model, one for each axis.
the number of iterations of the accelerated power algorithm.
the response matrix.
copy of the model argument to the
powerSVD.
copy of the opions_iter argument
to the smoothPRC.
an internal matrix: Intercept, optionally covariates and previous axes.
the number of user-defined row covariates using
formula_covariates with one for the "(Intercept)".
Boer, Martin P. 2023. Tensor Product P-Splines Using a Sparse Mixed Model Formulation. Statistical Modelling 23 (5-6): 465–79. doi:10.1177/1471082X231178591
Updates the parameters of solver (internal)
update_fit_state(object, state)update_fit_state(object, state)
object |
Object resulting for the fit. |
state |
Optional list of arguments for warm starts of the solver |