#######################
##### Cid Edson #######
#######################

Resultados previos

Tratamento

vs <-
  read_excel("completos.xlsx", sheet = 1) # maio a outubro de 2017/2018
vs <- as.data.frame(vs)
vs$loc <- as.factor(vs$loc)
vs$trat <- as.factor(vs$trat)
vs$m <- as.factor(vs$m)
vs$ano <- as.factor(vs$ano)
vs <- na.exclude(vs)

vr <-
  subset(vs,
         m %in% c("maio", "junho", "julho", "agosto", "setembro", "outubro"))

levels(vr$trat) <- c("Controle", "Tratamento")
levels(vr$loc) <- c("Camaca", "Cepec", "Itubera")


mi <- lapply(split(vr, f = vr$loc), aov, formula = fsa ~ trat)
lapply(mi, summary)

## $Camaca
##              Df Sum Sq Mean Sq F value Pr(>F)
## trat          1     41   40.84   2.533  0.112
## Residuals   958  15447   16.12               
## 
## $Cepec
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## trat          1    107  106.67   13.58 0.000242 ***
## Residuals   958   7526    7.86                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Itubera
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## trat          1   2637  2636.8    56.3 1.42e-13 ***
## Residuals   958  44869    46.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mi$Camaca$residuals %>% shapiro.test

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.68946, p-value < 2.2e-16

mi$Cepec$residuals %>% shapiro.test

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.80016, p-value < 2.2e-16

mi$Itubera$residuals %>% shapiro.test

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.86817, p-value < 2.2e-16

caw <- subset(vr, loc %in% "Camaca")
cew <- subset(vr, loc %in% "Cepec")
itw <- subset(vr, loc %in% "Itubera")


with(caw,
     kruskal(
       fsa,
       trat,
       p.adj = "bon",
       group = T,
       main = "Camaca",
       console = T
     ))

## 
## Study: Camaca
## Kruskal-Wallis test's
## Ties or no Ties
## 
## Critical Value: 6.279454
## Degrees of freedom: 1
## Pvalue Chisq  : 0.01221457 
## 
## trat,  means of the ranks
## 
##                 fsa   r
## Controle   458.4781 480
## Tratamento 502.5219 480
## 
## Post Hoc Analysis
## 
## P value adjustment method: bonferroni
## t-Student: 1.962443
## Alpha    : 0.05
## Minimum Significant Difference: 34.39699 
## 
## Treatments with the same letter are not significantly different.
## 
##                 fsa groups
## Tratamento 502.5219      a
## Controle   458.4781      b

with(cew,
     kruskal(
       fsa,
       trat,
       p.adj = "bon",
       group = T,
       main = "Cepec",
       console = T
     ))

## 
## Study: Cepec
## Kruskal-Wallis test's
## Ties or no Ties
## 
## Critical Value: 17.52555
## Degrees of freedom: 1
## Pvalue Chisq  : 2.834726e-05 
## 
## trat,  means of the ranks
## 
##                 fsa   r
## Controle   443.9698 480
## Tratamento 517.0302 480
## 
## Post Hoc Analysis
## 
## P value adjustment method: bonferroni
## t-Student: 1.962443
## Alpha    : 0.05
## Minimum Significant Difference: 33.95197 
## 
## Treatments with the same letter are not significantly different.
## 
##                 fsa groups
## Tratamento 517.0302      a
## Controle   443.9698      b

with(itw,
     kruskal(
       fsa,
       trat,
       p.adj = "bon",
       group = T,
       main = "Itubera",
       console = T
     ))

## 
## Study: Itubera
## Kruskal-Wallis test's
## Ties or no Ties
## 
## Critical Value: 86.35989
## Degrees of freedom: 1
## Pvalue Chisq  : 0 
## 
## trat,  means of the ranks
## 
##                 fsa   r
## Controle   397.5833 480
## Tratamento 563.4167 480
## 
## Post Hoc Analysis
## 
## P value adjustment method: bonferroni
## t-Student: 1.962443
## Alpha    : 0.05
## Minimum Significant Difference: 33.42318 
## 
## Treatments with the same letter are not significantly different.
## 
##                 fsa groups
## Tratamento 563.4167      a
## Controle   397.5833      b

