Análise Ricardo

# Definindo a trilha de dados
setwd("C:/Users/Windows/Experimental Analytics Corporation/Experimental Analytics Corporation - Empresa/Customers/Ricardo_Miranda")

# Leitura dos dados
dados<-read.table("dados.csv", sep=",", h=T)
head(dados)

##   Talhao Cultivar  Combinacao Dias   Dose PAHP   PPE     PPC  P10C   THP  T10C
## 1     15 TMG 4377 N_combinado   16 2.5088 41.0  3.20  883.80  30.2 27.50 24.80
## 2     18 TMG 4377 N_combinado   14 2.5088 59.6 14.20  981.60  65.8 26.29 24.36
## 3     24 TMG 4377 N_combinado   22 2.5166 63.2 30.60 1446.80 150.0 26.81 23.87
## 4     27 TMG 4377   Combinado   15 1.5010 49.2  8.61  791.54  31.8 26.40 24.76
## 5     28 TMG 4377   Combinado   14 0.9882 49.2  8.61  791.34  31.8 26.27 24.75
## 6     33 TMG 4377   Combinado   11 0.9805 55.0  6.00 1517.40 259.8 26.70 23.86
##   PROD Avariados Stand Ciclo EMERG
## 1 3396      1.53 11.80    98     6
## 2 3648      1.67 11.90   104     6
## 3  907     55.91 11.23   119     6
## 4 3775      1.18 11.60   102     7
## 5 3800      1.39 11.90   101     7
## 6 2115     38.35  8.00   118     8

Análise para produtividade

# Identificação e remoção de outlier
library(ggplot2)

variavel<-"PROD"
dados1<-dados[complete.cases(dados$PROD),]

p <- ggplot(data = dados1, aes(x = PROD)) + geom_histogram() + facet_grid(Cultivar~Combinacao) +
  labs(y="Distribuição da Produtividade (kg/ha)", title='Distribuição da Produtividade') + theme(plot.title = element_text(hjust=0.5, size=20, face='bold')) 
p

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Removendo outliers para produtividade
temp<-list()
for (j in unique(dados1$Cultivar))
{
  dados2<-dados1[dados1$Cultivar==j,]
  
  for(k in unique(dados2$Combinacao))
  {
    dados3<-dados2[dados2$Combinacao==k,]
    
    print(paste("################",j," ",k," ################"))
    
    # via boxplot
    out <- boxplot.stats(dados3[,"PROD"])$out
    out_ind <- which(dados3[,"PROD"] %in% c(out))
    
    if(length(out_ind)>0)
    {
      row.names(dados3)<-seq(1:nrow(dados3))
      dados4<-dados3[-out_ind,]
      print(dados3[out_ind,])
      temp[[paste(j,k)]]<-dados4
    } else {
      temp[[paste(j,k)]]<-dados3
    }
  }
}

## [1] "################ TMG 4377   N_combinado  ################"
## [1] "################ TMG 4377   Combinado  ################"
##    Talhao Cultivar Combinacao Dias   Dose PAHP  PPE    PPC  P10C   THP  T10C
## 3      33 TMG 4377  Combinado   11 0.9805 55.0  6.0 1517.4 259.8 26.70 23.86
## 4      36 TMG 4377  Combinado    7 0.5928  9.6  4.4 1716.4 166.4 27.22 23.87
## 17      6 TMG 4377  Combinado    5 0.5919 12.4 33.6 1516.2 188.8 25.65 24.04
## 22     21 TMG 4377  Combinado   24 4.9710 60.8  9.8 1407.6 182.0 26.34 23.90
##    PROD Avariados Stand Ciclo EMERG
## 3  2115     38.35  8.00   118     8
## 4  1945     43.53 11.23   119     6
## 17 1254     58.11 11.30   121     7
## 22 1494     51.64 11.30   116     7
## [1] "################ TMG 4182   N_combinado  ################"
##   Talhao Cultivar  Combinacao Dias   Dose PAHP  PPE    PPC  P10C   THP  T10C
## 1     16 TMG 4182 N_combinado   12 2.5088 32.2  8.8 1490.2 296.4 27.26 24.16
## 5     95 TMG 4182 N_combinado   18 2.5070 24.4 65.6 1455.2 434.8 26.66 24.05
##   PROD Avariados Stand Ciclo EMERG
## 1 1326     57.86  11.3   123     5
## 5 3900      2.64  12.0   119     6
## [1] "################ TMG 4182   Combinado  ################"
##   Talhao Cultivar Combinacao Dias   Dose PAHP PPE   PPC P10C   THP  T10C PROD
## 3     51 TMG 4182  Combinado   16 1.5042   56 1.2 993.6 62.4 26.71 24.43 2096
##   Avariados Stand Ciclo EMERG
## 3      1.08  10.4   114     5

