From multivariate regression to partial redundancy analysis

Define the model:

Y = QA\(_{0}\) + Z\(_{0}\)M + E

with Y the response matrix, Q a matrix with concomitant variables, Z\(_{0}\) a matrix of predictor, E an error matrix, and A\(_{0}\) and M matrices of unknown regression coefficients that must be estimated from the data Y, Q and Z\({_0}\).

In partial RDA, the focus lies on a low-dimensional representation of the effects of Z\({_0}\) on the response, that is on M, with the effects of Q partial out. To this aim and with \(\Pi_{Q}\) the orthogonal projector on Q, define Z = Z\(_{0}\) - \(\Pi_{Q}\)Z\(_{0}\), the matrix of residuals of the regression of Z\(_{0}\) on to Q. Now, Q and Z are orthogonal and span the same space as Q and Z\(_{0}\). Moreover, the regression coefficients associated with Z are the same as those with Z\(_{0}\), while the regression coefficients associated with Q change. The original model can thus be rewritten to

Y = QB + ZM + E

The low-dimensional representation is obtained by a rank-restriction on M. With the rank-1 restriction M = cb\(^T\), where c is a vector of canonical coefficients and b is a vector of loadings (species scores),

ZM = Zcb\(^T\) = xb\(^T\), which is a partial PCA model with constraint x = Zc with x the constrained site scores. Partial redundancy analysis is thus both reduced rank regression with concomitant variables and a partial PCA with constraints. Note that the matrix of regression coefficient B in the rewritten formula can be obtained by regression Y on Q only. Note also that the F- or LR-test on nullity of M is also a comparison of the full model with both Q and Z0 and the null model with Q only.

The transition formulas of classical partial redundancy analysis

The transition formulas specify the eigen problem of partial RDA and are useful to solve for the first eigen vector by iteration. The convergence of the iteration is accelerated by blocking and Anderson iteration. The transition formulas of partial redundancy analysis are:

x\(^*\) = Yb

c = (Z\(^T\)Z)\(^{-1}\)Z\(^T\)x\(^*\)

x = Zc

\(\lambda\)b=Y\(^T\)x

This set of equation forms the eigen problem

\(\lambda\)b=Y\(^T\)Z(Z\(^T\)Z)\(^{-1}\)Z\(^T\)Yb

This eigen equation is obtained by substitution of symbols, starting with the last transition formula and ending with substitution of x\(^*\) to Yb of the first equation.

In classical PRC, Q is an indicator matrix defined by the factor time and Z\(_{0}\) is an indicator matrix defined by the interaction of the factors time and treatment. In words, the main effects of the factor time are partialled out of the interaction effects and the reduced rank (dimension reduction) is applied to the time-dependent treatment effects. In classical PRC, x = Zc with Z orthogonal to Q and c obtained from regression of x\(^*\) on to Z. In smoothPRC below, Z is not explicitly easily available, whereas Z\(_{0}\) is. But, the regressions of x\(^*\) on Q and Z\(_{0}\) and of x\(^*\) on Q and Z yield identical fitted values, say y\(_{1}\).

Consequently,

x = Zc = y\(_{1}\) - Qa

where y\(_{1}\) is obtained from a regression on Q and Z\(_{0}\) and Qa from a regresion of x\(^*\) on Q alone. In pseudo R-code:

x = fitted(lm(x\(^*\) ~ Q + Z\(_{0}\))) - fitted(lm(x\(^*\) ~ Q)).

By consequence, Zc can be computed as a difference of the fitted values of two models, the full model with Q and Z\(_{0}\) and a null model with Q alone; in PRC, the model with time and treatment and the model with time only, respectively.

The transition formulas of smooth PRC