# res <- vr[, c(1,2,8,4,5,6,7)]
#
# rez = ea2(res,
#              design = 9,
#              list = TRUE,
#              plot = 3)
#
# rez


# Frutos sadios
bwplot(
  fsa ~ trat | loc * ano,
  vr,
  ylab = "Prod",
  aspect = "fill",
  main = "Frutos sadios",
  do.out = FALSE,
  par.settings = bw_ps
)

bwplot(
  fvb ~ trat | loc * ano,
  vr,
  ylab = "Prod",
  aspect = "fill",
  main = "Frutos com Vassoura",
  do.out = FALSE,
  par.settings = bw_ps
)

fsa0 <- aggregate(fsa ~ trat + loc + m +  ano, vr, "sum")
fvb0 <- aggregate(fvb ~ trat + loc + m +  ano, vr, "sum")
fpp0 <- aggregate(fpp ~ trat + loc + m + ano, vr, "sum")

exper <- data.frame(fsa0, fvb0$fvb, fpp0$fpp)

names(exper) <-
  c("Tratamento",
    "Local",
    "Mes",
    "Ano",
    "Sadios",
    "V. Bruxa",
    "Podridao P.")

exper$Mes <- factor(exper$Mes, c("maio", "junho", "julho", "agosto", "setembro", "outubro")) 

exper %>% kable

Tratamento	Local	Mes	Ano	Sadios	V. Bruxa	Podridao P.
Controle	Camaca	agosto	2017	94	7	4
Tratamento	Camaca	agosto	2017	54	3	4
Controle	Cepec	agosto	2017	137	6	3
Tratamento	Cepec	agosto	2017	147	6	1
Controle	Itubera	agosto	2017	187	52	11
Tratamento	Itubera	agosto	2017	320	92	9
Controle	Camaca	julho	2017	64	12	2
Tratamento	Camaca	julho	2017	64	4	0
Controle	Cepec	julho	2017	115	5	6
Tratamento	Cepec	julho	2017	132	4	6
Controle	Itubera	julho	2017	214	31	12
Tratamento	Itubera	julho	2017	338	92	0
Controle	Camaca	junho	2017	57	8	3
Tratamento	Camaca	junho	2017	58	3	0
Controle	Cepec	junho	2017	117	145	54
Tratamento	Cepec	junho	2017	186	180	35
Controle	Itubera	junho	2017	297	64	12
Tratamento	Itubera	junho	2017	461	48	49
Controle	Camaca	maio	2017	230	15	2
Tratamento	Camaca	maio	2017	227	39	3
Controle	Cepec	maio	2017	69	122	15
Tratamento	Cepec	maio	2017	227	77	32
Controle	Itubera	maio	2017	338	44	74
Tratamento	Itubera	maio	2017	456	60	23
Controle	Camaca	outubro	2017	111	59	3
Tratamento	Camaca	outubro	2017	73	48	1
Controle	Cepec	outubro	2017	154	42	2
Tratamento	Cepec	outubro	2017	155	40	1
Controle	Itubera	outubro	2017	175	15	6
Tratamento	Itubera	outubro	2017	313	36	5
Controle	Camaca	setembro	2017	128	37	0
Tratamento	Camaca	setembro	2017	78	18	0
Controle	Cepec	setembro	2017	181	15	1
Tratamento	Cepec	setembro	2017	189	24	1
Controle	Itubera	setembro	2017	97	13	0
Tratamento	Itubera	setembro	2017	162	7	0
Controle	Camaca	agosto	2018	47	20	1
Tratamento	Camaca	agosto	2018	122	56	14
Controle	Cepec	agosto	2018	13	9	0
Tratamento	Cepec	agosto	2018	26	8	5
Controle	Itubera	agosto	2018	160	30	5
Tratamento	Itubera	agosto	2018	350	58	3
Controle	Camaca	julho	2018	76	25	25
Tratamento	Camaca	julho	2018	231	35	63
Controle	Cepec	julho	2018	17	8	1
Tratamento	Cepec	julho	2018	31	7	0
Controle	Itubera	julho	2018	238	91	16
Tratamento	Itubera	julho	2018	344	140	5
Controle	Camaca	junho	2018	187	56	1
Tratamento	Camaca	junho	2018	211	10	7
Controle	Cepec	junho	2018	21	9	0
Tratamento	Cepec	junho	2018	31	6	0
Controle	Itubera	junho	2018	350	64	3
Tratamento	Itubera	junho	2018	480	37	2
Controle	Camaca	maio	2018	238	31	0
Tratamento	Camaca	maio	2018	234	14	0
Controle	Cepec	maio	2018	21	4	0
Tratamento	Cepec	maio	2018	27	12	0
Controle	Itubera	maio	2018	374	117	6
Tratamento	Itubera	maio	2018	424	34	4
Controle	Camaca	outubro	2018	20	16	0
Tratamento	Camaca	outubro	2018	31	17	8
Controle	Cepec	outubro	2018	38	12	0
Tratamento	Cepec	outubro	2018	47	0	0
Controle	Itubera	outubro	2018	62	4	3
Tratamento	Itubera	outubro	2018	215	46	13
Controle	Camaca	setembro	2018	31	18	2
Tratamento	Camaca	setembro	2018	98	6	0
Controle	Cepec	setembro	2018	42	17	0
Tratamento	Cepec	setembro	2018	47	10	0
Controle	Itubera	setembro	2018	113	63	6
Tratamento	Itubera	setembro	2018	333	129	8