dados5<-data.frame(do.call(rbind,temp))

p <- ggplot(data = dados5, aes(x = PROD)) + geom_histogram() + facet_grid(Cultivar~Combinacao) +
  labs(y="Distribuição da Produtividade (kg/ha)", title='Distribuição da Produtividade') + theme(plot.title = element_text(hjust=0.5, size=20, face='bold')) 
p

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Montando o modelo misto para retirar os efeitos de pluviosidade, umidade e dose

library(lme4)

## Carregando pacotes exigidos: Matrix

modelo <- lmer(PROD ~ Cultivar + Combinacao 
               + Cultivar:Combinacao
               + (1 | PAHP) + (1 | PPE) + (1 | PPC)
               + (1 | P10C) + (1 | THP) + (1 | T10C)
               + (1 | Stand) + (1 | Ciclo) + (1 | EMERG),
          data = dados5)

## boundary (singular) fit: see help('isSingular')

predicao<-predict(modelo, dados5)

Anova com covariável

# Criando novo arquivo para rodar anova com covariável
dados6<-cbind(dados5[2:5], predicao)
colnames(dados6)<-c("Cultivar","Combinacao", "Dias", "Dose", "PROD")
dados6$Dose<-round(dados6$Dose,1)

# rodando anova com covariável
dados6$Cultivar<-as.factor(dados6$Cultivar)
dados6$Combinacao<-as.factor(dados6$Combinacao)
dados6$Dose<-as.factor(dados6$Dose)
aovAncova <- aov(PROD ~ Cultivar*Combinacao*Dose + Dias, data=dados6)
summary(aovAncova)

##                     Df  Sum Sq Mean Sq F value   Pr(>F)    
## Cultivar             1 1571148 1571148   8.070  0.00719 ** 
## Combinacao           1 8416971 8416971  43.235 9.29e-08 ***
## Dose                 7 2099056  299865   1.540  0.18345    
## Dias                 1  330175  330175   1.696  0.20065    
## Cultivar:Combinacao  1 1342845 1342845   6.898  0.01237 *  
## Cultivar:Dose        6 2424614  404102   2.076  0.07907 .  
## Combinacao:Dose      3  605713  201904   1.037  0.38715    
## Residuals           38 7397767  194678                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova do fatorial duplo

library(easyanova)
dados7<-data.frame(dados6[,c(1,2,5)])
dados7$Cultivar<-gsub( " ", "_", dados7$Cultivar)

resultado<-ea2(data = dados7, design = 1, alpha = 0.05)

print(resultado$`Adjusted means (factor 1)`)

##   factor_1 adjusted.mean standard.error tukey snk duncan t scott_knott
## 1 TMG_4182      3507.395       111.0493     a   a      a a           a
## 2 TMG_4377      3190.577        81.9059     b   b      b b           b

print(resultado$`Adjusted means (factor 2)`)

##      factor_2 adjusted.mean standard.error tukey snk duncan t scott_knott
## 1   Combinado      3694.105        76.9372     a   a      a a           a
## 2 N_combinado      3003.867       114.5478     b   b      b b           b

print(resultado$`Adjusted means (factor 1 in levels of factor 2)`)

## $`factor_1 in  Combinado`
##            treatment adjusted.mean standard.error tukey snk duncan t
## 2 TMG_4377.Combinado      3732.421       108.8057     a   a      a a
## 1 TMG_4182.Combinado      3655.789       108.8057     a   a      a a
##   scott_knott
## 2           a
## 1           a
## 
## $`factor_1 in  N_combinado`
##              treatment adjusted.mean standard.error tukey snk duncan t
## 3 TMG_4182.N_combinado      3359.000       193.6211     a   a      a a
## 4 TMG_4377.N_combinado      2648.733       122.4568     b   b      b b
##   scott_knott
## 3           a
## 4           b

print(resultado$`Adjusted means (factor 2 in levels of factor 1)`)

## $`factor_2 in  TMG_4182`
##              treatment adjusted.mean standard.error tukey snk duncan t
## 1   TMG_4182.Combinado      3655.789       108.8057     a   a      a a
## 3 TMG_4182.N_combinado      3359.000       193.6211     a   a      a a
##   scott_knott
## 1           a
## 3           a
## 
## $`factor_2 in  TMG_4377`
##              treatment adjusted.mean standard.error tukey snk duncan t
## 2   TMG_4377.Combinado      3732.421       108.8057     a   a      a a
## 4 TMG_4377.N_combinado      2648.733       122.4568     b   b      b b
##   scott_knott
## 2           a
## 4           b

