Script do R Módulo 10

Script limpo

Aqui apresento o scrip na íntegra sem os textos ou outros comentários. Você pode copiar e colar no R para executa-lo. Lembre de remover os # ou ## caso necessite executar essas linhas.

## rm(list=ls(all=TRUE)) ##LIMPA A MEMORIA
## cat("\014") #limpa o console
## dev.off() #apaga os graficos, se houver algum
## install.packages("openxlsx")
## install.packages("tidyverse")
## install.packages("vegan")
## install.packages("gplots")
## install.packages("psych")
## install.packages("gt")
## install.packages("ggplot2")
library(tidyverse)
## getwd()
## setwd("C:/Seu/Diretório/De/Trabalho")
library(openxlsx)
ppbio <- read.xlsx("D:/Elvio/OneDrive/Disciplinas/_EcoNumerica/5.Matrizes/ppbio06p.xlsx",
                   rowNames = T,
                   colNames = T,
                   sheet = "Sheet1")
ppbio_a <- read.xlsx("D:/Elvio/OneDrive/Disciplinas/_EcoNumerica/5.Matrizes/ppbio06h.xlsx",
                   rowNames = T,
                   colNames = T,
                   sheet = "Sheet1")
#ppbio
#str(ppbio)
#ppbio_ma <- as.matrix(ppbio) #lê ppbio como uma matrix
#ppbio_ma
#str(ppbio_ma)
## getwd()
## ppbio <- read.xlsx(file.choose(),
##                    rowNames = T, colNames = T,
##                    sheet = "Sheet1")
#Lista as colunas
colnames(ppbio)
#Escolher quais colunas usar por nome
colnames(ppbio)[rev(order(colSums(ppbio)))] #ordena por maior soma
#Usar a função subset()
m_part <- subset(ppbio_a[, c("a.veloc", "a.temp", "a.do", "a.transp")])
m_part <- subset(ppbio_a[, 18:26])
m_trab <- (ppbio)
#m_trab <- (m_part)
#m_trab
#Dendrograma
library(vegan)
m_trns <- asin(sqrt(decostand(m_trab,
                               method="total", MARGIN = 2)))
#m_trns <- sqrt(m_trab)
vegdist <- vegdist(m_trns, method = "bray",
                   diag = TRUE,
                   upper = FALSE)
cluster_uas <- hclust(vegdist, method = "average")
plot (cluster_uas, main = "Cluster Dendrogram - Bray-Curtis da Matriz Comunitária",
      hang = 0.1) #testar com -.01
rect.hclust(cluster_uas, k = 3, h = NULL) 
#h = 0.8 fornece os grupos formados na altura h
as.matrix(vegdist)[1:6, 1:6]
#Heatmap
library("gplots")
heatdist <- as.matrix(vegdist)
col <- rev(heat.colors(999)) #rev() reverte as cores do heatmap
heatmap.2(x=(as.matrix(vegdist)), #objetos x objetos
          Rowv = as.dendrogram(cluster_uas),
          Colv = as.dendrogram(cluster_uas),
          key = T, tracecol = NA, revC = T,
          col = heat.colors,  #dissimilaridade = 1 - similaridade
          density.info = "none",
          xlab = "UA´s", ylab = "UA´s",
          mar = c(6, 6) + 0.2)
cluster_spp <- hclust((vegdist(t(m_trns), method = "bray",
                            diag = TRUE,
                            upper = FALSE)), method = "average")
plot (cluster_spp, main = "Dendrograma dos atributos")
heatmap.2(t(as.matrix(m_trns)), #objetos x atributos
          Colv = as.dendrogram(cluster_uas),
          Rowv = as.dendrogram(cluster_spp),
          key = T, tracecol = NA, revC = T,
          col = col,
          density.info = "none",
          xlab = "Unidades amostrais", ylab = "Espécies",
          mar = c(6, 6) + 0.1)  # adjust margin size
