# 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," ",variavel,"################"))
    # 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,]
      temp[[paste(j,k)]]<-dados4
    } else {
      temp[[paste(j,k)]]<-dados3
    }
  }
}

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 | Dose)
               + (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:4], predicao)
colnames(dados6)<-c("Cultivar","Combinacao", "Dias", "PROD")

# rodando anova com covariável
dados6$Cultivar<-as.factor(dados6$Cultivar)
dados6$Combinacao<-as.factor(dados6$Combinacao)
aovAncova <- aov(PROD ~ Cultivar*Combinacao + Dias, data=dados6)
summary(aovAncova)
##                     Df   Sum Sq Mean Sq F value   Pr(>F)    
## Cultivar             1  1571148 1571148   6.616   0.0129 *  
## Combinacao           1  8416971 8416971  35.444 2.02e-07 ***
## Dias                 1    12664   12664   0.053   0.8182    
## Cultivar:Combinacao  1  1819911 1819911   7.664   0.0077 ** 
## Residuals           54 12823620  237474                     
## ---
## 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[,-3])
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       113.0775     a   a      a a           a
## 2 TMG_4377      3190.577        83.4019     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        78.3424     a   a      a a           a
## 2 N_combinado      3003.867       116.6399     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       110.7929     a   a      a a
## 1 TMG_4182.Combinado      3655.789       110.7929     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       197.1574     a   a      a a
## 4 TMG_4377.N_combinado      2648.733       124.6933     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       110.7929     a   a      a a
## 3 TMG_4182.N_combinado      3359.000       197.1574     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       110.7929     a   a      a a
## 4 TMG_4377.N_combinado      2648.733       124.6933     b   b      b b
##   scott_knott
## 2           a
## 4           b

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 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," ",variavel,"################"))
    # 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,]
      temp[[paste(j,k)]]<-dados4
    } else {
      temp[[paste(j,k)]]<-dados3
    }
  }
}

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 | Dose)
               + (1 | PAHP) + (1 | PPE) + (1 | PPC)
               + (1 | P10C) + (1 | THP) + (1 | T10C)
               + (1 | Stand) + (1 | Ciclo) + (1 | EMERG),
          data = dados5)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
predicao<-predict(modelo, dados5)

Anova com covariável

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

# rodando anova com covariável
dados6$Cultivar<-as.factor(dados6$Cultivar)
dados6$Combinacao<-as.factor(dados6$Combinacao)
aovAncova <- aov(Avariados ~ Cultivar*Combinacao + Dias, data=dados6)
summary(aovAncova)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Cultivar             1   1158  1158.3  11.679 0.00122 ** 
## Combinacao           1   2364  2364.5  23.841   1e-05 ***
## Dias                 1      1     1.4   0.014 0.90718    
## Cultivar:Combinacao  1    989   989.4   9.976 0.00262 ** 
## Residuals           53   5256    99.2                    
## ---
## 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[,-3])
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