Modelos Jerárquicos de Comunidades de Especies en R

Hierarchical Models of Species Communities (HMSC)

Diego Nieto Lugilde

Universidad de Córdoba (España)

Cargamos algunos datos

library(openxlsx)
env <- read.xlsx("data/Wilson_Latorre/Data_multi_env_R_course.xlsx",
                 sheet = 1, 
                 rowNames = TRUE)
com <- read.xlsx("data/Wilson_Latorre/Data_multi_bio_R_course.xlsx",
                 sheet = 1, 
                 rowNames = TRUE)
sites <- read.xlsx("data/Wilson_Latorre/Data_multi_site_R_course.xlsx",
                   sheet = 1,
                   rowNames = TRUE)

Preparamos un poco los datos

comm <- as.matrix(com)

envm <- scale(env)
envm <- as.data.frame(envm)

Definimos el diseño experimental

library(Hmsc)
studyDesign <- data.frame(sitio = as.factor(sites$Sitio),
                          parcela = as.factor(sites$Parcela),
                          año = as.factor(sites$Año),
                         estacion = as.factor(sites$Estacion))
sitioL <- HmscRandomLevel(units = levels(sites$Sitio))
parcelaL <- HmscRandomLevel(units = levels(sites$Parcela))
añoL <- HmscRandomLevel(units = levels(sites$Año))
estacionL <- HmscRandomLevel(units = levels(sites$Estacion))


sitioL$nfMin <- 2
parcelaL$nfMin <- 2
añoL$nfMin <- 2
estacionL$nfMin <- 2

Especificamos la estructura del modelo

m0 <- Hmsc(Y = comm, 
         XData = envm, 
         XFormula = ~1,
         studyDesign = studyDesign, 
         ranLevels = list("sitio" = sitioL, 
                          "parcela" = parcelaL, 
                          "año" = añoL, 
                          "estacion" = estacionL),
         distr = "lognormal poisson")

m1 <- Hmsc(Y = comm, 
           XData = envm,
           XFormula = ~`N-NO3-` + `N-NH3` + `N-NO2-` +
             `PO43-` + pH + DO + EC + NTU + `T°C` +
             Current.velocity + Average.depth + Average.width,
           distr = "lognormal poisson")


m_env <- Hmsc(Y = comm, 
          XData = envm,
           XFormula = ~`N-NO3-` + `N-NH3` + `N-NO2-` +
            `PO43-` + pH + DO + EC + NTU + `T°C` +
             Current.velocity + Average.depth + Average.width,
          studyDesign = studyDesign, 
         ranLevels = list("sitio" = sitioL, 
                          "parcela" = parcelaL, 
                          "año" = añoL, 
                          "estacion" = estacionL),
           distr = "lognormal poisson")

Definimos los parámetros para las Cadenas de Markov Monte Carlo

samples <- 100 # Este valor debería ser más alto. Del orden de 10000 hacia arriba.
thin <- 1 # Este valor también suele ser más alto. Del orden de 50-100 hacia arriba.
transient <- 500*thin
nChains <- 4

Se calibran y estiman los parámetros del modelos con las MCMC

m0 <- sampleMcmc(m0, 
               thin = thin, 
               samples = samples, 
               transient = transient, 
               nChains = nChains, 
               nParallel = nChains)
m1 <- sampleMcmc(m1, 
               thin = thin, 
               samples = samples, 
               transient = transient, 
               nChains = nChains, 
               nParallel = nChains)
m_env <- sampleMcmc(m_env, 
               thin = thin, 
               samples = samples, 
               transient = transient, 
               nChains = nChains, 
               nParallel = nChains)

Convertimos los resultados en objetos CODA, son un tipo de datos de análisis bayesianos en R

mpost_m0 <- convertToCodaObject(m0)
mpost_m1 <- convertToCodaObject(m1)
mpost_m_env <- convertToCodaObject(m_env)

Calculamos dos estadísticos de calibración de MCMC

es_m0 <- effectiveSize(mpost_m0$Beta)
es_m1 <- effectiveSize(mpost_m1$Beta)
es_m_env <- effectiveSize(mpost_m_env$Beta)
gd_m0 <- gelman.diag(mpost_m0$Beta, multivariate=FALSE)$psrf
gd_m1 <- gelman.diag(mpost_m1$Beta, multivariate=FALSE)$psrf
gd_m_env <- gelman.diag(mpost_m_env$Beta, multivariate=FALSE)$psrf

Dibujamos las gráficas de esos dos parámetros

