Caderno de aulas
  • APRESENTAÇÃO
  • Aula 1
  • Aula 2
  • Aula 3
  • Aula 4
  • Aula 5
  • Aula 6
  • Aula de Mapa
library(MASS)
insects <- InsectSprays

m1 <- lm(count ~ spray, data = insects)

library(DHARMa)
plot(simulateResiduals(m1))

#não atendeu ao pressuposto de homogeneidade então fizemos transformação dos dados primeiro por raiz quadrada. 

m1 <- lm(sqrt(count) ~ spray, data = insects)
m1

Call:
lm(formula = sqrt(count) ~ spray, data = insects)

Coefficients:
(Intercept)       sprayB       sprayC       sprayD       sprayE       sprayF  
     3.7607       0.1160      -2.5158      -1.5963      -1.9512       0.2579  
library(DHARMa)
plot(simulateResiduals(m1))

#depois fizemos transformação por log - +1 para evitar log 0
m1_log <- lm(log(count + 1) ~ spray, data = insects)
m1_log

Call:
lm(formula = log(count + 1) ~ spray, data = insects)

Coefficients:
(Intercept)       sprayB       sprayC       sprayD       sprayE       sprayF  
     2.6967       0.0598      -1.7441      -0.9834      -1.2705       0.1189  
library(DHARMa)
plot(simulateResiduals(m1_log))

boxcox(lm(insects$count + 0.1 ~ 1))

b<- boxcox(lm(insects$count + 0.1 ~ 1))
lambda <- b$x[which.max(b$y)]
lambda
[1] 0.4242424
library(tidyverse)
insects <- insects |>
  mutate(count2 = count^lambda-1/lambda) |>
  mutate(count3 = sqrt(count))

insects$count2<- (insects$count ^ lambda - 1) / lambda
hist(insects$count2)

Aqui carregamos o conjunto de dados InsectSprays e ajustamos um modelo linear para verificar o efeito do tipo de inseticida sobre a contagem de insetos mortos. Em seguida, utilizamos o pacote DHARMa para avaliar os resíduos do modelo e verificamos possíveis desvios da normalidade. Depois aplicamos a transformação de Box-Cox, para melhorar a adequação do modelo e identificando o melhor valor de lambda. Também testamos a transformação por raiz quadrada para comparação.
library(gsheet)

estande <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=401662555#gid=401662555")

library(ggplot2)
 estande |> ggplot(aes(trat, nplants)) + geom_point() + geom_smooth(method = "lm", se = FALSE) +
   facet_wrap( ~ exp) +
   theme_minimal()

Aqui importamos e fizemos a visualização dos dados relacionados à emergência de plantas em função da porcentagem de inóculo nas sementes. utilizamos o ggplot2 para criar gráficos de dispersão entre o tratamento (trat) e o número de plantas emergidas (nplants), separando os gráficos por experimento (exp). Também adicionamos uma linha de tendência linear para facilitar a visualização do comportamento dos dados em cada experimento.
exp1 <- estande |>
  filter(exp == 1)
m_exp1 <- lm(nplants ~ trat, data = exp1)
summary(m_exp1)

Call:
lm(formula = nplants ~ trat, data = exp1)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.500  -6.532   1.758   8.573  27.226 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  52.5000     4.2044  12.487 1.84e-11 ***
trat         -0.2419     0.1859  -1.301    0.207    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15 on 22 degrees of freedom
Multiple R-squared:  0.07148,   Adjusted R-squared:  0.02928 
F-statistic: 1.694 on 1 and 22 DF,  p-value: 0.2066
exp2 <- estande |>
  filter(exp == 2)
m_exp2 <- lm(nplants ~ trat, data = exp2)
summary(m_exp2)

Call:
lm(formula = nplants ~ trat, data = exp2)

Residuals:
     Min       1Q   Median       3Q      Max 
-25.7816  -7.7150   0.5653   8.1929  19.2184 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  60.9857     3.6304  16.798 4.93e-14 ***
trat         -0.7007     0.1605  -4.365 0.000247 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.95 on 22 degrees of freedom
Multiple R-squared:  0.4641,    Adjusted R-squared:  0.4398 
F-statistic: 19.05 on 1 and 22 DF,  p-value: 0.0002473
exp3 <- estande |>
  filter(exp == 3)
m_exp3 <- lm(nplants ~ trat, data = exp3)
summary(m_exp3)

Call:
lm(formula = nplants ~ trat, data = exp3)

Residuals:
     Min       1Q   Median       3Q      Max 
