User Tools

Site Tools


en:expl_var_rscript
# Comparison of R2 and adjR2 on simulated data
# The simulation is run in parallel on 8 cores - if you don't have that many cores, decrease the value in "makeCluster" appropriately (e.g. cl <- makeCluster (4) for four cores)
 
library (simcom)  # devtools::install_github ('zdealveindy/simcom')
 
samp.vars <- as.matrix (expand.grid (samp = seq (10, 100, by = 10), vars = 1:5))
simcom <- simul.comm ()
 
library (parallel)
cl <- makeCluster (8)
clusterExport (cl, c('samp.vars', 'simcom'))
 
expl.var.0 <- parApply (cl, samp.vars, 1, FUN = function (sava)
{
  no.samp <- sava[1]
  no.vars <- sava[2]
  replicate (1000, expr =
  {
    com <- simcom:::sample.comm (simul.comm = simcom, Np = no.samp)
    rand.var <- matrix (rnorm (no.samp * no.vars), ncol = no.vars)
    rsq <- vegan:::RsquareAdj (vegan:::rda (com$a.mat ~ rand.var))
    rsq
  }, simplify = 'matrix')
})
 
r2.mean <- unlist (lapply (expl.var.0, FUN = function (x) mean (unlist(x[1,]))))
adj.r2.mean <- unlist (lapply (expl.var.0, FUN = function (x) mean (unlist(x[2,]))))
 
samp.vars.r2 <- as.data.frame (cbind (samp.vars, r2.mean, adj.r2.mean))
 
# relationship of R2 and number of samples
png ('R2-vs-adjR2-simulcomm.png', width = 8, height = 4, units = 'in', res = 300, pointsize = 8)
par (mfrow = c(1,2))
par (mar = c(5.1, 5.1, 4.1, 2.1))
plot (r2.mean ~ samp, data = samp.vars.r2[1:10,], cex = 2, pch = 21, col = 'black', bg = 'black', ylim = c(0,.15), type = 'b', xlab = list ('Number of samples', cex = 2), ylab = list (expression (italic (R)^2), cex = 2), las = 1, lwd = 2)
points (adj.r2.mean ~ samp, data = samp.vars.r2[1:10,], cex = 2, pch = 21, col = 'black', lwd = 2, bg = 'white', type = 'b', lty = 'dashed')
legend ('topright', legend = c(expression (italic (R)^2, italic (R)[adj]^2)), cex = 2, pch = c(21,21), col = c('black', 'black'), pt.bg = c('black', 'white'), lty = c('solid', 'dashed'), bty = 'n', lwd = 2)
 
 
plot (r2.mean ~ vars, data = samp.vars.r2[samp.vars.r2$samp == 100,], cex = 2, pch = 21, col = 'black', bg = 'black', ylim = c(0,.15), type = 'b', xlab = list ('Number of variables', cex = 2), ylab = list (expression (italic (R)^2), cex = 2), las = 1, lwd = 2)
points (adj.r2.mean ~ vars, data = samp.vars.r2[samp.vars.r2$samp == 100,], cex = 2, pch = 21, col = 'black', bg = 'white', type = 'b', lty = 'dashed', lwd = 2)
legend ('topright', legend = c(expression (italic (R)^2, italic (R)[adj]^2)), cex = 2, pch = c(21,21), col = c('black', 'black'), pt.bg = c('black', 'white'), lty = c('solid', 'dashed'), bty = 'n', lwd = 2)
dev.off ()
# comparison of variation explained by soil pH on Vltava data vs variation represented by the first axes of unconstrained ordination
 
library (vegan)
library (weimea)
data (vltava)
 
spe.hell <- decostand (log1p (vltava$spe), method = 'hell')
tbRDA <- rda (spe.hell ~ pH, vltava$env)
tbPCA <- rda (spe.hell)
 
tbRDA
tbPCA
 
tot.var <- tbPCA$CA$tot.chi
tbRDA_RDA1 <- tbRDA$CCA$eig[1]
tbRDA_PCA1 <- tbRDA$CA$eig[1]
tbPCA_PCA1 <- tbPCA$CA$eig[1]
tbPCA_PCA2 <- tbPCA$CA$eig[2]
 
png ('tbRDA-vs-tbPCA.png', width = 8, height = 4, units = 'in', res = 300)
par (mfrow = c(1,2))
ordiplot (tbRDA, xlab = paste ('RDA1 (', round (tbRDA_RDA1/tot.var*100, 1), '%)', sep = ''), ylab = paste ('PC1 (', round (tbRDA_PCA1/tot.var*100, 1), '%)', sep = ''), main = 'tb-RDA')
ordiplot (tbPCA, xlab = paste ('PC1 (', round (tbPCA_PCA1/tot.var*100, 1), '%)', sep = ''), ylab = paste ('PC2 (', round (tbPCA_PCA2/tot.var*100, 1), '%)', sep = ''), main = 'tb-PCA')
dev.off ()
en/expl_var_rscript.txt · Last modified: 2019/02/10 12:58 by David Zelený