xtabs(fsa + fvb + fpp ~ loc + ano, vr) %>% kable

	2017	2018
Camaca	1513	1951
Cepec	2632	469
Itubera	4113	4330

exper1 <- aggregate(Sadios ~ Tratamento + Local + Ano, exper, "sum")
exper2 <- aggregate(`V. Bruxa` ~ Tratamento + Local + Ano, exper, "sum")
exper3 <- data.frame(exper1,exper2[4])
names(exper3) <- c("Tratamento","Bloco","Ano","Sadios","Doentes")
exper3 %>% kable

Tratamento	Bloco	Ano	Sadios	Doentes
Controle	Camaca	2017	684	138
Tratamento	Camaca	2017	554	115
Controle	Cepec	2017	773	335
Tratamento	Cepec	2017	1036	331
Controle	Itubera	2017	1308	219
Tratamento	Itubera	2017	2050	335
Controle	Camaca	2018	599	166
Tratamento	Camaca	2018	927	138
Controle	Cepec	2018	152	59
Tratamento	Cepec	2018	209	43
Controle	Itubera	2018	1297	369
Tratamento	Itubera	2018	2146	444

# Frutos sadios

barchart(
  Sadios ~ Tratamento | Local,
  groups = Ano,
  auto.key = list(space = "right", title = "Ano"),
  exper,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos sadios"
)

barchart(
  Sadios ~  Local | Ano,
  exper,
  groups = Tratamento,
  stack = F,
  aspect = "fill",
  auto.key = list(space = 'right', title = "Tratamento"),
  ylab = "Prod",
  main = "Frutos sadios"
)

barchart(
  `V. Bruxa` ~ Tratamento | Local,
  groups = Ano,
  auto.key = list(space = "right", title = "Ano"),
  exper,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos com Vassoura"
)

barchart(
  `V. Bruxa` ~  Local | Ano,
  exper,
  groups = Tratamento,
  stack = F,
  aspect = "fill",
  auto.key = list(space = 'right', title = "Tratamento"),
  ylab = "Prod",
  main = "Frutos com Vassoura"
)

barchart(
  `Podridao P.` ~ Tratamento | Local,
  groups = Ano,
  auto.key = list(space = "right", title = "Ano"),
  exper,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos com Podridao"
)

barchart(
  `Podridao P.` ~  Local | Ano,
  exper,
  groups = Tratamento,
  stack = F,
  aspect = "fill",
  auto.key = list(space = 'right', title = "Tratamento"),
  ylab = "Prod",
  main = "Frutos com Podridao"
)

# Dados

gf <- xtabs(fsa + fvb + fpp ~ loc, vr)
gf <- as.data.frame(gf)
gp <- round(prop.table(gf$Freq) * 100, 2)
names(gp) <- c("Camaca", "Cepec", "Itubera")
gw <- data.frame(gf[-1], gp)
names(gw) <- c("Total", "(%) Frutos")
pie(gp, col = rainbow(5))

gw %>% kable

	Total	(%) Frutos
Camaca	3464	23.08
Cepec	3101	20.66
Itubera	8443	56.26