In smoothPRC time is treated as a quantitative variable. The model is allowed to be highly non-linear by transformation to a B-spline basis. The single variable time is thereby replaced by a set of variables (B-splines), like time is replaced by a set of indicator variables in classical PRC. As there can and should be many variables in the B-spline basis, the model is over-parametrized and the partial RDA becomes an unconstrained partial PCA. That is why penalized partial RDA is used in smoothPRC. With a treatment with a quantitative basis, a dose, smoothPRC uses a two-dimensional tensor P-spline (Boer 2023), that is, a Kronecker product of B-splines with two penalty parameters, one for the time and the other for treatment. In smoothPRC, the main effect of time is a one-dimensional smooth in time. Both models are fitted using the R package LMMsolver (Boer 2023). The main effects of time cannot be partialled out as in classical PRC as the smooths are not projectors. Therefore, the formula in classical PRC, in pseudo R-code

x = Zc = fitted(x\(^*\) ~ Q + Z\(_{0}\))) - fitted(lm(x\(^*\) ~ Q))

is replaced in smoothPRC by the difference of the fitted values of a model with both time and treatment and the fitted values of a null model with only time. In the default smoothPRC model, it is the difference of the fitted values of a 2d-spline of time and treatment and the fitted values of a 1d-spline of time.

In summary, the transition formulas of smoothPRC are the three equations, of which the second is motivated above,

x\(^*\) = Yb

x = fitted(smooth of time and treatment) - fitted(smooth of time)

\(\lambda\)b=Y\(^T\)x

After convergence

After convergence, x\(^*\) is replaced by x\(^*\) - fitted(smooth of time), so a to remove the time effect from the ‘unconstrained’ site scores. This is the \(x\_star\) entry in the output of smoothPRC.

Also, the PRC graph takes one level or value of the treatment as reference level or value. After convergence, the values of the PRC curves are obtained from the full model by subtracting from the fitted values the predictions for the data set derived from the original data with all treatment levels or values replaced by the reference level or value (0). This is the \(PRC\) entry in the output of smoothPRC.

By default, all scores are rescaled such that the mean square of the species scores \(b\) is 1. This makes smoothPRC behave as doPRC with scaling = "ms".

Example of the pyrifos data

The example uses the data from the chlorpyrifos experiment (Van den Brink and Ter Braak 1999) available in the vegan package, with the design of the expriment: 12 cosms with 5 treatments inclusing the control (no chlorpyrifos) measured at 11 time points. The names Time and Treatment are expected in smoothPRC.

#the chlorpyrifos experiment from van den Brink & ter Braak 1999
data(pyrifos, package = "vegan") #log-transformed species data from package vegan
Y <- pyrifos
Design <- data.frame(Time=gl(11, 12, labels=c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24)),
                     Treatment=factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11)),
                     cosm = gl(12, 1, length=132))

The default model is set up by:

library(PRC)
smoothPRC_model  <- set_smoothPRC_model(data = Design)

The model specifications for the null model and alternative models are:

smoothPRC_model$splineH0
#> ~spl1D(time, nseg = 6, degree = 3, pord = 2)
#> <environment: 0x6450a826bb90>
smoothPRC_model$spline
#> ~spl2D(time, dose, nseg = c(6, 6), degree = 3, pord = 2)
#> <environment: 0x6450a826bb90>