Effective size debería estar entorno al número de muestras de las MCMC (samples * nChains). El diagnóstico de Gelman debería estar entorno a 1.

par(mfrow=c(3,2))

hist(es_m0, main="ess(beta)")
hist(gd_m0, main="psrf(beta)")

hist(es_m1, main="ess(beta)")
hist(gd_m1, main="psrf(beta)")

hist(es_m_env, main="ess(beta)")
hist(gd_m_env, main="psrf(beta)")

… lo mismo pero con los valores de Omega, que son los coeficientes de interacciones entre especies

par(mfrow=c(2, 2))

hist(effectiveSize(mpost_m0$Omega[[1]]),
     main="ess(omega)")
hist(gelman.diag(mpost_m0$Omega[[1]], 
                 multivariate=FALSE)$psrf, 
     main="psrf(omega)")

hist(effectiveSize(mpost_m_env$Omega[[1]]),
     main="ess(omega)")
hist(gelman.diag(mpost_m_env$Omega[[1]],
                 multivariate=FALSE)$psrf, 
     main="psrf(omega)")

Tambien podemos dibujar algunas de las cadenas para visulizar su forma y ver si se han estabilizado

plot(mpost_m0$Beta)

Para evaluar el modelo en términos de la distribución/abundancia de las especies, debemos hacer predicciones primero.

preds_m0 <- computePredictedValues(m0, expected = FALSE)
evaluateModelFit(hM = m0, predY = preds_m0)
$RMSE
 [1]  2.1794495  0.9185587  1.0458250  0.3952847  1.0155048  1.9764235
 [7]  7.4885329  4.1250000  1.6863422  0.3952847  6.8248168 10.7172466
[13] 15.0052074  0.6731456  0.5994789  0.3535534  0.8291562 11.1139341
[19]  0.8592329  1.6105123 16.9822978  6.9664374 12.6336604  0.8003905
[25]  1.5155445  1.5103807  0.9457966  2.1176343  5.2055409  1.7139137
[31]  1.9685020  2.0991070  2.3091192 25.0341173  0.3644345  1.4523688
[37] 14.0562486  2.7556760  1.6770510 16.1086400  0.3535534  1.0307764
[43]  8.6250000  1.6298006  1.1792476  0.5728220 44.3586237

$SR2
 [1]  0.0301857013            NA  0.1118592883            NA            NA
 [6]            NA            NA            NA  0.6235870792 -0.0016398273
[11]  0.2916430541  0.3267520979            NA            NA            NA
[16]            NA            NA  0.3146942826  0.0254949354            NA
[21]  0.3228305490  0.3882955307  0.0191680859            NA  0.0005199556
[26]  0.0451046746  0.0027330455            NA  0.2212734720  0.1553161641
[31]  0.0525808215  0.2052764856  0.0265110112  0.5357610485  0.0074156470
[36]            NA  0.0559829460  0.1576894121  0.1354294172  0.0044618159
[41]            NA  0.0965918639            NA            NA            NA
[46]            NA  0.2151137191

$O.AUC
 [1] 0.5564784 0.4972678 0.7131410 0.4649123 0.6338798 0.9017857 0.8182957
 [8] 0.8262626 0.9539970 0.4353448 0.7617647 0.8376623 0.9000000 0.5273224
[15] 0.5343434 0.6406780 0.5448718 0.8815988 0.6728814 0.8750000 0.8454545
[22] 0.7924805 0.7312312 0.5471698 0.6655592 0.6232955 0.6322034 0.8474576
[29] 0.7591991 0.8555556 0.6394486 0.7778473 0.7109635 0.9140625 0.5431034
[36] 0.7562500 0.6787440 1.0000000 0.9335840 0.6892713 0.7474576 0.5785256
[43] 0.9365079 1.0000000 0.9435484 0.4437500 0.8160200

$O.TjurR2
 [1]  0.0094407530 -0.0007923497  0.0472275641 -0.0030639098  0.0086065574
 [6]  0.0679910714  0.0728571429  0.0931464646  0.2625603318 -0.0048419540
[11]  0.1827107843  0.2306764069  0.0565084746  0.0027322404  0.0059797980
[16]  0.0175423729  0.0064102564  0.2012330317  0.0118898305  0.0538750000
[21]  0.1181931818  0.1831250000  0.0467267267  0.0005617496  0.0206118494
[26]  0.0156818182  0.0140677966  0.0546694915  0.1046374459  0.1815648148
[31]  0.0247879109  0.0988548185  0.0414451827  0.1466964286  0.0048994253
[36]  0.0502083333  0.0323973430  0.2691250000  0.2697744361  0.0407945344
[41]  0.0171864407  0.0154647436  0.0409126984  0.0606854839  0.0793548387
[46] -0.0017500000  0.1069524406