tf <- xtabs(fsa + fvb + fpp ~ loc + ano, vr)
tf <- as.data.frame(tf)
tp1 <- round(prop.table(tf[1:3,]$Freq) * 100, 2)
twp <- data.frame(tf[1:3,][-2], tp1)
names(twp) <- c("Local", "Frutos 2017", "(%) Frutos")
names(tp1) <- c("Camaca", "Cepec", "Itubera")
pie(tp1, col = rainbow(5))

twp %>% kable  # total em 2017

Local	Frutos 2017	(%) Frutos
Camaca	1513	18.32
Cepec	2632	31.87
Itubera	4113	49.81

tf2 <- xtabs(fsa + fvb + fpp ~ loc + ano, vr)
tf2 <- as.data.frame(tf2)
tp12 <- round(prop.table(tf2[4:6,]$Freq) * 100, 2)
twp2 <- data.frame(tf2[4:6,][-2], tp12)
names(twp2) <- c("Local", "Frutos 2018", "(%) Frutos")
names(tp12) <- c("Camaca", "Cepec", "Itubera")
pie(tp12, col = rainbow(5))

twp2 %>% kable  # total em 2018

	Local	Frutos 2018	(%) Frutos
4	Camaca	1951	28.90
5	Cepec	469	6.95
6	Itubera	4330	64.15

ca <- subset(exper, Local %in% "Camaca")
ce <- subset(exper, Local %in% "Cepec")
it <- subset(exper, Local %in% "Itubera")

ca17 <- ca[1:12, 5:7]
ca18 <- ca[13:24, 5:7]
ce17 <- ce[1:12, 5:7]
ce18 <- ce[13:24, 5:7]
it17 <- it[1:12, 5:7]
it18 <- it[13:24, 5:7]

ca17 <- as.data.frame(ca17)
pca17 <- prop.table(ca17[1,] + ca17[2,]) * 100
pie(t(pca17)[1:3], col = rainbow(5), main = "Camaca 2017")

round(t(pca17), 2) %>% kable

	1
Sadios	89.16
V. Bruxa	6.02
Podridao P.	4.82

ca18 <- as.data.frame(ca18)
pca18 <- prop.table(ca18[1,] + ca18[2,]) * 100
pie(t(pca18)[1:3], col = rainbow(5), main = "Camaca 2018")

round(t(pca18), 2) %>% kable

	37
Sadios	65.00
V. Bruxa	29.23
Podridao P.	5.77

ce17 <- as.data.frame(ce17)
pce17 <- prop.table(ce17[1,] + ce17[2,]) * 100
pie(t(pce17)[1:3], col = rainbow(5), main = "Cepec 2017")

round(t(pce17), 2) %>% kable

	3
Sadios	94.67
V. Bruxa	4.00
Podridao P.	1.33

ce18 <- as.data.frame(ce18)
pce18 <- prop.table(ce18[1,] + ce18[2,]) * 100
pie(t(pce18)[1:3], col = rainbow(5), main = "Cepec 2018")

round(t(pce18), 2) %>% kable

	39
Sadios	63.93
V. Bruxa	27.87
Podridao P.	8.20

it17 <- as.data.frame(it17)
pit17 <- prop.table(it17[1,] + it17[2,]) * 100
pie(t(pit17)[1:3], col = rainbow(5), main = "Itubera 2017")

round(t(pit17), 2) %>% kable

	5
Sadios	75.56
V. Bruxa	21.46
Podridao P.	2.98

it18 <- as.data.frame(it18)
pit18 <- prop.table(it18[1,] + it18[2,]) * 100
pie(t(pit18)[1:3], col = rainbow(5), main = "Itubera 2018")

round(t(pit18), 2) %>% kable

	41
Sadios	84.16
V. Bruxa	14.52
Podridao P.	1.32

sf <- xtabs(fsa ~ loc, vr)
sf <- as.data.frame(sf)
sp <- round(prop.table(sf$Freq) * 100, 2)
sw <- data.frame(sf, sp)
uf <- xtabs(fvb + fpp ~ loc, vr)
uf <- as.data.frame(uf)
up <- round(prop.table(uf$Freq) * 100, 2)
uw <- data.frame(uf, up)
doe <- round(uf$Freq / sf$Freq * 100, 2)
taq <- data.frame(sf, uw$Freq, doe)
names(taq) <- c("Local", "Sadios", "Doentes", "Doentes (%)")
taq %>% kable

