R Under development (unstable) (2024-05-17 r86566 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(mefa) mefa 3.2-9 2024-05-19 > demo(mefa) demo(mefa) ---- ~~~~ > wait <- function(vign=0) { + if (vign!=0) { + if (vign==1) { + ANSWER <- readline("Do you want to open the vignette? (y/n) ") + if (substr(tolower(ANSWER), 1, 1) == "y") + vignette("mefa", package = "mefa")} + if (vign==2) { + ANSWER <- readline("Do you want to check out mefa website? (y/n) ") + if (substr(tolower(ANSWER), 1, 1) == "y") + mefaweb()} + } else { + ANSWER <- readline("Please press ENTER to continue ... ") + } + } > cat("## Demo for the 'mefa' package") ## Demo for the 'mefa' package > wait(1) Do you want to open the vignette? (y/n) > cat("## Load the package and the example data set") ## Load the package and the example data set > library(mefa) > data(dol.count, dol.samp, dol.taxa) > str(dol.count) 'data.frame': 297 obs. of 4 variables: $ samp : Factor w/ 24 levels "DQ1","DQ2","DQ3",..: 7 7 7 7 7 7 7 10 8 8 ... $ taxa : Factor w/ 32 levels "aacu","amin",..: 1 2 2 13 17 25 25 32 1 2 ... $ count: num 1 1 2 1 1 13 15 0 1 1 ... $ segm : Factor w/ 3 levels "fresh","broken",..: 1 2 1 2 2 1 2 3 1 2 ... > str(dol.samp) 'data.frame': 24 obs. of 2 variables: $ microhab: Factor w/ 4 levels "dead.wood","litter",..: 2 2 2 2 2 1 1 1 1 1 ... $ method : Factor w/ 2 levels "time","quadrat": 2 2 1 2 1 2 1 2 1 2 ... > str(dol.taxa) 'data.frame': 121 obs. of 4 variables: $ species: Factor w/ 121 levels "Acanthinula aculeata",..: 3 2 81 82 16 17 92 93 72 32 ... $ author : Factor w/ 65 levels "(Alder, 1830)",..: 26 53 42 19 44 47 31 14 47 43 ... $ familia: Factor w/ 21 levels "Aciculidae","Bradybaenidae",..: 1 1 14 14 6 6 17 17 17 5 ... $ size : num 3.4 5.5 16 16 2.2 2.3 17 8 12 7.5 ... > wait() Please press ENTER to continue ... > cat("## Object classes") ## Object classes > cat("## 'stcs'") ## 'stcs' > x1 <- stcs(dol.count) > str(x1) Classes 'stcs' and 'data.frame': 297 obs. of 4 variables: $ samp : Factor w/ 24 levels "DQ1","DQ2","DQ3",..: 1 1 1 1 1 1 1 2 2 2 ... $ taxa : Factor w/ 29 levels "aacu","amin",..: 2 2 9 11 14 28 28 2 2 3 ... $ count: num 2 2 2 2 1 2 7 1 3 1 ... $ segm : Factor w/ 3 levels "fresh","broken",..: 1 2 1 1 1 2 1 2 1 2 ... - attr(*, "call")= language stcs(dframe = dol.count) - attr(*, "expand")= logi FALSE - attr(*, "zero.count")= logi TRUE - attr(*, "zero.pseudo")= chr "zero.pseudo" > unique(x1$count) [1] 2 1 7 3 13 4 5 15 8 10 16 17 25 0 6 9 19 12 22 26 > wait() Please press ENTER to continue ... > x2 <- stcs(dol.count, expand = TRUE) > str(x2) Classes 'stcs' and 'data.frame': 732 obs. of 4 variables: $ samp : Factor w/ 24 levels "DQ1","DQ2","DQ3",..: 1 1 1 1 1 1 1 1 1 1 ... $ taxa : Factor w/ 29 levels "aacu","amin",..: 2 2 2 2 9 9 11 11 14 28 ... $ count: num 1 1 1 1 1 1 1 1 1 1 ... $ segm : Factor w/ 3 levels "fresh","broken",..: 1 1 2 2 1 1 1 1 1 1 ... - attr(*, "call")= language stcs(dframe = dol.count, expand = TRUE) - attr(*, "expand")= logi TRUE - attr(*, "zero.count")= logi TRUE - attr(*, "zero.pseudo")= chr "zero.pseudo" > sum(x2$count) [1] 731 > unique(x2$count) [1] 1 0 > wait() Please press ENTER to continue ... > x3 <- stcs(dol.count, drop.zero = TRUE) > str(x3) Classes 'stcs' and 'data.frame': 296 obs. of 4 variables: $ samp : Factor w/ 23 levels "DQ1","DQ2","DQ3",..: 7 7 7 7 7 7 7 8 8 8 ... $ taxa : Factor w/ 28 levels "aacu","amin",..: 1 2 2 11 15 23 23 1 2 11 ... $ count: num 1 1 2 1 1 13 15 1 1 1 ... $ segm : Factor w/ 2 levels "fresh","broken": 1 2 1 2 2 1 2 1 2 1 ... - attr(*, "call")= language stcs(dframe = dol.count, drop.zero = TRUE) - attr(*, "expand")= logi FALSE - attr(*, "zero.count")= logi FALSE - attr(*, "zero.pseudo")= chr "not.defined" > unique(x3$count) [1] 1 2 13 15 5 8 3 10 16 17 25 7 4 12 22 26 6 9 19 > wait() Please press ENTER to continue ... > cat("## 'mefa'") ## 'mefa' > m1 <- mefa(x1) > m1 An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 2 (non-nested) segments: fresh, broken, $ samp: table for samples not provided, $ taxa: table for taxa not provided. > m1$xtab["LT1", ] aacu amin apur bbip bcan ccer clam cort ctri dbre dper druf eful estr ffau hobv 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 iiso mbor mobs odol ogla pinc ppyg pvic tuni vcos vicr vidi 0 0 0 0 0 0 0 0 0 0 0 0 > wait() Please press ENTER to continue ... > str(m1$xtab) num [1:24, 1:28] 0 0 0 0 0 0 1 1 0 0 ... - attr(*, "dimnames")=List of 2 ..$ samp: chr [1:24] "DQ1" "DQ2" "DQ3" "DT1" ... ..$ taxa: chr [1:28] "aacu" "amin" "apur" "bbip" ... > str(m1$segm) List of 2 $ fresh : num [1:24, 1:28] 0 0 0 0 0 0 1 1 0 0 ... ..- attr(*, "dimnames")=List of 2 .. ..$ samp: chr [1:24] "DQ1" "DQ2" "DQ3" "DT1" ... .. ..$ taxa: chr [1:28] "aacu" "amin" "apur" "bbip" ... $ broken: num [1:24, 1:28] 0 0 0 0 0 0 0 0 0 0 ... ..- attr(*, "dimnames")=List of 2 .. ..$ samp: chr [1:24] "DQ1" "DQ2" "DQ3" "DT1" ... .. ..$ taxa: chr [1:28] "aacu" "amin" "apur" "bbip" ... > wait() Please press ENTER to continue ... > mefa(x1, nested = TRUE) An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 2 nested segments: fresh, fresh-broken, $ samp: table for samples not provided, $ taxa: table for taxa not provided. > wait() Please press ENTER to continue ... > mefa(x1, segment = FALSE) An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 1 (all inclusive) segment, $ samp: table for samples not provided, $ taxa: table for taxa not provided. > wait() Please press ENTER to continue ... > mefa(x1, dol.samp) An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 2 (non-nested) segments: fresh, broken, $ samp: table for samples provided (2 variables), $ taxa: table for taxa not provided. > mefa(x1, , dol.taxa) An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 2 (non-nested) segments: fresh, broken, $ samp: table for samples not provided, $ taxa: table for taxa provided (4 variables). > wait() Please press ENTER to continue ... > m2 <- mefa(x1, dol.samp, dol.taxa) > m2 An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 2 (non-nested) segments: fresh, broken, $ samp: table for samples provided (2 variables), $ taxa: table for taxa provided (4 variables). > str(m2$xtab) num [1:24, 1:28] 0 0 0 0 0 0 1 1 0 0 ... - attr(*, "dimnames")=List of 2 ..$ samp: chr [1:24] "DQ1" "DQ2" "DQ3" "DT1" ... ..$ taxa: chr [1:28] "aacu" "amin" "apur" "bbip" ... > str(m2$samp) 'data.frame': 24 obs. of 2 variables: $ microhab: Factor w/ 4 levels "dead.wood","litter",..: 1 1 1 1 1 1 2 2 2 2 ... $ method : Factor w/ 2 levels "time","quadrat": 2 2 2 1 1 1 2 2 2 1 ... > str(m2$taxa) 'data.frame': 28 obs. of 4 variables: $ species: Factor w/ 121 levels "Acanthinula aculeata",..: 1 4 6 10 14 35 37 38 17 40 ... $ author : Factor w/ 65 levels "(Alder, 1830)",..: 42 57 1 41 27 49 41 38 47 15 ... $ familia: Factor w/ 21 levels "Aciculidae","Bradybaenidae",..: 18 21 21 4 4 4 4 4 6 21 ... $ size : num 2 9 5 18 18 18 17 13 2.3 NA ... > wait() Please press ENTER to continue ... > m2.sub <- mefa(x1, dol.samp[-c(1:5), ], dol.taxa[-c(1:80), ], xtab.fixed = FALSE) > m2.sub An object of class 'mefa' containing $ xtab: 228 individuals of 9 taxa in 19 samples, $ segm: 2 (non-nested) segments: fresh, broken, $ samp: table for samples provided (2 variables), $ taxa: table for taxa provided (4 variables). > str(m2.sub$xtab) num [1:19, 1:9] 4 4 3 1 2 3 0 3 7 13 ... - attr(*, "dimnames")=List of 2 ..$ samp: chr [1:19] "DQ1" "DQ2" "DQ3" "DT1" ... ..$ taxa: chr [1:9] "amin" "apur" "estr" "ffau" ... > str(m2.sub$samp) 'data.frame': 19 obs. of 2 variables: $ microhab: Factor w/ 4 levels "dead.wood","litter",..: 1 1 1 1 1 1 2 4 4 4 ... $ method : Factor w/ 2 levels "time","quadrat": 2 2 2 1 1 1 1 2 2 2 ... > str(m2.sub$taxa) 'data.frame': 9 obs. of 4 variables: $ species: Factor w/ 121 levels "Acanthinula aculeata",..: 4 6 48 24 51 57 76 79 99 $ author : Factor w/ 65 levels "(Alder, 1830)",..: 57 1 13 48 42 55 42 53 15 $ familia: Factor w/ 21 levels "Aciculidae","Bradybaenidae",..: 21 21 12 12 12 12 12 12 12 $ size : num 9 5 18 20 15 11 16 15.5 8 > wait() Please press ENTER to continue ... > mefalogo() > wait() Please press ENTER to continue ... > cat("## S3 methods") ## S3 methods > dim(m2) [1] 24 28 2 > dimnames(m2) $samp [1] "DQ1" "DQ2" "DQ3" "DT1" "DT2" "DT3" "LQ1" "LQ2" "LQ3" "LT1" "LT2" "LT3" [13] "RQ1" "RQ2" "RQ3" "RT1" "RT2" "RT3" "WQ1" "WQ2" "WQ3" "WT1" "WT2" "WT3" $taxa [1] "aacu" "amin" "apur" "bbip" "bcan" "ccer" "clam" "cort" "ctri" "dbre" [11] "dper" "druf" "eful" "estr" "ffau" "hobv" "iiso" "mbor" "mobs" "odol" [21] "ogla" "pinc" "ppyg" "pvic" "tuni" "vcos" "vicr" "vidi" $segm [1] "fresh" "broken" > wait() Please press ENTER to continue ... > summary(m2) Call: mefa(xtab = x1, samp = dol.samp, taxa = dol.taxa) Summary Total sum 731 Matrix fill (%) 26 Number of samples 24 Number of taxa 28 Number of segments 2 Segments (non-nested): fresh, broken s.rich s.abu t.occ t.abu Min. 0.000000 0.00000 1.000000 1.00000 1st Qu. 5.000000 9.75000 2.000000 4.00000 Median 7.000000 20.50000 4.500000 10.50000 Mean 7.333333 30.45833 6.285714 26.10714 3rd Qu. 10.000000 44.75000 10.000000 27.75000 Max. 16.000000 97.00000 23.000000 245.00000 > wait() Please press ENTER to continue ... > names(summary(m2)) [1] "s.rich" "s.abu" "t.occ" "t.abu" "ntot" [6] "mfill" "nsamp" "ntaxa" "nsegm" "segment" [11] "call" "nested" "drop.zero" "xtab.fixed" > summary(m2)$s.rich DQ1 DQ2 DQ3 DT1 DT2 DT3 LQ1 LQ2 LQ3 LT1 LT2 LT3 RQ1 RQ2 RQ3 RT1 RT2 RT3 WQ1 WQ2 5 7 4 7 10 5 5 5 7 0 2 5 10 16 13 11 15 10 9 10 WQ3 WT1 WT2 WT3 8 5 4 3 > summary(m2)$mfill [1] 0.2619048 > wait() Please press ENTER to continue ... > opar <- par(mfrow = c(1, 2)) > plot(m2, 1, main="A") > plot(m2, 4, type="rank", trafo="log", main="B") > par(opar) > wait() Please press ENTER to continue ... > opar <- par(mfrow = c(1, 2)) > boxplot(m2, 2, main = "A") > boxplot(m2, 3, main = "B") > par(opar) > wait() Please press ENTER to continue ... > molten <- melt(m2, "method") > str(molten) Classes 'stcs' and 'data.frame': 177 obs. of 4 variables: $ samp : Factor w/ 24 levels "DQ1","DQ2","DQ3",..: 1 1 1 1 1 2 2 2 2 2 ... $ taxa : Factor w/ 29 levels "aacu","amin",..: 2 9 11 14 28 2 3 9 11 12 ... $ count: num 4 2 2 1 9 4 1 20 2 2 ... $ segm : Factor w/ 3 levels "time","quadrat",..: 2 2 2 2 2 2 2 2 2 2 ... - attr(*, "call")= language melt.mefa(x = m2, segm.var = "method") - attr(*, "expand")= logi FALSE - attr(*, "zero.count")= logi TRUE - attr(*, "zero.pseudo")= chr "zero.pseudo" > wait() Please press ENTER to continue ... > m3 <- mefa(molten, dol.samp, dol.taxa) > m3 An object of class 'mefa' containing $ xtab: 731 individuals of 28 taxa in 24 samples, $ segm: 2 (non-nested) segments: time, quadrat, $ samp: table for samples provided (2 variables), $ taxa: table for taxa provided (4 variables). > wait() Please press ENTER to continue ... > opar <- par(mfrow = c(1, 3)) > image(m3, trafo = "log", sub = "all segments", main="A") > for (i in 1:2) + image(m3, segm = i, trafo = "log", + sub = dimnames(m3)$segm[i], main = LETTERS[i + 1]) > par(opar) > wait() Please press ENTER to continue ... > ex1 <- m2[1:20, 11:15, "fresh"] > dim(ex1) [1] 20 5 1 > dim(ex1$samp) [1] 20 2 > dim(ex1$taxa) [1] 5 4 > wait() Please press ENTER to continue ... > ex2 <- m2[m2$samp$method == "time"] > levels(ex2$samp$method) [1] "time" "quadrat" > ex3 <- m2[m2$samp$method == "time", drop = TRUE] > levels(ex3$samp$method) [1] "time" > wait() Please press ENTER to continue ... > size.5 <- as.factor(is.na(m3$taxa$size) | m3$taxa$size < 5) > levels(size.5) <- c("large", "small") > m4 <- aggregate(m3, "microhab", size.5) > t(m4$xtab) dead.wood litter live.wood rock large 53 33 80 210 small 50 112 99 94 > lapply(m4$segm, t) $time dead.wood litter live.wood rock large 32 12 16 118 small 5 0 1 3 $quadrat dead.wood litter live.wood rock large 21 21 64 92 small 45 112 98 91 > wait() Please press ENTER to continue ... > cat("## Writing reports") ## Writing reports > #set.seed(1234) > #m5 <- m2[ , sample(1:dim(m2)[2], 10)] > #report(m5, "report.tex", tex = TRUE, > # segment = TRUE, taxa.name = 1, author.name = 2, > # drop.redundant = 1) > > cat("## see mefadocs(\"SampleReport\") also") ## see mefadocs("SampleReport") also > wait() Please press ENTER to continue ... > cat("## Data analysis") ## Data analysis > mod.amin <- glm(m2$xtab[, "amin"] ~ ., + data = m2$samp, family = poisson) > summary(mod.amin) Call: glm(formula = m2$xtab[, "amin"] ~ ., family = poisson, data = m2$samp) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.6475 0.2858 2.266 0.02346 * microhablitter -0.3483 0.3770 -0.924 0.35559 microhablive.wood 0.3448 0.3170 1.088 0.27668 microhabrock 0.6633 0.2985 2.222 0.02630 * methodquadrat 0.6758 0.2281 2.963 0.00305 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 49.377 on 23 degrees of freedom Residual deviance: 28.515 on 19 degrees of freedom AIC: 106.99 Number of Fisher Scoring iterations: 5 > wait() Please press ENTER to continue ... > library(MASS) > mod.abu <- glm.nb(summary(m2)$s.abu ~ .^2, + data = m2$samp) > summary(mod.abu) Call: glm.nb(formula = summary(m2)$s.abu ~ .^2, data = m2$samp, init.theta = 4.730555855, link = log) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.5123 0.3122 8.046 8.54e-16 *** microhablitter -1.1260 0.5013 -2.246 0.02469 * microhablive.wood -0.7777 0.4762 -1.633 0.10245 microhabrock 1.1849 0.4198 2.823 0.00476 ** methodquadrat 0.5787 0.4279 1.352 0.17622 microhablitter:methodquadrat 1.8267 0.6441 2.836 0.00457 ** microhablive.wood:methodquadrat 1.6756 0.6237 2.687 0.00722 ** microhabrock:methodquadrat -0.1650 0.5812 -0.284 0.77643 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for Negative Binomial(4.7306) family taken to be 1) Null deviance: 94.216 on 23 degrees of freedom Residual deviance: 26.815 on 16 degrees of freedom AIC: 198.21 Number of Fisher Scoring iterations: 1 Theta: 4.73 Std. Err.: 1.71 2 x log-likelihood: -180.214 > wait() Please press ENTER to continue ... > prop.fr <- cbind(summary(m2[ , , "fresh"])$s.abu, summary(m2)$s.abu) > mod.fr <- glm(prop.fr ~ .^2, + data = m2$samp, family = binomial) > summary(mod.fr) Call: glm(formula = prop.fr ~ .^2, family = binomial, data = m2$samp) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.14518 0.24141 -0.601 0.547573 microhablitter -0.03714 0.49154 -0.076 0.939771 microhablive.wood -0.49081 0.47772 -1.027 0.304230 microhabrock -1.09526 0.30840 -3.551 0.000383 *** methodquadrat -0.35559 0.31373 -1.133 0.257036 microhablitter:methodquadrat -0.16278 0.55175 -0.295 0.767977 microhablive.wood:methodquadrat 0.29843 0.53561 0.557 0.577402 microhabrock:methodquadrat 1.17404 0.38609 3.041 0.002359 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 33.970 on 22 degrees of freedom Residual deviance: 13.825 on 15 degrees of freedom AIC: 117.2 Number of Fisher Scoring iterations: 4 > wait() Please press ENTER to continue ... > if (require(vegan)) { + m6 <- m2[summary(m2)$s.abu != 0, , ] + m6.ado <- adonis2(m6$xtab ~ .^2, + data = m6$samp, permutations = 100) + m6.ado + } Loading required package: vegan Loading required package: permute Loading required package: lattice This is vegan 2.6-6 Permutation test for adonis under reduced model Terms added sequentially (first to last) Permutation: free Number of permutations: 100 adonis2(formula = m6$xtab ~ .^2, data = m6$samp, permutations = 100) Df SumOfSqs R2 F Pr(>F) microhab 3 1.6124 0.24614 2.9314 0.019802 * method 1 1.1796 0.18007 6.4336 0.009901 ** microhab:method 3 1.0086 0.15397 1.8337 0.019802 * Residual 15 2.7502 0.41983 Total 22 6.5509 1.00000 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > wait() Please press ENTER to continue ... > if (require(vegan)) { + m2.cca <- cca(m2$segm[["fresh"]] ~ ., data=m2$samp, + subset=rowSums(m2$segm[["fresh"]]) > 0) + plot(m2.cca) + } > wait() Please press ENTER to continue ... > m.list <- list() > n1 <- rep(c("time", "quadrat"), each = 2) > n2 <- rep(c("fresh", "broken"), 2) > n3 <- paste(n1, n2, sep=".") > for (i in 1:4) { + m.list[[n3[i]]] <- + aggregate(m2[m2$samp$method == n1[i], , n2[i]], "microhab") + } > wait() Please press ENTER to continue ... > opar <- par(mfrow=c(2, 2)) > for (i in 1:4) { + tmp <- hclust(dist(m.list[[i]]$xtab), "ward") + plot(tmp, main = LETTERS[i], sub = names(m.list)[i], xlab = "") + } The "ward" method has been renamed to "ward.D"; note new "ward.D2" The "ward" method has been renamed to "ward.D"; note new "ward.D2" The "ward" method has been renamed to "ward.D"; note new "ward.D2" The "ward" method has been renamed to "ward.D"; note new "ward.D2" > par(opar) > cat("## End of mefa demo\n") ## End of mefa demo > wait(2) Do you want to check out mefa website? (y/n) > > proc.time() user system elapsed 1.23 0.23 1.45