$O.RMSE
 [1] 0.4873509 0.3623502 0.4591959 0.3651411 0.3623604 0.3907788 0.3585382
 [8] 0.3705943 0.3847182 0.4589676 0.4421826 0.4062763 0.3313888 0.3715266
[15] 0.4126335 0.3580388 0.4429125 0.3395069 0.5181089 0.3053666 0.4376108
[22] 0.4449304 0.4739357 0.4080402 0.5075016 0.5125173 0.5590210 0.3249818
[29] 0.4635378 0.4177758 0.5045804 0.4432973 0.5048105 0.2989291 0.4687441
[36] 0.3545790 0.4548103 0.4307380 0.3621189 0.4737365 0.3374812 0.4549806
[43] 0.2187560 0.3390662 0.2749975 0.3854659 0.4332120

$C.SR2
 [1]  8.156160e-03            NA  2.273929e-01            NA            NA
 [6]            NA            NA            NA -1.296155e-02            NA
[11]  3.159638e-01  4.660477e-01            NA            NA            NA
[16]            NA            NA  8.019433e-02  7.894737e-01            NA
[21]  8.264464e-03  4.042758e-01 -5.030541e-03            NA -3.155814e-03
[26]  6.465963e-03            NA            NA  2.904638e-01 -3.204494e-03
[31] -7.336096e-05  2.715683e-02  1.187001e-01  4.265523e-01            NA
[36]            NA  4.513932e-02 -4.000000e-02  4.275148e-01 -9.205181e-04
[41]            NA  2.217295e-01            NA            NA            NA
[46]            NA  1.006601e-01

$C.RMSE
 [1]  3.4610243  4.2426407  1.8743802  1.1952286  4.6904158  5.5901699
 [7] 22.6431951 11.0000000  3.3703118  1.2247449  9.1628990 17.8923355
[13] 53.6842621  3.1091264  1.5986105  1.2649111  1.9148542 12.1237279
[19]  0.6594585  6.4420494 29.9635844  9.4978650 19.4354088  1.9306146
[25]  2.0321763  1.9015680  0.7244884  7.5762788  8.4970240  3.9427107
[31]  2.7598916  3.4605381  3.3934652 26.6039467  0.9920365  5.8094750
[37] 26.4877150 10.4211514  4.6824083 20.5508599  1.2649111  2.2642910
[43] 69.0000000  9.2195445  6.6708320  2.2912878 85.7393132
preds_m1 <- computePredictedValues(m1, expected = FALSE)
evaluateModelFit(hM = m1, predY = preds_m1)
$RMSE
 [1]  2.2638463  0.9270248  1.2183493  0.3952847  1.0155048  1.9764235
 [7]  7.2929761  4.1268935  1.2437343  0.3750000  6.4529547 11.0474827
[13] 15.0052074  0.5994789  0.5994789  0.3535534  0.5728220 11.7659466
[19]  0.4145781  1.6105123 16.9204250  5.9621200 10.8836173  0.8003905
[25]  1.6441183  1.5309311  0.2795085  2.1176343  5.1158455  1.6535946
[31]  2.1794495  2.4173850  2.3881216 31.0164397  0.3061862  1.0383280
[37] 13.9306317  2.9348392  1.8708287  9.7207092  0.2795085  1.0915156
[43]  8.6250000  1.6298006  0.8385255  0.5728220 43.8897767

$SR2
 [1]  0.1178687785 -0.0007801073  0.0130607455            NA            NA
 [6]            NA  0.4567862715 -0.0025780348  0.4532401896            NA
[11]  0.3685336275  0.3873977967            NA  0.1585379280            NA
[16]            NA  0.2570729498  0.2640855867            NA            NA
[21]  0.2818839329  0.3570321517  0.2010412974            NA -0.0074093676
[26]  0.2971395727            NA            NA  0.5432342664  0.3360096347
[31]  0.0405365673 -0.0055690144  0.0853207221  0.3076218991            NA
[36]  0.7615085413  0.2443503626  0.0978884223 -0.0039454481  0.2852344437
[41]  0.2133333333  0.0875089517            NA            NA  0.9997480474
[46]            NA  0.2597392613