Local	Sadios	Doentes	Doentes (%)
Camaca	2764	700	25.33
Cepec	2170	931	42.90
Itubera	6801	1642	24.14

bwplot(
  Sadios ~ Tratamento | Mes*Ano,
  auto.key = list(space = "right", title = "Ano"),
  exper,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos sadios"
)

bwplot(
  `V. Bruxa` ~ Tratamento | Mes*Ano,
  auto.key = list(space = "right", title = "Ano"),
  exper,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos com Vassoura"
)

bwplot(
  `Podridao P.` ~ Tratamento | Mes*Ano,
  auto.key = list(space = "right", title = "Ano"),
  exper,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos com Podridao"
)

xyplot(Sadios ~ Tratamento, 
       group=Local,
       data=exper,
       type=c("p","a"),
       auto.key = list(space = 'right'),
       pch=19,
       jitter.x=T,
       ylab = "Frutos",
       xlab = "Tratamento",
       cex=1,
       main = "Frutos Sadios",
       scales = list(x = list(rot = 45)))

xyplot(Sadios ~ Tratamento|Ano, 
       group=Local,
       data=exper,
       type=c("p","a"),
       auto.key = list(space = 'right'),
       pch=19,
       jitter.x=T,
       ylab = "Frutos",
       xlab = "Tratamento",
       main = "Frutos Sadios",
       cex=1,
       scales = list(x = list(rot = 45)))

xyplot(Sadios ~ Tratamento|Mes*Ano, 
       group=Local,
       data=exper,
       type=c("p","a"),
       auto.key = list(space = 'right'),
       pch=19,
       jitter.x=T,
       ylab = "Frutos",
       xlab = "Tratamento",
       cex=1,
       main = "Frutos Sadios",
       scales = list(x = list(rot = 45)))

xyplot(`V. Bruxa` ~ Tratamento|Mes*Ano, 
       group=Local,
       data=exper,
       type=c("p","a"),
       auto.key = list(space = 'right'),
       pch=19,
       jitter.x=T,
       main = "Frutos com Vassoura",
       ylab = "Frutos",
       xlab = "Tratamento",
       cex=1,
       scales = list(x = list(rot = 45)))

itubera <- subset(exper, Local %in% "Itubera" & Ano %in% "2018" & Tratamento %in% "Tratamento" & Mes %in% "maio")
itubera2 <- subset(exper, Local %in% "Itubera" & Ano %in% "2018" & Tratamento %in% "Controle" & Mes %in% "maio")

xyplot(`Podridao P.` ~ Tratamento|Ano, 
       group=Local,
       data=exper,
       type=c("p","a"),
       auto.key = list(space = 'right'),
       pch=19,
       jitter.x=T,
       ylab = "Frutos",
       main = "Frutos com Podridao",
       xlab = "Tratamento",
       cex=1,
       scales = list(x = list(rot = 45)))

exper17 <- exper[1:36, ]
exper18 <- exper[37:72, ]

exper17s <- subset(exper17, Tratamento %in% "Controle")
exper17c <- subset(exper17, Tratamento %in% "Tratamento")

# cor(exper17s[, 5:7])
# cor(exper17c[, 5:7])

exper18s <- subset(exper18, Tratamento %in% "Controle")
exper18c <- subset(exper18, Tratamento %in% "Tratamento")

# cor(exper18s[, 5:7])
# cor(exper18c[, 5:7])

expers <- subset(exper, Tratamento %in% "Controle")
experc <- subset(exper, Tratamento %in% "Tratamento")

# cor(expers[, 5:7])
# cor(experc[, 5:7])


plot(
  bpca(exper[, 5:7]),
  type = "bp",
  var.f = .8,
  var.cex = 1.2,
  var.color = c("blue", "red", "green", "black"),
  obj.cex = 1,
  obj.names = F,
  obj.id = 12,
  var.id = 1
)

source("E:/Google/UESC/R/Scritps/Stephane/AA/ggbiplot2.R")

exper.pca = prcomp(exper[, 5:7], scale = T)

