relatorio.Rmd

---
title: "Relatório Estatística Computacional"
author: 
- Cristian Pessatti dos Anjos
- Thor Marinho
- Rafael Alisson Schipiura
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: pdf_document
header-includes: 
  - \renewcommand{\and}{\\}
bibliography: cite.bib
csl: associacao-brasileira-de-normas-tecnicas-note.csl
nocite: '@*'
---

```{r, echo = FALSE, include = FALSE}
library("tidyverse")
library("magrittr")

worms <- read_csv2("./worms.csv")
worms %<>% mutate(across(where(is.character), factor))
legenda <- read_csv2("./legenda.csv")
rownames(legenda) <- legenda$var
```

# Introdução

A base de dados usada no presente estudo buscou coletar
variáveis relacionadas às minhocas, ao solo, ao ambiente
e ao manejo avaliadas em 124 estudos publicados que
quantificaram a abundância de minhocas (> 7300 amostras)
em 765 locais com diferentes tipos de clima, solos, uso da terra
e sistemas de manejo no Brasil. O dataset propõe um conjunto de
variáveis que devem ser avaliadas ao estudar as comunidades de
minhocas e outros grupos de macrofauna para promover
o uso de minhocas como bioindicadores da qualidade do solo.
As variáveis resposta são 'densidade',
número de indivíduos por metro cúbico,
e 'massa', peso dos indíviduos normalizado para um metro cúbico.

# Análise Univariada
Analisando cada variável individualmente. Primeiro as variáveis quantitativas, todas contínuas. Mesmo a densidade (em indivíduos/m³) é uma variável contínua devido à normalização realizada, pois não se coleta um metro cúbico de solo para contar os indivíduos.

```{r, echo = FALSE, include = TRUE}
numeric_vars <- colnames(worms %>% select_if(is.numeric))
for (aplot in numeric_vars) {
	symplot <- ensym(aplot)
	n <- worms %>% select(all_of(symplot)) %>% drop_na(all_of(symplot)) %>% nrow()
	sturges <- log(n,2)+1
	print(worms %>% drop_na(all_of(symplot)) %>% ggplot(aes(x = !!symplot)) + geom_histogram(bins = sturges) +  labs(x = legenda[aplot,][2], y="frequência"))
}
```

E as variáveis qualitativas, todas nominais.

```{r, echo = FALSE, include = TRUE}
factor_vars <- colnames(worms %>% select_if(is.factor))
for (aplot in factor_vars) {
	symplot <- ensym(aplot)
	print(worms %>% drop_na(all_of(symplot)) %>% ggplot(aes(x = !!symplot)) + geom_bar() + theme(axis.text.x = element_text(angle = 90)) + labs(x = legenda[aplot,][2], y = "frequência"))
}
```

# Análise Bivariada

Analisando as variáveis comparadas umas com as outras. Primeiro quantitativas, em pares.

```{r, echo = FALSE, include = TRUE}
duas_numeric <- combn(numeric_vars, 2)
for (aplot in 1:ncol(duas_numeric)) {
	primeira <- duas_numeric[,aplot][1]
	segunda <- duas_numeric[,aplot][2]
	symprim <- sym(primeira)
	symseg <- sym(segunda)
	print(worms %>% drop_na(all_of(symprim), all_of(symseg)) %>% ggplot(aes(x = !!symprim, y = !!symseg)) + geom_point() + geom_smooth(method = "lm",formula = y ~ x) + labs(x = legenda[primeira,][2], y = legenda[segunda,][2]) )
}
```

E os pares de variáveis qualitativas.

```{r, echo = FALSE, include = TRUE}
duas_factor <- combn(factor_vars, 2)
for (aplot in 1:ncol(duas_factor)) {
	primeira <- duas_factor[,aplot][1]
	symprim <- sym(primeira)
	segunda <- duas_factor[,aplot][2]
	symseg <- sym(segunda)
	print(worms %>% drop_na(all_of(symprim), all_of(symseg)) %>% ggplot(aes(x = !!symprim, y = !!symseg)) + geom_count() + labs(x = legenda[primeira,][2], y = legenda[segunda,][2]))
}
```
E, por fim, combinações de uma variável quatitativa e outra qualitativa:


```{r, echo = FALSE, include = TRUE}
duas_mista <- expand_grid(factor_vars, numeric_vars)
for (aplot in 1:nrow(duas_mista)) {
	primeira <- duas_mista[aplot,][[1]]
	symprim <- sym(primeira)
	segunda <- duas_mista[aplot,][[2]]
	symseg <- sym(segunda)
	print(worms %>% drop_na(all_of(symprim), all_of(symseg)) %>% ggplot(aes(x = !!symprim, y = !!symseg)) + geom_boxplot() + theme(axis.text.x = element_text(angle = 90)) + labs(x = legenda[primeira,][2], y = legenda[segunda,][2]))
}
```

# Tabela de Correlação