$O.AUC
 [1] 0.8028793 0.9262295 0.6241987 0.7167920 0.8224044 0.8772321 0.8558897
 [8] 0.8262626 0.9396682 0.6839080 0.7740196 0.8484848 0.7559322 0.6803279
[15] 0.7383838 0.4491525 0.8100962 0.8499246 0.8525424 0.7375000 0.8409091
[22] 0.8461914 0.7917918 0.7795883 0.6788483 0.7517045 0.4322034 0.5406780
[29] 0.8636364 0.7694444 0.6951220 0.6677096 0.7635659 0.8705357 0.8433908
[36] 0.8625000 0.8393720 0.9458333 0.8609023 0.7925101 0.8593220 0.7379808
[43] 0.9523810 0.9838710 1.0000000 0.9354167 0.9211514

$O.TjurR2
 [1]  0.116591916  0.191407104  0.053958333  0.051478697  0.068538251
 [6]  0.145133929  0.227788221  0.126444444  0.351142534  0.055977011
[11]  0.203504902  0.218885281  0.071491525  0.184480874  0.066934343
[16] -0.004838983  0.216282051  0.184736048  0.143576271  0.090291667
[21]  0.230318182  0.233515625  0.190072573  0.097002573  0.049202658
[26]  0.128738636 -0.033211864  0.016762712  0.270941558  0.143629630
[31]  0.078629374  0.055228411  0.161691584  0.137857143  0.131566092
[36]  0.352833333  0.202940821  0.190833333  0.166434837  0.205855263
[41]  0.212338983  0.157131410  0.110952381  0.194233871  0.539193548
[46]  0.153916667  0.212112015

$O.RMSE
 [1] 0.4275010 0.2319346 0.3939589 0.3056240 0.2261211 0.2990729 0.2829292
 [8] 0.3355488 0.2912931 0.2952528 0.4399098 0.3944131 0.2774447 0.2475974
[15] 0.3400918 0.2783891 0.3369767 0.3482809 0.2579922 0.2516979 0.3901677
[22] 0.4174201 0.4322659 0.3507257 0.4582077 0.4213516 0.3094332 0.2925541
[29] 0.3774607 0.3404257 0.4594969 0.4325231 0.4252330 0.2995680 0.2715636
[36] 0.2369100 0.3833776 0.2634592 0.2982129 0.4241788 0.2460155 0.3894278
[43] 0.1757617 0.1710472 0.2191764 0.2322201 0.3697992

$C.SR2
 [1]  0.25325581          NA -0.13173242          NA          NA          NA
 [7]  0.86538462          NA  0.29687625          NA  0.25863290  0.19297824
[13]          NA  0.75000000          NA          NA  0.40515202  0.07153903
[19]          NA          NA  0.02538134  0.27194107  0.14799928          NA
[25]          NA  0.13859506          NA          NA  0.13691492  0.66902355
[31]  0.13429489          NA  0.10743222  0.18067369          NA  0.64000000
[37]  0.14519975 -0.06666667          NA  0.11167650  1.00000000 -0.06576529
[43]          NA          NA  1.00000000          NA  0.06396005

$C.RMSE
 [1]  3.8289846  4.2426407  2.8054180  1.1952286  4.6904158  5.5901699
 [7] 21.5377635 11.0000000  2.4684240  1.2247449  8.5888239 18.6291965
[13] 53.6842621  2.4452861  1.5986105  1.2649111  1.2084911 12.7377889
[19]  1.4832397  6.4420494 29.8767353  8.0650526 16.4057429  1.9306146
[25]  2.8619008  2.5501895  1.0000000  7.5762788  8.5403186  3.7326446
[31]  3.5256052  4.6841410  3.8957210 32.8043011  1.0000000  3.6425548
[37] 26.0579962 11.6189500  5.6315438 12.1800394  0.9110262  2.4572232
[43] 69.0000000  9.2195445  4.0210677  2.2912878 84.5247484
preds_m_env <- computePredictedValues(m_env, expected = FALSE)
evaluateModelFit(hM = m_env, predY = preds_m_env)
$RMSE
 [1]  2.2741757  0.9185587  1.2183493  0.4145781  1.0155048  1.6393596
 [7]  7.1687255  4.1250000  1.1180340  0.3750000  6.4565858 10.8074410
