Aquacultura - Dados de rotiferos

Introdução

Scrip com funções simples em R para análise exploratória de dados das aulas PL de Pescas e Aquacultura.

Carregar dados de Rotiferos das PLs

dados.rot <- read.table("PL_dados_rotiferos_R.csv",header=TRUE,sep=",",dec=".",na.strings="NA")

Ver os dados

head(dados.rot, 5)
##   Turma      Data Grupo  Hora Amostra_obs_ml Diluicao Estirpe Rotiferos_N Obs
## 1  PL-4 20-Feb-24     5 10:30            0.5     1:05       S          77    
## 2  PL-4 20-Feb-24     5 10:30            0.5     1:05       S          78    
## 3  PL-4 20-Feb-24     5 10:30            0.5     1:05       S          73    
## 4  PL-4 20-Feb-24     5 10:30            0.5     1:05       L          87    
## 5  PL-4 20-Feb-24     5 10:30            0.5     1:05       L          35

Calcular os rotiferos obs por 1ml

dados.rot$Rotiferos_obs_1ml <- dados.rot$Rotiferos_N/dados.rot$Amostra_obs_ml

calcular os rotiferos totais para as várias diluições

dados.rot$Rotiferos_total_ml <- rep(0, nrow(dados.rot)) #criar a nova variavel (placeholder)
for (i in 1:nrow(dados.rot)){
  if (dados.rot$Diluicao[i]=="1:05") dados.rot$Rotiferos_total_ml[i] <- dados.rot$Rotiferos_obs_1ml[i]*5 
  if (dados.rot$Diluicao[i]=="1:02") dados.rot$Rotiferos_total_ml[i] <- dados.rot$Rotiferos_obs_1ml[i]*2 
}

Sumarios dos dados

tabelas sumario dos dados, para cada estirpe

library(reshape2)
aggregate(Rotiferos_total_ml ~ Estirpe, data =dados.rot, FUN=mean)
##   Estirpe Rotiferos_total_ml
## 1       L           798.3407
## 2       S           488.0709
aggregate(Rotiferos_total_ml ~ Estirpe, data =dados.rot, FUN=min)
##   Estirpe Rotiferos_total_ml
## 1       L                155
## 2       S                140
aggregate(Rotiferos_total_ml ~ Estirpe, data =dados.rot, FUN=max)
##   Estirpe Rotiferos_total_ml
## 1       L               2640
## 2       S               1368

Graficos exploratórios

Boxplots para descrever os dados e detectar outliers

boxplot(Rotiferos_total_ml ~ Estirpe, data=dados.rot)

Gráfico mais completo, usando a library ggplot2

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
ggplot(dados.rot, aes(x=Estirpe, y=Rotiferos_total_ml, fill=Estirpe)) + 
  geom_boxplot()+
  ylab("Rotiferos (ml)")

Exemplo de gráfico para visualizar as contagens da duas estirpes, em cada uma das turmas PL

ggplot(dados.rot, aes(x=Turma, y=Rotiferos_total_ml, fill=Turma)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8)+
  ylab("Rotiferos (ml)")+
  theme(legend.position="bottom")+
  facet_grid(~Estirpe)

Pescas - Dados de maturação

Dados

Carregar dados de Maturacao de peixes

dados.mat <- read.table("PAq_Dados_Maturacao_R.csv",header=TRUE,sep=",",dec=".",na.strings="NA")

Ver os dados

head(dados.mat, 5)
##   index         Especie CT_cm Estado.maturacao
## 1     1 Pagellus acarne  10.0          Juvenil
## 2     2 Pagellus acarne  10.5          Juvenil
## 3     3 Pagellus acarne  11.0          Juvenil
## 4     4 Pagellus acarne  11.2          Juvenil
## 5     5 Pagellus acarne  11.5          Juvenil

Criar variavel binomial para o modelo logistico (imaturo = 0, maturo = 1)

dados.mat$maturacao.bin <- rep(0, nrow(dados.mat)) #criar a nova variavel (placeholder)
for (i in 1:nrow(dados.mat)){
  if (dados.mat$Estado.maturacao[i]=="Juvenil") dados.mat$maturacao.bin[i] <- 0 
  if (dados.mat$Estado.maturacao[i]=="Adulto") dados.mat$maturacao.bin[i] <- 1 
}

Ajustar um modelo logistico

Modelo logistico

modelo <- glm(formula = maturacao.bin ~ CT_cm, data = dados.mat, family = "binomial")
summary(modelo)
## 
## Call:
## glm(formula = maturacao.bin ~ CT_cm, family = "binomial", data = dados.mat)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -23.0967     4.6542  -4.963 6.96e-07 ***
## CT_cm         1.3147     0.2657   4.949 7.47e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 170.441  on 122  degrees of freedom
## Residual deviance:  68.515  on 121  degrees of freedom
## AIC: 72.515
## 
## Number of Fisher Scoring iterations: 7

Ver M50 e M95 no modelo ajustado

Extrair valores de tamanhos de 1a maturação (M50) e tamanho quando 95% dos individuos estão maturos

M50 <- MASS::dose.p(modelo, p=c(0.5)); M50
##              Dose        SE
## p = 0.5: 17.56771 0.2310515
M95 <- MASS::dose.p(modelo, p=c(0.95)); M95
##              Dose        SE
## p = 0.95: 19.8073 0.5233862

Previsões e gráficos das ogivas de maturação

Previões do modelo para a gama de dados observados

dados.novos <- data.frame(CT_cm = seq(min(dados.mat$CT_cm), max(dados.mat$CT_cm), 0.1))
previsoes <- predict(modelo, dados.novos, se.fit=T, type = "response")
dados.previsoes <- data.frame(CT_cm = dados.novos$CT_cm,
                              fit = previsoes$fit,
                              se.fit = previsoes$se.fit)

Gráfico simples dos dados observados e ogiva de maturação

plot(maturacao.bin ~ CT_cm, dados.mat,
     main = "Ogiva de maturação",
     xlab="Proporção de maturos",
     ylab="Comprimento total (cm)")
lines(dados.previsoes$CT_cm, dados.previsoes$fit)

Gráfico mais completo da ogiva e dados obervações

ggplot(dados.previsoes, aes(x = CT_cm, y = fit)) +
  geom_line(linewidth = 2, col = "royalblue") +
  geom_point(aes(M50[1], 0.5), col = 'black', size = 4)+
  scale_y_continuous(name = "Proporção de maturos",  limits = c(0,1), breaks = seq(0,1,0.2))+
  scale_x_continuous(name = "Comprimento total (cm)",limits =c(min(dados.previsoes$CT_cm), max(dados.previsoes$CT_cm)))+
  geom_point(data = dados.mat, aes(CT_cm, maturacao.bin, color=factor(maturacao.bin)), alpha = .4, size = 3, show.legend = FALSE) +
  theme_bw()