```{r, echo = FALSE, include = TRUE, message=FALSE, warning=FALSE}
library(corrplot)
worms_numeric <- worms[, numeric_vars]
mycor<- function(x, ...){
  r<- apply(x, 2, function(j){
    apply(x, 2, function(i){
      as.numeric(cor.test(i,j, ...)$estimate)
    })
  })
  P<- apply(x, 2, function(j){
    apply(x, 2, function(i){
      as.numeric(cor.test(i,j, ...)$p.value)
    })
  })
  out<-c()
  out$P<- P
  out$r<- r
  return(out)
}

mycorDat <- mycor(worms_numeric)

corrplot(corr=mycorDat$r,
        p.mat = mycorDat$P,
        type="full", insig="pch", sig.level =.1, pch.cex = .9)

```

# Correlação das variáveis com a variável densidade

```{r, echo = FALSE, include = TRUE}
nomes_colunas_dens_corr <- c("Estatística", "Parâmetro", "p-valor", "Estimativa", "Valor nulo", "Alternativa", "Método", "Variável", "Intervalo de confiança")

dens_corr <- matrix(ncol = 9)
for(numeric_var in numeric_vars) {
  corr_result <- cor.test(worms_numeric$dens, worms_numeric[[numeric_var]])
  corr_result$data.name <- numeric_var
  dens_corr <- rbind(
    matrix(corr_result, ncol = 9),
    dens_corr
  )
}

dens_corr_df <- data.frame(dens_corr)
colnames(dens_corr_df) <- nomes_colunas_dens_corr
knitr::kable(dens_corr_df)
```

# Correlação das variáveis com a variável massa

```{r, echo = FALSE, include = TRUE}
nomes_colunas_massa_corr <- c("Estatística", "Parâmetro", "p-valor", "Estimativa", "Valor nulo", "Alternativa", "Método", "Variável", "Intervalo de confiança")

massa_corr <- matrix(ncol = 9)
for(numeric_var in numeric_vars) {
  corr_result <- cor.test(worms_numeric$massa, worms_numeric[[numeric_var]])
  corr_result$data.name <- numeric_var
  massa_corr <- rbind(
    matrix(corr_result, ncol = 9),
    massa_corr
  )
}

massa_corr_df <- data.frame(massa_corr)
colnames(massa_corr_df) <- nomes_colunas_massa_corr
knitr::kable(massa_corr_df)
```

# Teste de hipótese para verificar linearidade da variável em relação a variável dens
Ou seja, a hipótese nula está assumindo que não há relação linear das propostas variáveis
em contraste com a variável dens (Densidade de indivíduos de minhocas por cm³)

```{r, echo = FALSE, include = TRUE}
check_hypothesis <- function(alpha_value, df, t_statistic) {
  p.upper <- pt(abs(t_statistic), df = df, lower.tail = FALSE)
  p.lower <- pt(-abs(t_statistic), df = df, lower.tail = TRUE)
  p <- p.upper + p.lower
  return(p <= alpha_value)
}

hypothesis_result_dens_corr_df <- data.frame(`Variável` = character(), `Rejeita hipótese nula` = logical())

for (var in dens_corr_df$`Variável`) {
  df <- dens_corr_df[dens_corr_df$`Variável` == var,]$`Parâmetro`[[1]] - 2
  t_statistic <- dens_corr_df[dens_corr_df$`Variável` == var,]$`Estatística`[[1]]
  if (is.numeric(t_statistic)) {
    hypothesis_result_dens_corr <- check_hypothesis(alpha_value = 0.05, df = df, t_statistic = t_statistic)
    hypothesis_result_dens_corr_df[nrow(hypothesis_result_dens_corr_df) + 1,] <- c(var, hypothesis_result_dens_corr)
  }
}

knitr::kable(hypothesis_result_dens_corr_df)
```

# Teste de hipótese para verificar linearidade da variável em relação a variável massa
Ou seja, a hipótese nula está assumindo que não há relação linear das propostas variáveis
em contraste com a variável massa (Massa de minhocas em Kg/m³)

```{r, echo = FALSE, include = TRUE}
hypothesis_result_massa_corr_df <- data.frame(`Variável` = character(), `Rejeita hipótese nula` = logical())

for (var in massa_corr_df$`Variável`) {
  df <- massa_corr_df[massa_corr_df$`Variável` == var,]$`Parâmetro`[[1]] - 2
  t_statistic <- massa_corr_df[massa_corr_df$`Variável` == var,]$`Estatística`[[1]]
  if (is.numeric(t_statistic)) {
    hypothesis_result_massa_corr <- check_hypothesis(alpha_value = 0.05, df = df, t_statistic = t_statistic)
    hypothesis_result_massa_corr_df[nrow(hypothesis_result_massa_corr_df) + 1,] <- c(var, hypothesis_result_massa_corr)
  }
}

knitr::kable(hypothesis_result_massa_corr_df)
```
# Referências