-26.5887  -3.9597   0.7177   5.5806  19.8952 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  95.7500     2.9529  32.425  < 2e-16 ***
trat         -0.7634     0.1306  -5.847 6.97e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.53 on 22 degrees of freedom
Multiple R-squared:  0.6085,    Adjusted R-squared:  0.5907 
F-statistic: 34.19 on 1 and 22 DF,  p-value: 6.968e-06
Aqui realizamos análises para cada um dos experimentos da planilha importada anteriormente. Filtramos os dados por experimento e ajustamos modelos lineares para avaliar a relação entre o percentual de inóculo (trat) e o número de plantas emergidas (nplants). Para cada experimento, usamos a função lm() e analisamos os resultados por meio do summary(), permitindo identificar a significância e o efeito do tratamento em cada caso específico.
library(lme4)
m_misto <- lmer(nplants ~ trat + (1 | exp/bloco), data = estande)
summary(m_misto)
Linear mixed model fit by REML ['lmerMod']
Formula: nplants ~ trat + (1 | exp/bloco)
   Data: estande

REML criterion at convergence: 575.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.21697 -0.63351  0.04292  0.67094  1.92907 

Random effects:
 Groups    Name        Variance Std.Dev.
 bloco:exp (Intercept)  54.76    7.40   
 exp       (Intercept) 377.43   19.43   
 Residual              134.99   11.62   
Number of obs: 72, groups:  bloco:exp, 12; exp, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept) 69.74524   11.57191   6.027
trat        -0.56869    0.08314  -6.840

Correlation of Fixed Effects:
     (Intr)
trat -0.111
confint(m_misto)
                 2.5 %     97.5 %
.sig01       3.3332097 14.4218422
.sig02       7.2377419 47.8269818
.sigma       9.7314178 13.9359486
(Intercept) 43.4631239 96.0274587
trat        -0.7328972 -0.4044812
car::Anova(m_misto)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: nplants
      Chisq Df Pr(>Chisq)    
trat 46.788  1  7.909e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
estande |>
  ggplot(aes(trat, nplants, color = factor(exp))) +
  geom_point()+
  #geom_smooth(method = "lm", se = FALSE) +
  geom_abline(intercept = 69.74, 
              slope = -0.568, linewidth = 2) +
  geom_abline(intercept = 43,
              slope = -0.73, linetype ="dashed") +
    geom_abline(intercept = 96,
              slope = -0.40, linetype ="dashed")

Aqui realizamos uma análise com modelo misto para considerar variações entre experimentos e blocos. Fizemos isso ajustando um modelo misto linear com o pacote lme4, considerando o efeito fixo do tratamento (trat) e os efeitos aleatórios dos blocos dentro dos experimentos ((1 | exp/bloco)). Avaliamos o modelo com summary(), obtivemos os intervalos de confiança dos parâmetros com confint() e testamos a significância dos efeitos com a função Anova() do pacote car. Por fim, criamos um gráfico com os pontos dos dados coloridos por experimento e adicionamos retas com diferentes interceptos e inclinações para representar tendências gerais e específicas, comparando possíveis ajustes.
library(gsheet)
fungi <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1bq2N19DcZdtax2fQW9OHSGMR0X2__Z9T/edit?gid=465348652#gid=465348652")

fungi |>
  group_by(code, dose)  |>
  summarise(germination = mean(germination)) |> 
             ggplot(aes(dose, germination)) + geom_point() +
  geom_line() +
  facet_wrap(~ code)

FGT43 <- fungi |>
  group_by(code, dose) |>
  summarise(germination = mean(germination)) |>
  filter (code == "FGT43")

library(drc)
m43 <- drm(germination ~ dose, data = FGT43,
           fct = LL.3())
summary(m43)

Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms)

Parameter estimates:

               Estimate Std. Error t-value   p-value    
b:(Intercept)  1.219692   0.175081  6.9664  0.006069 ** 
d:(Intercept) 48.486911   1.456007 33.3013 5.952e-05 ***
e:(Intercept)  0.495895   0.060851  8.1494  0.003864 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 1.636105 (3 degrees of freedom)
AIC(m43)
[1] 26.7762
plot(m43)

ED(m43, 50)

Estimated effective doses

       Estimate Std. Error
e:1:50 0.495895   0.060851
library(ec50estimator)
df_ec50 = estimate_EC50(germination ~ dose,
                        data = fungi,
                        isolate_col = "code",
                        strata_col = "state",
                        interval = "delta",
                        fct = drc :: LL.3())
df_ec50 |>
  ggplot(aes(reorder(ID, Estimate), Estimate)) +
  geom_point()+
  coord_flip()

df_ec50 |>
  ggplot(aes(x = Estimate)) +
  geom_histogram(bins = 5, color = "white")

Aqui importamos os dados de germinação de esporos de diferentes isolados sob diferentes doses de tratamento, depois calculamos a média de germinação por isolado e dose e visualizamos os resultados com gráficos de dispersão e linhas. Selecionamos o isolado FGT43 e ajustamos um modelo não linear com o pacote drc, estimando a EC50 e visualizando o ajuste. Em seguida, utilizamos o pacote ec50estimator para estimar a EC50 de todos os isolados ao mesmo tempo, por estado. Por fim, apresentamos os resultados em gráfico de pontos ordenados e histograma das EC50s para avaliar a variabilidade entre os isolados.