[13] 14.9942697  0.5153882  0.5994789  0.3535534  0.6614378 10.7476378
[19]  0.4145781  1.5724682 16.7575776  4.7467094 11.0467755  0.8003905
[25]  1.6369751  1.5674920  0.2795085  2.1176343  5.1844871  1.1659224
[31]  2.2360680  2.2500000  2.3518610 27.1409803  0.3061862  1.0458250
[37] 13.7239811  2.9606798  1.8624581 10.8818226  0.3535534  1.0680005
[43]  8.6250000  1.5064963  0.8385255  0.5728220 42.9478245

$SR2
 [1]  0.067265023           NA           NA -0.001946117           NA
 [6]  0.016576588  0.463932917           NA  0.567369002           NA
[11]  0.475422030  0.320717512  0.048078189  0.677133106           NA
[16]           NA  0.231141642  0.439019589           NA  0.466655461
[21]  0.429837923  0.587319651  0.393356752           NA  0.000140046
[26]  0.200495848           NA           NA  0.477244396  0.500250523
[31]           NA  0.159558196  0.306250000  0.497289192           NA
[36]  0.761898570  0.468544710           NA           NA  0.363385529
[41]           NA  0.130232321           NA  0.250000000  0.999748047
[46]           NA  0.642030934

$O.AUC
 [1] 0.7657807 0.7650273 0.7860577 0.8020050 0.6830601 0.8738839 0.8847118
 [8] 0.8434343 0.9562594 0.6278736 0.8352941 0.9134199 0.8322034 0.8579235
[15] 0.6717172 0.5338983 0.8477564 0.9472097 0.5694915 0.9916667 0.8630682
[22] 0.9140625 0.8518519 0.7684391 0.7585825 0.7232955 0.5881356 0.9779661
[29] 0.8695887 0.8472222 0.6521739 0.9186483 0.8748616 0.9330357 0.8908046
[36] 0.8250000 0.8888889 0.9729167 0.7819549 0.8218623 0.5423729 0.8461538
[43] 0.9682540 0.9758065 1.0000000 0.9416667 0.9443054

$O.TjurR2
 [1]  0.115287929  0.065915301  0.086153846  0.094035088  0.093346995
 [6]  0.214687500  0.306096491  0.150217172  0.411048265  0.023836207
[11]  0.263367647  0.241590909  0.157372881  0.282964481  0.077782828
[16]  0.009855932  0.155480769  0.263371041  0.025059322  0.304750000
[21]  0.254840909  0.342578125  0.225305305  0.101042024  0.086683278
[26]  0.111534091  0.007584746  0.155635593  0.219063853  0.296129630
[31]  0.047680276  0.152859825  0.197956811  0.261964286  0.112945402
[36]  0.396958333  0.216799517  0.140875000  0.149085213  0.189099190
[41] -0.005059322  0.152804487  0.190357143  0.284596774  0.564838710
[46]  0.215500000  0.295619524

$O.RMSE
 [1] 0.4314218 0.2274015 0.3728042 0.3159321 0.2432111 0.3052175 0.2576981
 [8] 0.3248772 0.2885379 0.2991020 0.4099954 0.3822428 0.2913284 0.2238868
[15] 0.3332652 0.3053240 0.3594098 0.3159574 0.2835150 0.2213497 0.3791829
[22] 0.3731024 0.4129366 0.3511970 0.4405368 0.4285494 0.2706995 0.2564292
[29] 0.3887310 0.2971433 0.4655246 0.3908544 0.3984952 0.2687369 0.2756514
[36] 0.1918167 0.3848196 0.2501794 0.2918616 0.4356568 0.3048846 0.3524755
[43] 0.1871661 0.2012830 0.1741088 0.2347297 0.3466638

$C.SR2
 [1]  0.195348837           NA           NA           NA           NA
 [6]  0.437500000  0.715976331           NA  0.738033561           NA
[11]  0.327479056  0.420213604  0.125000000  0.250000000           NA
[16]           NA  0.071168843  0.140079855           NA -0.400000000
[21]  0.137381364  0.533785083  0.311856459           NA  0.009824102
[26]  0.283306752           NA           NA  0.056872571  0.317373161
[31]           NA  0.149962719  0.104685819  0.348511840           NA
[36]  0.640000000  0.412515496           NA           NA  0.417066885
[41]           NA -0.017007660           NA  1.000000000  1.000000000
[46]           NA  0.221145017

$C.RMSE
 [1]  3.874744  4.242641  2.813657  1.195229  4.690416  4.329483 21.161383
 [8] 11.000000  2.050051  1.224745  8.640500 18.013609 53.602610  1.951468