ggbiplot2(
  exper.pca,
  obs.scale = 1,
  var.scale = 1,
  ellipse = F,
  circle = T,
  grupos = exper$Local
)

# caw <- subset(vr, loc %in% "Camaca")
# cew <- subset(vr, loc %in% "Cepec" & fsa > 0)
# itw <- subset(vr, loc %in% "Itubera" & fsa > 0)
# 
# require(MASS)
# boxcox(fsa ~ trat, data=caw, plotit = F)
# 
# boxcox(fsa ~ trat, data=caw, lam=seq(-1, 1, 1/10))
# #caw$fsat <- (caw$fsa^(0.1) - 1)/0.1
# caw$fsat <- (rank(c(caw$fsa)))^1/2
# 
# caw.avt <- aov(log10(fsat) ~ trat, data=caw)
# summary(caw.avt)
# shapiro.test(caw.avt$residuals)
# shapiro.test(caw$fsat)
# caw
# 
# m0 <- aov(fsat ~ trat, caw)
# shapiro.test(m0$residuals)
# boxcox(m0, plotit = F)


# cap <- subset(exper, Local %in% "Camaca")
# cep <- subset(exper, Local %in% "Cepec")
# itp <- subset(exper, Local %in% "Itubera")
# 
# 
# as <- aov(Sadios ~ Tratamento*Ano, cap)
# summary(as)
# shapiro.test(as$residuals)
# qqnorm(as$residuals)
# qqline(as$residuals)
# 
# 
# 
# at <- aov(Sadios ~ Tratamento*Ano, cep)
# summary(at)
# shapiro.test(at$residuals)
# qqnorm(at$residuals)
# qqline(at$residuals)
# 
# 
# aq <- aov(Sadios ~ Tratamento*Ano, itp)
# summary(aq)
# shapiro.test(aq$residuals)
# qqnorm(aq$residuals)
# qqline(aq$residuals)
# 
#  wb <- aggregate(vp$fsa,
#                 by = list(l = vp$loc,
#                           t = vp$trat,
#                           m = vp$m,
#                           a = vp$ano),
#                 FUN = mean)
# 
# out<-with(wb,friedman(a,t,x,alpha=0.05, group=TRUE,console=TRUE,
#                          main="Tratamento"))
# plot(out,variation="IQR")
# boxplot(x ~ t*l,wb)



opts_chunk$set(
  results = "show",
  cache = FALSE,
  tidy = FALSE,
  fig.width = 7,
  fig.height = 6,
  fig.align = "center",
  dpi = 300,
  dev = "png",
  dev.args = list(png = list(family = "Ubuntu Light", type = "cairo"))
)
opts_chunk$set(fig.path = "fig/script01")
options(width = 90)
options(knitr.table.format = "pandoc")
options(knitr.kable.NA = '')

bw_ps <- list(
  plot.symbol = list(col = 1),
  box.dot = list(pch = "|"),
  box.umbrella = list(col = 1, lty = 1),
  box.rectangle = list(col = 1, fill = "seagreen")
)

Resultados previos

Tratamento

#exper3 <- read.table("niella.txt", header = T)
exper3$Ano <- as.factor(exper3$Ano)

exper3

##    Tratamento   Bloco  Ano Sadios Doentes
## 1    Controle  Camaca 2017    684     138
## 2  Tratamento  Camaca 2017    554     115
## 3    Controle   Cepec 2017    773     335
## 4  Tratamento   Cepec 2017   1036     331
## 5    Controle Itubera 2017   1308     219
## 6  Tratamento Itubera 2017   2050     335
## 7    Controle  Camaca 2018    599     166
## 8  Tratamento  Camaca 2018    927     138
## 9    Controle   Cepec 2018    152      59
## 10 Tratamento   Cepec 2018    209      43
## 11   Controle Itubera 2018   1297     369
## 12 Tratamento Itubera 2018   2146     444

ftable(Bloco~Ano+Tratamento, data=exper3)

##                 Bloco Camaca Cepec Itubera
## Ano  Tratamento                           
## 2017 Controle              1     1       1
##      Tratamento            1     1       1
## 2018 Controle              1     1       1
##      Tratamento            1     1       1

xyplot(Sadios+Doentes~Tratamento|Bloco, groups=Ano, data=exper3, type="o",auto.key = T)