print(resultado$`Residual analysis`$`residual analysis`)

##                                     values
## p.value Shapiro-Wilk test           0.0010
## p.value Bartlett test (factor_1)    0.0001
## p.value Bartlett test (factor_2)    0.0000
## p.value Bartlett test (treatments)  0.0000
## coefficient of variation (%)       13.9700
## first value most discrepant         3.0000
## second value most discrepant       15.0000
## third value most discrepant        13.0000

Análise para avariados

# Identificação e remoção de outlier
library(ggplot2)

variavel<-"Avariados"
dados1<-dados[complete.cases(dados$Avariados),]

p <- ggplot(data = dados1, aes(x = Avariados)) + geom_histogram() + facet_grid(Cultivar~Combinacao) +
  labs(y="Distribuição de Avariados)", title='Distribuição de Avariados') + theme(plot.title = element_text(hjust=0.5, size=20, face='bold')) 
p

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Removendo outliers para avariados
temp<-list()
for (j in unique(dados1$Cultivar))
{
  dados2<-dados1[dados1$Cultivar==j,]
  
  for(k in unique(dados2$Combinacao))
  {
    dados3<-dados2[dados2$Combinacao==k,]
    
    print(paste("################",j," ",k," ################"))
    # via boxplot
    out <- boxplot.stats(dados3[,"Avariados"])$out
    out_ind <- which(dados3[,"Avariados"] %in% c(out))
    
    if(length(out_ind)>0)
    {
      row.names(dados3)<-seq(1:nrow(dados3))
      dados4<-dados3[-out_ind,]
      print(dados3[out_ind,])
      temp[[paste(j,k)]]<-dados4
    } else {
      temp[[paste(j,k)]]<-dados3
    }
  }
}

## [1] "################ TMG 4377   N_combinado  ################"
## [1] "################ TMG 4377   Combinado  ################"
##    Talhao Cultivar Combinacao Dias   Dose PAHP  PPE    PPC  P10C   THP  T10C
## 3      33 TMG 4377  Combinado   11 0.9805 55.0  6.0 1517.4 259.8 26.70 23.86
## 4      36 TMG 4377  Combinado    7 0.5928  9.6  4.4 1716.4 166.4 27.22 23.87
## 17      6 TMG 4377  Combinado    5 0.5919 12.4 33.6 1516.2 188.8 25.65 24.04
## 22     21 TMG 4377  Combinado   24 4.9710 60.8  9.8 1407.6 182.0 26.34 23.90
##    PROD Avariados Stand Ciclo EMERG
## 3  2115     38.35  8.00   118     8
## 4  1945     43.53 11.23   119     6
## 17 1254     58.11 11.30   121     7
## 22 1494     51.64 11.30   116     7
## [1] "################ TMG 4182   N_combinado  ################"
##   Talhao Cultivar  Combinacao Dias   Dose PAHP PPE    PPC  P10C   THP  T10C
## 1     16 TMG 4182 N_combinado   12 2.5088 32.2 8.8 1490.2 296.4 27.26 24.16
##   PROD Avariados Stand Ciclo EMERG
## 1 1326     57.86  11.3   123     5
## [1] "################ TMG 4182   Combinado  ################"
##    Talhao Cultivar Combinacao Dias   Dose PAHP  PPE    PPC  P10C   THP  T10C
## 5      77 TMG 4182  Combinado    8 0.9924 49.2 12.4 1301.6 391.8 26.50 23.76
## 6      78 TMG 4182  Combinado    6 0.5922 21.4 27.4 1397.2 427.6 26.00 23.80
## 20      3 TMG 4182  Combinado    5 1.0295 27.0  8.6 1141.8  61.0 26.78 24.10
##    PROD Avariados Stand Ciclo EMERG
## 5  3886      3.96  12.5   117     5
## 6  3087     10.88  11.3   118     5
## 20 3969      0.46  10.9   106     8

dados5<-data.frame(do.call(rbind,temp))

p <- ggplot(data = dados5, aes(x = Avariados)) + geom_histogram() + facet_grid(Cultivar~Combinacao) +
  labs(y="Distribuição de Avariados", title='Distribuição de Avariados') + theme(plot.title = element_text(hjust=0.5, size=20, face='bold')) 
p

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Montando o modelo misto para retirar os efeitos de pluviosidade, umidade e dose

library(lme4)

modelo <- lmer(Avariados ~ Cultivar + Combinacao 
               + Cultivar:Combinacao
               + (1 | PAHP) + (1 | PPE) + (1 | PPC)
               + (1 | P10C) + (1 | THP) + (1 | T10C)
               + (1 | Stand) + (1 | Ciclo) + (1 | EMERG),
          data = dados5)