library(gt)
merge <- as.data.frame(cluster_uas$merge)
merge[nrow(merge)+1,] = c("0","0")
height <- as.data.frame(round(cluster_uas$height, 2))
height[nrow(height)+1,] = c("1.0")
fusoes <- data.frame(Cluster_uas = merge, Height = height)
colnames(fusoes) <- c("Cluster1", "Cluster2", "Height")
UA <- rownames_to_column(as.data.frame(m_trns[, 0]))
colnames(UA) <- c("UAs")
No.UA <- 1:nrow(fusoes)
fusoes <- cbind(No.UA, UA, fusoes)
fusoes$Histórico <- 1:nrow(fusoes)
#fusoes
gt(fusoes)
#Dendrograma
library(vegan)
#m_trns <- asin(sqrt(decostand(m_trab,
#                               method="total", MARGIN = 2)))
m_trns <- sqrt(m_part)
vegdist <- vegdist(m_trns, method = "bray",
                   diag = TRUE,
                   upper = FALSE)
cluster_uas <- hclust(vegdist, method = "average")
plot (cluster_uas, main = "Cluster Dendrogram - Bray-Curtis para a Matriz Ambiental",
      hang = 0.1) #testar com -.01
rect.hclust(cluster_uas, k = 3, h = NULL) 
#h = 0.8 fornece os grupos formados na altura h
as.matrix(vegdist)[1:6, 1:6]
#Heatmap
library("gplots")
heatdist <- as.matrix(vegdist)
col <- rev(heat.colors(999)) #rev() reverte as cores do heatmap
heatmap.2(x=(as.matrix(vegdist)), #objetos x objetos
          Rowv = as.dendrogram(cluster_uas),
          Colv = as.dendrogram(cluster_uas),
          key = T, tracecol = NA, revC = T,
          col = heat.colors,  #dissimilaridade = 1 - similaridade
          density.info = "none",
          xlab = "UA´s", ylab = "UA´s",
          mar = c(6, 6) + 0.2)
cluster_spp <- hclust((vegdist(t(m_trns), method = "bray",
                            diag = TRUE,
                            upper = FALSE)), method = "average")
plot (cluster_spp, main = "Dendrograma dos atributos")
heatmap.2(t(as.matrix(m_trns)), #objetos x atributos
          Colv = as.dendrogram(cluster_uas),
          Rowv = as.dendrogram(cluster_spp),
          key = T, tracecol = NA, revC = T,
          col = col,
          density.info = "none",
          xlab = "Unidades amostrais", ylab = "Espécies",
          mar = c(6, 6) + 0.1)  # adjust margin size
library(gt)
merge <- as.data.frame(cluster_uas$merge)
merge[nrow(merge)+1,] = c("0","0")
height <- as.data.frame(round(cluster_uas$height, 2))
height[nrow(height)+1,] = c("1.0")
fusoes <- data.frame(Cluster_uas = merge, Height = height)
colnames(fusoes) <- c("Cluster1", "Cluster2", "Height")
UA <- rownames_to_column(as.data.frame(m_trns[, 0]))
colnames(UA) <- c("UAs")
No.UA <- 1:nrow(fusoes)
fusoes <- cbind(No.UA, UA, fusoes)
fusoes$Histórico <- 1:nrow(fusoes)
#fusoes
gt(fusoes)
#Dendrograma
library(vegan)
library(gplots)
# Dendrograma para as UA's da matriz comunitária
m_trns_c <- asin(sqrt(decostand(ppbio,
                               method="total", MARGIN = 2)))
vegdist_c <- vegdist(m_trns_c, method = "bray",
                   diag = TRUE,
                   upper = FALSE)
cluster_uas_c <- hclust(vegdist_c, method = "average")
plot (cluster_uas_c, main = "Cluster Dendrogram - Bray-Curtis da Matriz Comunitária",
      hang = 0.1) #testar com -.01