with(exper3,
     interaction.plot(Bloco, Tratamento, Sadios,
                      ylab='Sadios',
                      xlab='Cultivar',
                      col=1:4,
                      lwd=2))                      

with(exper3,
     interaction.plot(Tratamento, Bloco, Sadios,
                      ylab='Sadios',
                      xlab='Cultivar',
                      col=1:4,
                      lwd=2))                      

bwplot(Sadios ~ Tratamento|Ano, data=exper3)

barchart(
  Sadios ~ Tratamento | Ano,
  groups = Bloco,
  auto.key = list(space = "right", title = "Ano"),
  exper3,
  aspect = "fill",
  ylab = "Prod",
  main = "Frutos sadios"
)

# opa <- aov(Sadios ~ Bloco + Ano*Tratamento + Error(Bloco/Ano), exper3)
# av.a <- with(exper3,sp.plot(Bloco,Ano,Tratamento,Sadios))
# 
# qqnorm(exper3$Sadios)
# qqline(exper3$Sadios)
# shapiro.test(exper3$Sadios)
# 
# ## Efeito principal: Ano
# 
# tk1 <- TukeyC(opa,
#               which='Ano',
#               error='Bloco:Ano',
#               sig.level=.05)
# 
# summary(tk1)
# 
# ## Efeito principal: Tratamento
# 
# tk2 <- TukeyC(opa,
#               which='Tratamento',
#               error='Within',
#               sig.level=0.05)
# summary(tk2)
# 
# 
# tkn3 <- TukeyC(opa,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=1)
# 
# summary(tkn3)
# 
# tkn4 <- TukeyC(opa,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=2)
# 
# summary(tkn4)
# 
# library(ScottKnott)
# 
# sk1 <- SK(opa,
#           which='Ano',
#           error='Bloco:Ano',
#           sig.level=.05)
# 
# summary(sk1)
# 
# 
# sk2 <- SK(opa,
#           which='Tratamento',
#           error='Within',
#           sig.level=0.05)
# summary(sk2)
# 
# 
# sk3 <- SK.nest(opa,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=1)
# 
# summary(sk3)
# 
# sk4 <- SK.nest(opa,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=2)
# 
# summary(sk4)
# 
# 
# opp <- aov(Doentes ~ Bloco + Ano*Tratamento + Error(Bloco/Ano), exper3)
# summary(opp)
# 
# qqnorm(exper3$Doentes)
# qqline(exper3$Doentes)
# shapiro.test(exper3$Doentes)
# 
# ## Efeito principal: Ano
# 
# tk5 <- TukeyC(opp,
#               which='Ano',
#               error='Bloco:Ano',
#               sig.level=.05)
# 
# summary(tk5)
# 
# ## Efeito principal: Tratamento
# 
# tk5 <- TukeyC(opp,
#               which='Tratamento',
#               error='Within',
#               sig.level=0.05)
# summary(tk5)
# 
# 
# tkn5 <- TukeyC(opp,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=1)
# 
# summary(tkn5)
# 
# tkn6 <- TukeyC(opp,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=2)
# 
# summary(tkn5)
# 
# library(ScottKnott)
# 
# sk4 <- SK(opp,
#           which='Ano',
#           error='Bloco:Ano',
#           sig.level=.05)
# 
# summary(sk4)
# 
# 
# sk5 <- SK(opp,
#           which='Tratamento',
#           error='Within',
#           sig.level=0.05)
# 
# summary(sk5)
# 
# 
# sk6 <- SK.nest(opp,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=1)
# 
# summary(sk6)
# 
# sk7 <- SK.nest(opp,
#                which='Ano:Tratamento',
#                error='Within',
#                fl1=2)
# 
# summary(sk7)




library(bpca)

#par(mfrow=c(1,2))  

vvvv <- subset(exper3, Ano %in% 2017)

rownames(vvvv) <- c("Camaca-C","Camaca-T","Cepec-C","Cepec-T", "Itubera-C", "Itubera-T")

plot(
  bpca(vvvv[, 4:5]),
  type = "bp",
  var.f = .8,
  var.cex = 1.2,
  var.color = c("blue", "red", "green", "black"),
  obj.cex = 1,
  obj.names = T,
  obj.id = 12,
  var.id = 1,
  main="2017"
)