[15]  1.598611  1.264911  1.396435 11.747394  1.483240  6.178818 29.494733
[22]  6.420966 16.639738  1.930615  2.836833  2.638448  1.000000  7.576279
[29]  8.550549  2.648784  3.730019  4.094269  3.800918 28.335322  1.000000
[36]  3.563396 25.456573 11.842719  5.631544 13.244206  1.264911  2.466114
[43] 69.000000  7.864421  3.928963  2.291288 82.210428

También es posible realizar evaluación por validación cruzada

partition_env <- createPartition(m_env, nfolds = 4)
preds_m_env = computePredictedValues(m_env, 
                               partition = partition_env,
                               nParallel = nChains, 
                               expected = FALSE)
Cross-validation, fold 1 out of 4
Cross-validation, fold 2 out of 4
Cross-validation, fold 3 out of 4
Cross-validation, fold 4 out of 4
evaluateModelFit(hM = m_env, predY = preds_m_env)
$RMSE
 [1]  2.3692364  0.9270248  1.1989579  0.4841229  1.0155048  1.9921722
 [7]  7.1556359  3.9863047  1.4415703  0.3750000  7.4908799 11.0101161
[13] 15.0057281  0.6959705  0.6281172  0.3750000  0.5590170 12.4318454
[19]  0.4330127  1.6345871 16.9831028  6.3087860 12.5151471  0.8196798
[25]  1.6345871  1.7589059  0.2500000  2.1250000  5.1588062  1.6583124
[31]  2.3659367  2.4302778  2.3551805 33.5030316  0.3061862  1.5051993
[37] 13.8379291  2.9869508  1.8874586 14.1140721  0.3535534  1.1336473
[43]  8.6268114  1.6298006  1.1319231  0.5728220 41.0925499

$SR2
 [1]  0.0148978333 -0.0007801073  0.0006487137 -0.0039540153            NA
 [6] -0.0022574956  0.4664876950  0.0609968401  0.2204927550            NA
[11]  0.2533349598  0.2412122513 -0.0013423269 -0.0015853793 -0.0080011391
[16] -0.0013437474  0.1726759601  0.1387219803 -0.0013430368 -0.0021472000
[21]  0.0987028632  0.3960691696  0.1399746531 -0.0066310908 -0.0002148796
[26] -0.0007494076  0.1873015873 -0.0027282431  0.3289912540  0.1280791408
[31]  0.0061158414 -0.0297328735  0.1885815015  0.2795779046            NA
[36] -0.0021472000  0.1218893149 -0.0021472000 -0.0060121115  0.1180849557
[41]            NA -0.0235504317 -0.0005120328            NA  0.1431421385
[46]            NA  0.2400568677

$O.AUC
 [1] 0.60797342 0.43169399 0.52804487 0.58521303 0.38797814 0.70424107
 [7] 0.81328321 0.65656566 0.89064857 0.72270115 0.78431373 0.74512987
[13] 0.60000000 0.09562842 0.47575758 0.45593220 0.71153846 0.78054299
[19] 0.58305085 0.67916667 0.73125000 0.80712891 0.72472472 0.60634648
[25] 0.57419712 0.51136364 0.55932203 0.47288136 0.81493506 0.74629630
[31] 0.57529162 0.53441802 0.74197121 0.82254464 0.70258621 0.71041667
[37] 0.75181159 0.62708333 0.69674185 0.68572874 0.69322034 0.44471154
[43] 0.38095238 0.81451613 0.91129032 0.49375000 0.82102628

$O.TjurR2
 [1]  0.058842746 -0.026994536  0.008477564  0.026384712 -0.040095628
 [6]  0.077767857  0.241679198  0.089464646  0.328024133  0.087859195
[11]  0.185274510  0.168625541  0.016601695 -0.126857923 -0.026333333
[16]  0.001915254  0.175945513  0.183080694  0.004974576  0.042375000
[21]  0.144409091  0.221953125  0.131536537  0.046702401  0.028186600
[26]  0.006795455  0.058957627 -0.036093220  0.228360390  0.162120370
[31]  0.019711029  0.011799124  0.152981728  0.214062500  0.081091954
[36]  0.117125000  0.146117150  0.024625000  0.086942356  0.107667004
[41]  0.046483051 -0.041522436 -0.045595238  0.090120968  0.243629032
[46] -0.007041667  0.229105131

$O.RMSE
 [1] 0.4629301 0.2762097 0.4182818 0.3379512 0.2821512 0.3398078 0.3076775
 [8] 0.3715797 0.3412119 0.3124148 0.4412411 0.4336597 0.3298927 0.2940314