The null model is a spline function of time only, whereas the alternative model uses a two-dimensional spline model in time and dose. The time variable is the quantitative version of the variable Time in the Design.The dose variable is the quantitative version of the variable Treatment with dose set to 0 for pre-application times (negative times by default and, in general, times before start_time. The warnings are from LMMsolve and are expected in early iterations of the power algorithm when the number of LMMsolve iterations is then limited. In smoothPRC, the maximum number of LMMsolve iterations is set to 50. If the number of power iterations exceeds 50 a warning is issued such as No convergence of the power algorithm for axis 2. Such non-convergence is expected when the eigenvalue is close in value to the next one, in the example axis 3.

smooth_PRC <- smoothPRC(Y, model = smoothPRC_model)
#> Warning: No convergence after 10 iterations
#> Warning: No convergence after 10 iterations
#> Warning: No convergence after 10 iterations
#> Warning: No convergence after 10 iterations
#> Warning: No convergence after 10 iterations
#>  No convergence of power algorithm for axis  2
#> Warning in smoothPRC(Y, model = smoothPRC_model):
summary(smooth_PRC$obj[[1]])
#> Table with effective dimensions and penalties: 
#> 
#>             Term Effective Model Nominal Ratio Penalty
#>      (Intercept)      1.00     1       1  1.00    0.00
#>  lin(time, dose)      3.00     3       3  1.00    0.00
#>          s(time)     11.38    81      77  0.15    0.09
#>          s(dose)      4.43    81      77  0.06    0.00
#>         residual    112.19   132     128  0.88    0.30
#> 
#>  Total Effective Dimension: 132
summary(smooth_PRC$objH0[[1]])
#> Table with effective dimensions and penalties: 
#> 
#>         Term Effective Model Nominal Ratio Penalty
#>  (Intercept)      1.00     1       1  1.00    0.00
#>    lin(time)      1.00     1       1  1.00    0.00
#>      s(time)      0.86     9       7  0.12   44.33
#>     residual    129.14   132     130  0.99    0.03
#> 
#>  Total Effective Dimension: 132

Illustration that smoothPRC extends classical PRC

classical_PRC_model <- set_smoothPRC_model(
  fixed= y ~ Treatment:Time,
  fixedH0 = y ~Time,
  spline = NULL,
  splineH0 = NULL,
  data = Design)
classical_PRC <- smoothPRC(Y, model = classical_PRC_model)
summary(classical_PRC$obj[[1]])
#> Table with effective dimensions and penalties: 
#> 
#>            Term Effective Model Nominal Ratio Penalty
#>     (Intercept)         1     1       1     1    0.00
#>  Treatment:Time        54    54      54     1    0.00
#>        residual        77   132      77     1    0.32
#> 
#>  Total Effective Dimension: 132
summary(classical_PRC$objH0[[1]])
#> Table with effective dimensions and penalties: 
#> 
#>         Term Effective Model Nominal Ratio Penalty
#>  (Intercept)         1     1       1     1    0.00
#>         Time        10    10      10     1    0.00
#>     residual       121   132     121     1    0.03
#> 
#>  Total Effective Dimension: 132

Classical PRC is obtained by:

mod_prc <- doPRC(pyrifos ~ Treatment:Time + Condition(Time),  data = Design)
#> 
#> Some constraints or conditions were aliased because they were redundant. This
#> can happen if terms are constant or linearly dependent (collinear):
#> 'Treatment44:Time-4', 'Treatment44:Time-1', 'Treatment44:Time0.1',
#> 'Treatment44:Time1', 'Treatment44:Time2', 'Treatment44:Time4',
#> 'Treatment44:Time8', 'Treatment44:Time12', 'Treatment44:Time15',
#> 'Treatment44:Time19', 'Treatment44:Time24'

The absolute correlation between the last two PRC models is 1, 0.068738, whereas the correlation between the smooth and classical PRC models is 0.97, 0.09, 0.05, 0.63

Here is the PRC diagram showing both models with species scores of the smooth version.

plotsmoothPRC(smooth_PRC, flip = FALSE, threshold = 1)
#> 122  species with names out of 178 species

Here is the plot for classical PRC via the plotting function of smoothPRC:

plotsmoothPRC(classical_PRC, flip = TRUE, title ="Classical PRC via smoothPRC")
#> 50  species with names out of 178 species

classical_PRC$eig/mod_prc$CCA$eig[seq_along(classical_PRC$eig)]
#>      eig1 eig2
#> [1,]    1    1

The first two eigenvalues agree up to numerical precision and the solid lines for smooth PRC overlay the dotted lines of classical PRC, illustrating that smooth PRC can reproduce the outcome of a classical PRC and thus is a real extension of classical PRC.

With Euclidean distance, distance-based smooth PRC can yield the same result:

D <- vegan::vegdist(Y, method = "euclidean") 
out_D<- smoothPRC(D, model = classical_PRC_model)
out_D$eig/mod_prc$CCA$eig[seq_along(classical_PRC$eig)]
#>      eig1 eig2
#> [1,]    1    1

Distance-based smooth PRC generalizes smooth PRC on Y by allowing other choices of distance.