zzzz <- subset(exper3, Ano %in% 2018)

rownames(zzzz) <- c("Camaca-C","Camaca-T","Cepec-C","Cepec-T", "Itubera-C", "Itubera-T")

plot(
  bpca(zzzz[, 4:5], scale = T),
  type = "bp",
  var.f = .8,
  var.cex = 1,
  var.color = c("blue", "red", "green", "black"),
  obj.cex = 1,
  obj.names = T,
  obj.id = 16,
  var.id = 1,
  main="2018"
)

#par(mfrow=c(1,2))  



resultado.pca<-rda(vvvv[,4:5], scale=T)

#Avaliando os resultados
summary(resultado.pca)

## 
## Call:
## rda(X = vvvv[, 4:5], scale = T) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total               2          1
## Unconstrained       2          1
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                         PC1    PC2
## Eigenvalue            1.568 0.4319
## Proportion Explained  0.784 0.2160
## Cumulative Proportion 0.784 1.0000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  1.778279 
## 
## 
## Species scores
## 
##           PC1     PC2
## Sadios  1.113 -0.5843
## Doentes 1.113  0.5843
## 
## 
## Site scores (weighted sums of species scores)
## 
##                PC1     PC2
## Camaca-C  -0.78277 -0.3018
## Camaca-T  -0.98927 -0.2920
## Cepec-C    0.15205  1.2033
## Cepec-T    0.34863  0.7620
## Itubera-C  0.07977 -0.5941
## Itubera-T  1.19158 -0.7774

#Gráfico
biplot(resultado.pca, main="2017")

###

#Diferenciando os ambientes no resultado gráfico:
simbolos<-c(15, 17)
biplot(resultado.pca, display = "species", xlim=c(-2, 2), ylim=c(-1,1), main="2017")
points(resultado.pca, pch=simbolos[vvvv$Tratamento])
legend("topleft", c("Controle","Tratamento"), cex=0.8, pch = c(15, 17), inset = .02)


#Com o gráfico aberto, podemos conectar os pontos:
ordihull(resultado.pca, vvvv$Tratamento)

###

#Salvando os componentes como variáveis:
componentes<-scores(resultado.pca, display="sites", choices=1:2)

pc1<-componentes[ , 1]
pc2<-componentes[ , 2]

boxplot(pc1 ~ vvvv$Tratamento)

boxplot(pc2 ~ vvvv$Tratamento)

resultado.pca<-rda(zzzz[,4:5], scale=T)

#Avaliando os resultados
summary(resultado.pca)

## 
## Call:
## rda(X = zzzz[, 4:5], scale = T) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total               2          1
## Unconstrained       2          1
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                          PC1     PC2
## Eigenvalue            1.9505 0.04952
## Proportion Explained  0.9752 0.02476
## Cumulative Proportion 0.9752 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  1.778279 
## 
## 
## Species scores
## 
##           PC1     PC2
## Sadios  1.242 -0.1979
## Doentes 1.242  0.1979
## 
## 
## Site scores (weighted sums of species scores)
## 
##               PC1     PC2
## Camaca-C  -0.2449  0.4043
## Camaca-T  -0.1375 -1.1226
## Cepec-C   -0.7435  0.2736
## Cepec-T   -0.7519 -0.1614
## Itubera-C  0.6210  1.1557
## Itubera-T  1.2568 -0.5497

#Gráfico
biplot(resultado.pca, main="2018")

###

#Diferenciando os ambientes no resultado gráfico:
simbolos<-c(15, 17)
biplot(resultado.pca, display = "species", xlim=c(-2, 2), ylim=c(-1,1), main="2018")
points(resultado.pca, pch=simbolos[zzzz$Tratamento])
legend("topleft", c("Controle","Tratamento"), cex=0.8,pch = c(15, 17), inset = .02)


#Com o gráfico aberto, podemos conectar os pontos:
ordihull(resultado.pca, vvvv$Tratamento)

###

#Salvando os componentes como variáveis:
componentes<-scores(resultado.pca, display="sites", choices=1:2)

pc1<-componentes[ , 1]
pc2<-componentes[ , 2]

boxplot(pc1 ~ zzzz$Tratamento)

boxplot(pc2 ~ zzzz$Tratamento)