For this model structure, just like for the univariate version seen
in the Starter Pack, we assume that species do not interact with one
another. This implies that species interactions don’t affect the
ecosystem function responses leaving only the species proportions as
fixed effects, along with any experimental structures. An example of
this model is as follows:
\[\large y_{mt}=\sum_{i=1}^{S=4} \beta_{it}
p_{im} + \epsilon_{mt}\]
Where \(y_{mt}\) represents the
recorded value, \(y\), of the ecosystem
function at the time point \(t\) (where
\(t\) goes from \(1\) to \(T=3\)) from the \(m^{th}\) experimental/observational study
unit (where \(m\) goes from \(1\) to \(M=384\)). The coefficient \(\beta_{it}\) scaled by the initial
proportion \(p\) of species \(i\) forms the ID effect of said species for
the ecosystem function response at time \(t\). The error of this model, \(\epsilon\) is assumed to follow a
multivariate Normal distribution with mean 0 and variance \(\Sigma^{*}\), where \(\Sigma^{*}\) is an \(MT \times MT\) block-diagonal matrix.
\[\Sigma^{*} =
\begin{bmatrix}
\Sigma_1 & & 0\\
& \ddots & \\
0 & & \Sigma_M
\end{bmatrix}\]
Each matrix along the main diagonal, \(\Sigma_{m}\) represents the
variance-covariance matrix of the time points. There are a few different
forms that this may take (these are called auto-correlation structures).
The first of which is the unstructured or general structure (UN), where
each element is allowed to differ from one another.
\[\Sigma_m =
\begin{bmatrix}
\sigma_1^2 & \sigma_{1,2} & ... & \sigma_{1,T}\\
\sigma_{1,2} & \sigma_{2}^2 & ... & \sigma_{2,T}\\
\vdots & \vdots & \ddots & \vdots \\
\sigma_{1,T} & \sigma_{2, T} & ... & \sigma_T^2
\end{bmatrix}\]
IDModel_UN <- DImulti(y = "Y", time = c("time", "UN"), unit_IDs = "plot",
prop = 2:5, data = simRM, DImodel = "ID", method = "REML")
Another of which is the compound symmetry (CS) structure, which
assumes that all variances are equal and all covariances are equal in
the matrix. So for our three time points, it would take the form
below.
\[\Sigma_m = \sigma^2
\begin{bmatrix}
1 & \rho & \rho\\
\rho & 1 & \rho\\
\rho & \rho & 1
\end{bmatrix} \]
IDModel_CS <- DImulti(y = "Y", time = c("time", "CS"), unit_IDs = "plot",
prop = 2:5, data = simRM, DImodel = "ID", method = "REML")
The final structure available for use in the R package is an
autoregressive model of order one (AR(1)), which assumes that all
variances are equal in the matrix and all covariances vary by a degree
() for ech measure of distance between readings, i.e. time point 1 and 2
will share a covariance with points 2 and 3, but time points 1 and 3
will have the same covariance times \(\rho\). Therefore, for our three time
points, it would take the form below.
\[\Sigma_m = \sigma^2
\begin{bmatrix}
1 & \rho & \rho^2\\
\rho & 1 & \rho\\
\rho^2 & \rho & 1
\end{bmatrix} \]
IDModel_AR1 <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = 2:5, data = simRM, DImodel = "ID", method = "REML")
As no fixed effects vary between these models, we can fit them using
the REML estimation method, however, when comparing this model with
another which has different fixed effects, we would use ML, as
below.
IDModel <- DImulti(y = "Y", time = c("time", "UN"), unit_IDs = "plot",
prop = 2:5, data = simRM, DImodel = "ID", method = "ML")
And just like with the univariate models, we can use summary() to
look at our model output.
summary(IDModel)
## Generalized least squares fit by maximum likelihood
## Model: Y ~ 0 + time:((p1_ID + p2_ID + p3_ID + p4_ID))
## Data: data
## AIC BIC logLik
## 1296.106 1365.445 -630.0529
##
## Correlation Structure: General
## Formula: ~0 | plot
## Parameter estimate(s):
## Correlation:
## 1 2
## 2 0.143
## 3 0.273 0.376
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~0 | time
## Parameter estimates:
## 1 2 3
## 1.0000000 1.2525790 0.8717619
##
## Coefficients:
## Value Std.Error t-value p-value
## time1:p1_ID 6.150533 0.4589786 13.400479 0
## time2:p1_ID 5.320499 0.5749069 9.254540 0
## time3:p1_ID 3.442754 0.4001200 8.604302 0
## time1:p2_ID 4.793913 0.4800486 9.986307 0
## time2:p2_ID 6.381084 0.6012988 10.612168 0
## time3:p2_ID 5.136560 0.4184881 12.274089 0
## time1:p3_ID 4.791368 0.4102820 11.678232 0
## time2:p3_ID 7.441564 0.5139106 14.480271 0
## time3:p3_ID 7.227109 0.3576682 20.206183 0
## time1:p4_ID 7.098678 0.4374651 16.226842 0
## time2:p4_ID 10.484815 0.5479596 19.134283 0
## time3:p4_ID 9.099688 0.3813654 23.860808 0
##
## Theta values: 1
##
##
## Correlation:
## t1:1_I t2:1_I t3:1_I t1:2_I t2:2_I t3:2_I t1:3_I t2:3_I t3:3_I
## time2:p1_ID 0.143
## time3:p1_ID 0.273 0.376
## time1:p2_ID -0.133 -0.019 -0.036
## time2:p2_ID -0.019 -0.133 -0.050 0.143
## time3:p2_ID -0.036 -0.050 -0.133 0.273 0.376
## time1:p3_ID -0.188 -0.027 -0.051 -0.194 -0.028 -0.053
## time2:p3_ID -0.027 -0.188 -0.071 -0.028 -0.194 -0.073 0.143
## time3:p3_ID -0.051 -0.071 -0.188 -0.053 -0.073 -0.194 0.273 0.376
## time1:p4_ID -0.301 -0.043 -0.082 -0.283 -0.040 -0.077 -0.090 -0.013 -0.024
## time2:p4_ID -0.043 -0.301 -0.113 -0.040 -0.283 -0.106 -0.013 -0.090 -0.034
## time3:p4_ID -0.082 -0.113 -0.301 -0.077 -0.106 -0.283 -0.024 -0.034 -0.090
## t1:4_I t2:4_I
## time2:p1_ID
## time3:p1_ID
## time1:p2_ID
## time2:p2_ID
## time3:p2_ID
## time1:p3_ID
## time2:p3_ID
## time3:p3_ID
## time1:p4_ID
## time2:p4_ID 0.143
## time3:p4_ID 0.273 0.376
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.20807103 -0.63836145 0.04442358 0.65743274 2.35308203
##
## Residual standard error: 1.493056
## Degrees of freedom: 348 total; 336 residual