## boundary (singular) fit: see help('isSingular')

predicao<-predict(modelo, dados5)

Anova com covariável

# Criando novo arquivo para rodar anova com covariável
dados6<-cbind(dados5[2:5], predicao)
colnames(dados6)<-c("Cultivar","Combinacao", "Dias","Dose", "Avariados")
dados6$Dose<-round(dados6$Dose,1)

# rodando anova com covariável
dados6$Cultivar<-as.factor(dados6$Cultivar)
dados6$Combinacao<-as.factor(dados6$Combinacao)
dados6$Dose<-as.factor(dados6$Dose)
aovAncova <- aov(Avariados ~ Cultivar*Combinacao*Dose + Dias, data=dados6)
summary(aovAncova)

##                     Df Sum Sq Mean Sq F value   Pr(>F)    
## Cultivar             1 1158.3  1158.3  14.737 0.000467 ***
## Combinacao           1 2364.5  2364.5  30.084 3.13e-06 ***
## Dose                 7 1006.0   143.7   1.829 0.110777    
## Dias                 1  201.9   201.9   2.568 0.117519    
## Cultivar:Combinacao  1  835.0   835.0  10.624 0.002398 ** 
## Cultivar:Dose        6 1033.8   172.3   2.192 0.065708 .  
## Combinacao:Dose      3  262.5    87.5   1.113 0.356121    
## Residuals           37 2908.0    78.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova do fatorial duplo

# Rodando sem a covariável
library(easyanova)
dados7<-data.frame(dados6[,c(1,2,5)])
dados7$Cultivar<-gsub( " ", "_", dados7$Cultivar)

resultado<-ea2(data = dados7, design = 1, alpha = 0.05)

print(resultado$`Adjusted means (factor 1)`)

##   factor_1 adjusted.mean standard.error tukey snk duncan t scott_knott
## 1 TMG_4377       12.0897         1.7039     a   a      a a           a
## 2 TMG_4182        2.2827         2.2154     b   b      b b           b

print(resultado$`Adjusted means (factor 2)`)

##      factor_2 adjusted.mean standard.error tukey snk duncan t scott_knott
## 1 N_combinado       12.6986         2.2581     a   a      a a           a
## 2   Combinado        1.6739         1.6469     b   b      b b           b

print(resultado$`Adjusted means (factor 1 in levels of factor 2)`)

## $`factor_1 in  Combinado`
##            treatment adjusted.mean standard.error tukey snk duncan t
## 2 TMG_4377.Combinado        2.1195         2.2635     a   a      a a
## 1 TMG_4182.Combinado        1.2282         2.3929     a   a      a a
##   scott_knott
## 2           a
## 1           a
## 
## $`factor_1 in  N_combinado`
##              treatment adjusted.mean standard.error tukey snk duncan t
## 4 TMG_4377.N_combinado       22.0600         2.5475     a   a      a a
## 3 TMG_4182.N_combinado        3.3371         3.7291     b   b      b b
##   scott_knott
## 4           a
## 3           b

print(resultado$`Adjusted means (factor 2 in levels of factor 1)`)

## $`factor_2 in  TMG_4182`
##              treatment adjusted.mean standard.error tukey snk duncan t
## 3 TMG_4182.N_combinado        3.3371         3.7291     a   a      a a
## 1   TMG_4182.Combinado        1.2282         2.3929     a   a      a a
##   scott_knott
## 3           a
## 1           a
## 
## $`factor_2 in  TMG_4377`
##              treatment adjusted.mean standard.error tukey snk duncan t
## 4 TMG_4377.N_combinado       22.0600         2.5475     a   a      a a
## 2   TMG_4377.Combinado        2.1195         2.2635     b   b      b b
##   scott_knott
## 4           a
## 2           b

print(resultado$`Residual analysis`$`residual analysis`)

##                                    values
## p.value Shapiro-Wilk test            0.00
## p.value Bartlett test (factor_1)     0.00
## p.value Bartlett test (factor_2)     0.00
## p.value Bartlett test (treatments)   0.00
## coefficient of variation (%)       137.75
## first value most discrepant          3.00
## second value most discrepant        14.00
## third value most discrepant          1.00

Análise Ricardo

Leonardo Peixoto

2022-06-29

Análise para produtividade

Montando o modelo misto para retirar os efeitos de pluviosidade, umidade e dose

Anova com covariável

Anova do fatorial duplo

Análise para avariados

Montando o modelo misto para retirar os efeitos de pluviosidade, umidade e dose

Anova com covariável

Anova do fatorial duplo