[15] 0.3832552 0.3131945 0.3575132 0.3682917 0.3063869 0.3016337 0.4441822
[22] 0.4301488 0.4629820 0.3832752 0.4755124 0.4875897 0.2998299 0.3347934
[29] 0.4018972 0.3555274 0.4969531 0.4536270 0.4341529 0.2994313 0.2891535
[36] 0.2783516 0.4236588 0.3047000 0.3362994 0.4727714 0.3028557 0.4406333
[43] 0.2321206 0.2340696 0.2363682 0.2886078 0.3928792

$C.SR2
 [1] -7.920865e-02            NA  1.008576e-01            NA            NA
 [6]            NA  6.523669e-01  5.047269e-01  2.287485e-01            NA
[11]  1.631133e-01  2.068807e-01            NA            NA            NA
[16]            NA  2.994989e-01  1.487011e-02            NA            NA
[21] -1.322377e-02  3.681269e-01 -1.192154e-05            NA  2.063061e-02
[26] -6.501338e-02            NA            NA  7.928032e-02  5.259160e-04
[31]  6.513959e-02            NA  8.725524e-02  1.811099e-01            NA
[36]            NA  1.248036e-01            NA            NA  4.435049e-02
[41]            NA            NA            NA            NA -1.000000e+00
[46]            NA  2.862973e-01

$C.RMSE
 [1]  4.1197201  4.2426407  2.6442242  1.1952286  4.6904158  5.5901699
 [7] 21.0988731 10.2553376  2.6372507  1.2247449 10.0675036 18.4885165
[13] 53.6842621  3.1091264  1.5986105  1.2649111  1.2507799 13.5975974
[19]  1.4832397  6.4420494 30.1115758  8.4895983 18.9323919  1.9306146
[25]  2.7966653  3.0598119  0.9636271  7.5762788  8.6106794  4.1500689
[31]  3.4040679  4.6841410  3.6909849 35.3374065  1.0000000  5.8094750
[37] 25.7846807 11.8427193  5.6315438 17.7649962  1.2649111  2.5331140
[43] 69.0000000  9.2195445  6.0723675  2.2912878 78.4833393

Una vez que tengamos un modelo sólido. Nos interesa ver que papel juegan las variables ambientales en las especies/comunidades.

postBeta <- getPostEstimate(m_env, parName = "Beta")
plotBeta(m_env, 
         post = postBeta, 
         param = "Mean", 
         supportLevel = 0.9)

Otro aspecto importante/novedoso es estimar las correlaciones entre especies

library(corrplot)
OmegaCor <- computeAssociations(m0)
supportLevel <- 0.95
toPlot <- ((OmegaCor[[1]]$support > supportLevel) +
           (OmegaCor[[1]]$support < (1-supportLevel)
            ) > 0) * OmegaCor[[1]]$mean
corrplot(toPlot)

library(corrplot)
OmegaCor <- computeAssociations(m_env)
supportLevel <- 0.95
toPlot <- ((OmegaCor[[1]]$support > supportLevel) +
           (OmegaCor[[1]]$support < (1-supportLevel)
            ) > 0) * OmegaCor[[1]]$mean
corrplot(toPlot)

Cuando especificamos un modelo sin variables ambientales, estamos haciendo un modelo “unconstrained” o indirecto. Por lo que se puede usar para generar gráficos de ordenación.

etaPost <- getPostEstimate(m_env, "Eta")
lambdaPost <- getPostEstimate(m_env, "Lambda")
biPlot(m_env, 
       etaPost = etaPost, 
       lambdaPost = lambdaPost, 
       factors = c(1,2), 
       "pH")

En este tipo de modelos, la importancia de cada variable se mide en términos de variación explicada y se calcula con análisis de particionado de la varianza

VP = computeVariancePartitioning(m_env)
plotVariancePartitioning(m_env, VP = VP)

Como en los Modelos de Distribución de Especies, podemos calcular curvas de respuesta de cada especie/comunidad a las variables ambientales

Gradient = constructGradient(m_env, focalVariable = "N-NO3-")
Gradient$XDataNew
       N-NO3-        N-NH3      N-NO2-        PO43-            pH           DO