rect.hclust(cluster_uas_c, k = 3, h = NULL) 
#h = 0.8 fornece os grupos formados na altura h
as.matrix(vegdist_c)[1:6, 1:6]
# Dendrograma para as variáveis morfologicas da matriz ambiental
m_trns_a <- sqrt(m_part)
cluster_mfg_a <- hclust((vegdist(t(m_trns_a),
                            method = "bray",
                            diag = TRUE,
                            upper = FALSE)), method = "average")
plot (cluster_mfg_a, main = "Dendrograma dos atributos")
# Heatmap Comunidade vs. Morfologia
heatmap.2(t(as.matrix(m_trns_a)), #objetos x atributos
          Colv = as.dendrogram(cluster_uas_c),
          Rowv = as.dendrogram(cluster_mfg_a),
          key = T, tracecol = NA, revC = T,
          col = col,
          density.info = "none",
          xlab = "Unidades amostrais", ylab = "Morfologia",
          mar = c(6, 9) + 0.1)  # adjust margin size
plot(m_trns_a$"m.elev", m_trns_a$"m.river")
plot(scale(m_trns_a$"m.elev"), scale(m_trns_a$"m.river"))
plot((m_trns_a$"m.elev" - mean(m_trns_a$"m.elev")) / sd(m_trns_a$"m.elev"))
library(psych)
pairs.panels(m_trns_a[,1:9], 
             method = "pearson", # correlation method
             scale = FALSE, lm = FALSE,
             hist.col = "#00AFBB", pch = 19,
             density = TRUE,  # show density plots
             ellipses = TRUE, # show correlation ellipses
             alpha = 0.5
             )
hell <- decostand(ppbio, method = "hellinger")
m_trns_a <- sqrt(m_part)
rda <- rda(hell ~ ., m_trns_a) # variavel de resposta ~ variavel explanatória 
summary(rda)
str(rda)
modelo_rda <- rda$call
modelo_rda
#Gráfico triplot basico
plot(rda, choices = c(1,2))
plot (rda, display = c("sp","bp"))
#Ajustes no gráfico triplot
model <- ordiplot(rda, type = "none", scaling = 2, xlab = "RDA1", ylab = "RDA2") # type="points" ou "text"
points(rda, col = "black", pch = 21, bg = "gray", cex = 1)
text(rda, col = "black", cex = 0.7, pos = 1)
points(rda, dis = "bp", col = "red", pch = 21, cex = 1)
text(rda, dis = "bp", col = "red", cex = 0.9, pos = 2)
points(rda, dis = "sp", col = "blue", pch = 3, cex = 1)
text(rda, dis = "sp", col = "blue", cex = 0.8, pos = 1)
# Show only the first 5 species with higher values
#high_species <- order(abs(scores(rda, display = "species")[, 1]), decreasing = TRUE)[1:5]
# Get the indices of the first 5 species with higher values on the first axis
#text(rda, display = "species", select = high_species, col = "darkblue", pos = 4)  # Display the selected species
RsquareAdj(rda)
anova.cca(rda) #significância global
anova.cca(rda, by="axis") #sign. dos eixos da RDA
anova.cca(rda, by="terms") #sign. das var. explanatórias
#Função variance inflation factors
vif.cca(rda)
# remover valores com colinearidade maior que 20
# GLM
rda1 <- step(rda, scope = formula(rda), scale = 0, direction = "both", steps = 1000, test = "perm")
summary(rda1)
modelo_rda1 <- rda1$call
modelo_rda1
#|results: hold
#(matriz de dados comunitária ~ variáveis da matriz ambiental)
rda2 <- rda(formula = hell ~ m.elev + m.river + m.stream + m.distsource + m.distmouth, data = m_trns_a)
#(variáveis de resposta ~ variaveis explanatórias sugeridas pelo modelo stepwise regression)
summary(rda2)
str(rda2)
#Colinearidade
vif.cca(rda2)
#ANOVAS
anova(rda2, perm.max = 999)
anova(rda2, by="axis", perm.max = 999)
anova(rda2, by="terms", perm.max = 999)
anova(rda2, by="margin", perm.max = 999)
summary(rda2)
#Triplot
simplified_model <- ordiplot(rda2, type = "points", scaling = 2)