rorschach <- read.csv("Rorschach.csv")
print(rorschach)
## X Fear Anger Depression Love Ambition Security
## 1 Bat 33 10 18 1 2 6
## 2 Bear 0 0 2 0 0 0
## 3 Blood 10 5 2 1 0 0
## 4 Boots 0 1 2 0 0 0
## 5 Bridge 1 0 0 0 0 0
## 6 Butterfly 0 2 1 26 5 18
## 7 Cave 7 0 13 1 4 2
## 8 Clouds 2 9 30 4 1 6
## 9 Fire 5 9 1 2 1 1
## 10 Fur 0 3 4 5 5 21
## 11 Hair 0 1 1 2 0 0
## 12 Island 0 0 0 1 0 0
## 13 Mask 3 2 6 2 2 3
## 14 Mountains 2 1 4 1 18 2
## 15 Rocks 0 4 2 1 2 2
## 16 Smoke 1 6 1 0 1 0
# 1. PREPARE DATA
# ----------------
# assume X has columns: Card, Fear, Anger, Depression, Love, Ambition, Security
rownames(rorschach) <- rorschach$X
# species table: four primary emotions
rorsch_species <- rorschach[, c("Fear","Anger","Depression","Love")]
# drop any rows (cards) with all zeros
rorsch_species <- rorsch_species[rowSums(rorsch_species) > 0, ]
# environmental table: secondary variables
rorsch_env <- rorschach[rownames(rorsch_species), c("Ambition","Security")]
# 2. CORRESPONDENCE ANALYSIS
# --------------------------
# run CA
rorsch_ca <- cca(rorsch_species)
# basic plot
plot(rorsch_ca, type = "n", xlim = c(-2, 2), ylim = c(-2, 2))
# add site labels
text(rorsch_ca, display = "sites", labels = rownames(rorsch_species), cex = 0.8)
# add points
points(rorsch_ca, pch = 21, col = "red", bg = "yellow", cex = 1.2)
# add species labels
text(rorsch_ca, display = "species", col = "blue", cex = 0.8)
# fit env vars
fit_ca <- envfit(rorsch_ca, rorsch_env, permutations = 999)
plot(fit_ca, col = "black", lwd = 3)
fit_ca # shows which vectors are significant
##
## ***VECTORS
##
## CA1 CA2 r2 Pr(>r)
## Ambition 0.94241 -0.33446 0.0853 0.627
## Security 0.99855 -0.05376 0.5889 0.070 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
# 3. DETRENDED CORRESPONDENCE ANALYSIS
# ------------------------------------
rorsch_dca <- decorana(rorsch_species)
plot(rorsch_dca, type = "n")
text(rorsch_dca, display = "sites", labels = rownames(rorsch_species), cex = 0.8)
points(rorsch_dca, pch = 21, col = "red", bg = "yellow", cex = 0.8)
text(rorsch_dca, display = "species", col = "blue", cex = 0.8)
# optional: tweak number of segments
rorsch_dca2 <- decorana(rorsch_species, mk = 5)
plot(rorsch_dca2, type = "n")
text(rorsch_dca2, display = "sites", labels = rownames(rorsch_species), cex = 0.8)
points(rorsch_dca2, pch = 21, col = "red", bg = "yellow", cex = 0.8)
text(rorsch_dca2, display = "species", col = "blue", cex = 0.8)
# 4. NON-METRIC MDS (Euclidean distance)
# --------------------------------------
# find a 2-D solution
rorsch_mds2 <- metaMDS(rorsch_species, k = 2, distance = "euclidean")
## Wisconsin double standardization
## Run 0 stress 0.1051141
## Run 1 stress 0.1051141
## ... New best solution
## ... Procrustes: rmse 2.828322e-06 max resid 4.821115e-06
## ... Similar to previous best
## Run 2 stress 0.08455158
## ... New best solution
## ... Procrustes: rmse 0.1841046 max resid 0.5189631
## Run 3 stress 0.1593132
## Run 4 stress 0.1051141
## Run 5 stress 0.08455158
## ... Procrustes: rmse 5.64791e-06 max resid 1.573309e-05
## ... Similar to previous best
## Run 6 stress 0.1051141
## Run 7 stress 0.1176256
## Run 8 stress 0.1176256
## Run 9 stress 0.08455158
## ... Procrustes: rmse 1.070532e-05 max resid 2.753626e-05
## ... Similar to previous best
## Run 10 stress 0.1176256
## Run 11 stress 0.08455159
## ... Procrustes: rmse 4.944646e-05 max resid 0.0001366444
## ... Similar to previous best
## Run 12 stress 0.08455158
## ... Procrustes: rmse 6.579794e-06 max resid 1.785269e-05
## ... Similar to previous best
## Run 13 stress 0.1051141
## Run 14 stress 0.1176256
## Run 15 stress 0.08455159
## ... Procrustes: rmse 2.02057e-05 max resid 5.416518e-05
## ... Similar to previous best
## Run 16 stress 0.1176256
## Run 17 stress 0.1051141
## Run 18 stress 0.1852348
## Run 19 stress 0.08455158
## ... Procrustes: rmse 6.146712e-06 max resid 1.700302e-05
## ... Similar to previous best
## Run 20 stress 0.1051141
## *** Best solution repeated 6 times
# quick base plot
plot(rorsch_mds2, type = "t", cex = 1.1)
# more control with ordiplot
# build the empty ordiplot
fig <- ordiplot(rorsch_mds2, type = "none", cex = 1.1)
# add the species labels
text(fig, "species", col = "red", cex = 1.1)
# add the site labels
text(fig, "sites", col = "blue", cex = 0.8)
# overlay envfit
fit_nmds <- envfit(rorsch_mds2, rorsch_env, permutations = 999)
plot(fit_nmds, col = "black", lwd = 3)
fit_nmds
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Ambition 0.23899 -0.97102 0.0148 0.849
## Security 0.99775 -0.06699 0.1025 0.460
## Permutation: free
## Number of permutations: 999
# 5. SURFACE FITTING WITH SPLINES
# -------------------------------
fig <- ordiplot(rorsch_mds2, type = "none", cex = 1.1, main = "NMDS – Rorschach Data")
text(fig, "species", col = "red", cex = 0.7)
text(fig, "sites", col = "black", cex = 0.7)
# draw env vectors
plot(fit_nmds)
# add spline surfaces
ordisurf(rorsch_mds2, rorsch_env$Ambition, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 0 total = 1
##
## REML score: 45.13682
ordisurf(rorsch_mds2, rorsch_env$Security, add = TRUE, col = "green4")
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 1.69 total = 2.69
##
## REML score: 50.51879
# 6. THREE-DIMENSIONAL NMDS (Euclidean)
# -------------------------------------
# compute 3D NMDS
rorsch_mds3 <- metaMDS(rorsch_species, k = 3, distance = "euclidean")
## Wisconsin double standardization
## Run 0 stress 0
## Run 1 stress 9.482384e-05
## ... Procrustes: rmse 0.008857267 max resid 0.01753358
## Run 2 stress 9.720553e-05
## ... Procrustes: rmse 0.01142188 max resid 0.03005223
## Run 3 stress 9.955727e-05
## ... Procrustes: rmse 0.0133849 max resid 0.03637471
## Run 4 stress 9.514302e-05
## ... Procrustes: rmse 0.01233941 max resid 0.0263258
## Run 5 stress 9.80779e-05
## ... Procrustes: rmse 0.01346913 max resid 0.03685275
## Run 6 stress 9.975342e-05
## ... Procrustes: rmse 0.009873619 max resid 0.02235394
## Run 7 stress 8.945628e-05
## ... Procrustes: rmse 0.01012406 max resid 0.01982804
## Run 8 stress 9.758098e-05
## ... Procrustes: rmse 0.01295526 max resid 0.03132428
## Run 9 stress 0.07304579
## Run 10 stress 0.0003645748
## ... Procrustes: rmse 0.01471358 max resid 0.0403137
## Run 11 stress 0.0003947278
## ... Procrustes: rmse 0.01478234 max resid 0.04167848
## Run 12 stress 9.85796e-05
## ... Procrustes: rmse 0.009441807 max resid 0.0213194
## Run 13 stress 0.0002221011
## ... Procrustes: rmse 0.01363597 max resid 0.03860673
## Run 14 stress 9.368478e-05
## ... Procrustes: rmse 0.009233189 max resid 0.02065845
## Run 15 stress 9.977209e-05
## ... Procrustes: rmse 0.01163461 max resid 0.03046311
## Run 16 stress 9.572783e-05
## ... Procrustes: rmse 0.009561074 max resid 0.02294383
## Run 17 stress 0.07304633
## Run 18 stress 0.07304593
## Run 19 stress 9.796775e-05
## ... Procrustes: rmse 0.01238378 max resid 0.03475706
## Run 20 stress 0.000353623
## ... Procrustes: rmse 0.01440625 max resid 0.03960154
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 13: stress < smin
## 3: stress ratio > sratmax
## Warning in metaMDS(rorsch_species, k = 3, distance = "euclidean"): stress is
## (nearly) zero: you may have insufficient data
# fit environmental vectors for dimensions 1–3
fit_mds3 <- envfit(rorsch_mds3, rorsch_env, permutations = 999, choices = 1:3)
# 3D ordination plot with arrows
pl3 <- ordiplot3d(rorsch_mds3, envfit = fit_mds3, pch = 19, main = "3D NMDS (Euclidean)")
points(pl3, "points", pch = 16, col = "red", cex = 0.7)
text(pl3, "arrows", col = "blue", pos = 3)
# add species labels in 3D
sp3 <- scores(rorsch_mds3, choices = 1:3, display = "species")
text(pl3$xyz.convert(sp3), rownames(sp3), cex = 0.7, xpd = TRUE)
# compare stress of 2D vs 3D
stressplot(rorsch_mds2, main = "Stressplot 2D")
stressplot(rorsch_mds3, main = "Stressplot 3D")
# 7. NMDS USING BRAY–CURTIS DISTANCE
# ----------------------------------
# default Bray–Curtis NMDS in 2 dimensions
rorsch_mds_bc <- metaMDS(rorsch_species, k = 2)
## Wisconsin double standardization
## Run 0 stress 0.1116559
## Run 1 stress 0.1570399
## Run 2 stress 0.1116559
## ... New best solution
## ... Procrustes: rmse 8.102761e-06 max resid 1.909446e-05
## ... Similar to previous best
## Run 3 stress 0.1111066
## ... New best solution
## ... Procrustes: rmse 0.01974439 max resid 0.06116106
## Run 4 stress 0.06688559
## ... New best solution
## ... Procrustes: rmse 0.1911017 max resid 0.5480701
## Run 5 stress 0.06688561
## ... Procrustes: rmse 0.0002634186 max resid 0.0006529957
## ... Similar to previous best
## Run 6 stress 0.06688553
## ... New best solution
## ... Procrustes: rmse 0.0001126197 max resid 0.0002927313
## ... Similar to previous best
## Run 7 stress 0.15704
## Run 8 stress 0.06688554
## ... Procrustes: rmse 2.592285e-05 max resid 6.879284e-05
## ... Similar to previous best
## Run 9 stress 0.06688555
## ... Procrustes: rmse 5.104569e-05 max resid 0.0001352803
## ... Similar to previous best
## Run 10 stress 0.1111066
## Run 11 stress 0.1111066
## Run 12 stress 0.2068189
## Run 13 stress 0.1372835
## Run 14 stress 0.15704
## Run 15 stress 0.06688559
## ... Procrustes: rmse 0.0001013386 max resid 0.0002659198
## ... Similar to previous best
## Run 16 stress 0.1111066
## Run 17 stress 0.1372835
## Run 18 stress 0.1116559
## Run 19 stress 0.06688555
## ... Procrustes: rmse 6.76235e-05 max resid 0.0001858452
## ... Similar to previous best
## Run 20 stress 0.1116559
## *** Best solution repeated 5 times
# quick scree of stress across dims 1–4
stress_vals <- sapply(1:4, function(d)
metaMDS(rorsch_species, k = d, distance = "bray")$stress
)
## Wisconsin double standardization
## Run 0 stress 0.2039544
## Run 1 stress 0.2021874
## ... New best solution
## ... Procrustes: rmse 0.01756611 max resid 0.05089031
## Run 2 stress 0.2513238
## Run 3 stress 0.3125744
## Run 4 stress 0.3289399
## Run 5 stress 0.5349923
## Run 6 stress 0.330129
## Run 7 stress 0.2030572
## Run 8 stress 0.5159673
## Run 9 stress 0.2018833
## ... New best solution
## ... Procrustes: rmse 0.01046816 max resid 0.02768177
## Run 10 stress 0.3029342
## Run 11 stress 0.3387221
## Run 12 stress 0.3029171
## Run 13 stress 0.3121038
## Run 14 stress 0.3216512
## Run 15 stress 0.3207417
## Run 16 stress 0.2513238
## Run 17 stress 0.3920812
## Run 18 stress 0.2543006
## Run 19 stress 0.3982969
## Run 20 stress 0.378103
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: stress ratio > sratmax
## 19: scale factor of the gradient < sfgrmin
## Wisconsin double standardization
## Run 0 stress 0.1116559
## Run 1 stress 0.06688553
## ... New best solution
## ... Procrustes: rmse 0.1887123 max resid 0.547947
## Run 2 stress 0.1111066
## Run 3 stress 0.06688557
## ... Procrustes: rmse 8.931332e-05 max resid 0.0002338569
## ... Similar to previous best
## Run 4 stress 0.06688563
## ... Procrustes: rmse 0.0001451382 max resid 0.0003548361
## ... Similar to previous best
## Run 5 stress 0.0668856
## ... Procrustes: rmse 0.0001157369 max resid 0.0003130171
## ... Similar to previous best
## Run 6 stress 0.06688564
## ... Procrustes: rmse 0.0001768978 max resid 0.0004626383
## ... Similar to previous best
## Run 7 stress 0.06688557
## ... Procrustes: rmse 6.39868e-05 max resid 0.0001770611
## ... Similar to previous best
## Run 8 stress 0.06688556
## ... Procrustes: rmse 6.816399e-05 max resid 0.0001659013
## ... Similar to previous best
## Run 9 stress 0.1570399
## Run 10 stress 0.06688554
## ... Procrustes: rmse 3.126223e-05 max resid 8.290704e-05
## ... Similar to previous best
## Run 11 stress 0.06688555
## ... Procrustes: rmse 5.755777e-05 max resid 0.0001559243
## ... Similar to previous best
## Run 12 stress 0.06688555
## ... Procrustes: rmse 2.199084e-05 max resid 5.566325e-05
## ... Similar to previous best
## Run 13 stress 0.1394338
## Run 14 stress 0.06688554
## ... Procrustes: rmse 1.989108e-05 max resid 5.419757e-05
## ... Similar to previous best
## Run 15 stress 0.06688565
## ... Procrustes: rmse 0.0001441273 max resid 0.0003739288
## ... Similar to previous best
## Run 16 stress 0.1111066
## Run 17 stress 0.1111066
## Run 18 stress 0.06688555
## ... Procrustes: rmse 7.684716e-05 max resid 0.000204921
## ... Similar to previous best
## Run 19 stress 0.1116559
## Run 20 stress 0.06688553
## ... Procrustes: rmse 3.727463e-05 max resid 9.932557e-05
## ... Similar to previous best
## *** Best solution repeated 13 times
## Wisconsin double standardization
## Run 0 stress 0.02868572
## Run 1 stress 0.02868656
## ... Procrustes: rmse 0.0004083708 max resid 0.001188444
## ... Similar to previous best
## Run 2 stress 0.02873228
## ... Procrustes: rmse 0.004220213 max resid 0.01285496
## Run 3 stress 0.02872498
## ... Procrustes: rmse 0.003677402 max resid 0.01118975
## Run 4 stress 0.02866305
## ... New best solution
## ... Procrustes: rmse 0.006842169 max resid 0.02092068
## Run 5 stress 0.02877989
## ... Procrustes: rmse 0.01401928 max resid 0.04301267
## Run 6 stress 0.02873238
## ... Procrustes: rmse 0.01106778 max resid 0.03389779
## Run 7 stress 0.02865938
## ... New best solution
## ... Procrustes: rmse 0.00131359 max resid 0.004013101
## ... Similar to previous best
## Run 8 stress 0.028702
## ... Procrustes: rmse 0.007406573 max resid 0.02266563
## Run 9 stress 0.02869406
## ... Procrustes: rmse 0.006667264 max resid 0.0203044
## Run 10 stress 0.02866153
## ... Procrustes: rmse 0.001861082 max resid 0.00569946
## ... Similar to previous best
## Run 11 stress 0.02867817
## ... Procrustes: rmse 0.004942178 max resid 0.01513631
## Run 12 stress 0.02870996
## ... Procrustes: rmse 0.00808319 max resid 0.02473073
## Run 13 stress 0.02867066
## ... Procrustes: rmse 0.003838365 max resid 0.01175131
## Run 14 stress 0.02866031
## ... Procrustes: rmse 0.000531133 max resid 0.001628663
## ... Similar to previous best
## Run 15 stress 0.02865973
## ... Procrustes: rmse 0.0002690586 max resid 0.0008197836
## ... Similar to previous best
## Run 16 stress 0.02871989
## ... Procrustes: rmse 0.008861811 max resid 0.02711497
## Run 17 stress 0.02867525
## ... Procrustes: rmse 0.004552367 max resid 0.01393937
## Run 18 stress 0.02876407
## ... Procrustes: rmse 0.01124805 max resid 0.03444532
## Run 19 stress 0.02865942
## ... Procrustes: rmse 0.000680353 max resid 0.00186622
## ... Similar to previous best
## Run 20 stress 0.02866039
## ... Procrustes: rmse 0.0005975065 max resid 0.001817526
## ... Similar to previous best
## *** Best solution repeated 6 times
## Wisconsin double standardization
## Run 0 stress 0.01261575
## Run 1 stress 0.01905545
## Run 2 stress 0.0127123
## ... Procrustes: rmse 0.01289386 max resid 0.04102407
## Run 3 stress 0.01190344
## ... New best solution
## ... Procrustes: rmse 0.06503753 max resid 0.2028914
## Run 4 stress 0.01112145
## ... New best solution
## ... Procrustes: rmse 0.02700836 max resid 0.07069003
## Run 5 stress 0.01903925
## Run 6 stress 0.01271357
## Run 7 stress 0.01120128
## ... Procrustes: rmse 0.004608802 max resid 0.01119371
## Run 8 stress 0.0190282
## Run 9 stress 0.01276499
## Run 10 stress 0.01121247
## ... Procrustes: rmse 0.007575308 max resid 0.02098673
## Run 11 stress 0.01904605
## Run 12 stress 0.01261145
## Run 13 stress 0.01108631
## ... New best solution
## ... Procrustes: rmse 0.005151215 max resid 0.01314577
## Run 14 stress 0.01904811
## Run 15 stress 0.01903931
## Run 16 stress 0.01328132
## Run 17 stress 0.01108844
## ... Procrustes: rmse 0.0007368829 max resid 0.002104249
## ... Similar to previous best
## Run 18 stress 0.01164383
## Run 19 stress 0.01263937
## Run 20 stress 0.01264341
## *** Best solution repeated 1 times
plot(1:4, stress_vals, type = "b", lwd = 2,
xlab = "Dimensions", ylab = "Stress",
main = "Stress vs Dimensions (Bray–Curtis)")
# ordination plot
# 1) make a blank ordination plot
plot(rorsch_mds_bc, type = "n", main = "NMDS (Bray–Curtis)")
# 2) extract coordinates
site_scores <- scores(rorsch_mds_bc, display = "sites")
species_scores <- scores(rorsch_mds_bc, display = "species")
# 3) add points
points(site_scores, pch = 21, bg = "yellow", cex = 1)
points(species_scores, pch = 19, cex = 1)
# 4) add labels
text(site_scores, labels = rownames(site_scores), col = "blue", cex = 0.8, pos = 3)
text(species_scores, labels = rownames(species_scores), col = "red", cex = 1.1, pos = 3)
# 5) overlay envfit
fit_bc <- envfit(rorsch_mds_bc, rorsch_env, permutations = 999)
plot(fit_bc, col = "black", lwd = 3)
# 8. SPLINE SURFACES ON BRAY–CURTIS ORDINATION
# --------------------------------------------
fig_bc2 <- ordiplot(rorsch_mds_bc, type = "none",
main = "NMDS (Bray–Curtis) with Splines")
text(fig_bc2, "species", col = "red", cex = 0.8)
text(fig_bc2, "sites", col = "black", cex = 0.8)
plot(fit_bc)
# add ordisurf contours
ordisurf(rorsch_mds_bc, rorsch_env$Ambition, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 0.0001 total = 1
##
## REML score: 45.13681
ordisurf(rorsch_mds_bc, rorsch_env$Security, add = TRUE, col = "green4")
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 1.42 total = 2.42
##
## REML score: 50.57659