1  -0.8670817 -0.124331499 -0.27542798  0.107613107 -0.0116785746  0.094728129
2  -0.6626868 -0.095023153 -0.21050205  0.082245745 -0.0089256141  0.072398110
3  -0.4582918 -0.065714807 -0.14557611  0.056878382 -0.0061726536  0.050068090
4  -0.2538969 -0.036406461 -0.08065018  0.031511020 -0.0034196931  0.027738071
5  -0.0495019 -0.007098115 -0.01572425  0.006143658 -0.0006667326  0.005408052
6   0.1548930  0.022210231  0.04920168 -0.019223704  0.0020862279 -0.016921968
7   0.3592880  0.051518577  0.11412762 -0.044591066  0.0048391884 -0.039251987
8   0.5636830  0.080826923  0.17905355 -0.069958429  0.0075921489 -0.061582006
9   0.7680779  0.110135269  0.24397948 -0.095325791  0.0103451094 -0.083912026
10  0.9724729  0.139443615  0.30890541 -0.120693153  0.0130980699 -0.106242045
11  1.1768678  0.168751962  0.37383134 -0.146060515  0.0158510304 -0.128572065
12  1.3812628  0.198060308  0.43875728 -0.171427877  0.0186039909 -0.150902084
13  1.5856577  0.227368654  0.50368321 -0.196795240  0.0213569514 -0.173232103
14  1.7900527  0.256677000  0.56860914 -0.222162602  0.0241099119 -0.195562123
15  1.9944476  0.285985346  0.63353507 -0.247529964  0.0268628724 -0.217892142
16  2.1988426  0.315293692  0.69846101 -0.272897326  0.0296158329 -0.240222161
17  2.4032375  0.344602038  0.76338694 -0.298264688  0.0323687934 -0.262552181
18  2.6076325  0.373910384  0.82831287 -0.323632050  0.0351217540 -0.284882200
19  2.8120274  0.403218730  0.89323880 -0.348999413  0.0378747145 -0.307212219
20  3.0164224  0.432527076  0.95816473 -0.374366775  0.0406276750 -0.329542239
            EC         NTU          T°C Current.velocity Average.depth
1  -0.40452804 -0.19248398 -0.064716439     0.0141457403   0.023325882
2  -0.30916968 -0.14711022 -0.049460999     0.0108112011   0.017827331
3  -0.21381132 -0.10173647 -0.034205558     0.0074766620   0.012328781
4  -0.11845296 -0.05636271 -0.018950118     0.0041421228   0.006830231
5  -0.02309460 -0.01098896 -0.003694677     0.0008075837   0.001331680
6   0.07226376  0.03438480  0.011560764    -0.0025269554  -0.004166870
7   0.16762211  0.07975855  0.026816204    -0.0058614946  -0.009665421
8   0.26298047  0.12513231  0.042071645    -0.0091960337  -0.015163971
9   0.35833883  0.17050606  0.057327085    -0.0125305729  -0.020662522
10  0.45369719  0.21587982  0.072582526    -0.0158651120  -0.026161072
11  0.54905555  0.26125357  0.087837967    -0.0191996512  -0.031659622
12  0.64441391  0.30662733  0.103093407    -0.0225341903  -0.037158173
13  0.73977226  0.35200109  0.118348848    -0.0258687294  -0.042656723
14  0.83513062  0.39737484  0.133604288    -0.0292032686  -0.048155274
15  0.93048898  0.44274860  0.148859729    -0.0325378077  -0.053653824
16  1.02584734  0.48812235  0.164115170    -0.0358723469  -0.059152374
17  1.12120570  0.53349611  0.179370610    -0.0392068860  -0.064650925
18  1.21656406  0.57886986  0.194626051    -0.0425414252  -0.070149475
19  1.31192241  0.62424362  0.209881491    -0.0458759643  -0.075648026
20  1.40728077  0.66961737  0.225136932    -0.0492105034  -0.081146576
   Average.width
1    0.101588442
2    0.077641259
3    0.053694075
4    0.029746892
5    0.005799708
6   -0.018147475
7   -0.042094658
8   -0.066041842
9   -0.089989025
10  -0.113936209
11  -0.137883392
12  -0.161830575
13  -0.185777759
14  -0.209724942
15  -0.233672126
16  -0.257619309
17  -0.281566492
18  -0.305513676
19  -0.329460859
20  -0.353408043
predY <- predict(m_env, 
                XData = Gradient$XDataNew, 
                studyDesign = Gradient$studyDesignNew,
                ranLevels = Gradient$rLNew, 
                expected = FALSE)
plotGradient(m_env, Gradient, pred=predY, measure="S")

[1] 0.815
plotGradient(m_env, Gradient, pred=predY, measure="Y", index = 34)

[1] 0.525