Análise multivariada
Juan Pablo Edwards Molina
## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE
##
## [[9]]
## [1] TRUE
##
## [[10]]
## [1] TRUE
##
## [[11]]
## [1] TRUE
##
## [[12]]
## [1] TRUEdatxy<-read.csv("dat_xy_isa.csv", dec=",",sep=",", header=TRUE)
attach(datxy);str(datxy) ## 'data.frame': 84 obs. of 6 variables:
## $ exp : int 1 1 1 1 1 1 2 2 2 2 ...
## $ esp : Factor w/ 2 levels "M. frut","M. laxa": 1 1 1 2 2 2 1 1 1 2 ...
## $ iso : Factor w/ 6 levels "Mf1","Mf2","Mf3",..: 1 2 3 4 5 6 1 2 3 4 ...
## $ var : Factor w/ 7 levels "inc","incub",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ var_id: Factor w/ 7 levels "Incidence","Incubation period",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ val : num 90 72.8 75 80 68.6 ...datxy <- transform(datxy, exp=factor(exp)); str(datxy) ## 'data.frame': 84 obs. of 6 variables:
## $ exp : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 2 2 2 2 ...
## $ esp : Factor w/ 2 levels "M. frut","M. laxa": 1 1 1 2 2 2 1 1 1 2 ...
## $ iso : Factor w/ 6 levels "Mf1","Mf2","Mf3",..: 1 2 3 4 5 6 1 2 3 4 ...
## $ var : Factor w/ 7 levels "inc","incub",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ var_id: Factor w/ 7 levels "Incidence","Incubation period",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ val : num 90 72.8 75 80 68.6 ...levels(datxy$esp)[1]="M. fruticola"
full_dataset <- reshape(datxy[,c(1,2,3,5,6)],
timevar = "var_id",
idvar = c("exp", "esp", "iso"),
direction = "wide")
names(full_dataset)<-c( "Experiment","Specie","Isolate","Incidence","Incubation_period","Latency_period",
"Lesion_size","Sporulation","In_vitro_Growth_rate", "In_vitro_Sporulation" )datxy5<-read.csv("dat_xy5.csv", dec=",",sep=",", header=TRUE)
exp.dat <- reshape(datxy5[,c(1,2,3,4,6)],
timevar = "var",
idvar = c("exp", "esp", "iso"),
direction = "wide")
expdat_graph <- reshape(datxy5[,c(1,2,3,5,6)],
timevar = "var_id",
idvar = c("exp", "esp", "iso"),
direction = "wide")
names(expdat_graph)<-c( "Experiment","Specie","Isolate","Incidence","Incubation_period","Latency_period", "Lesion_size","In_vitro_Growth_rate")- Reduced dataset (12obs x 5var -> without sporulation)
dat = cast(datxy5, iso*esp~ var, mean, value = 'val') # mean values wide
dat.graph=cast(datxy5, iso*esp~ var_id, mean, value = 'val') # mean values wide for graphics
datp = cast(datxy5, iso*esp~ var, mean, value = 'val.p') # Standarized mean values wide
datp.graph=cast(datxy5, iso*esp~ var_id, mean, value = 'val.p')
dap=cbind(datp.graph[,3:7], row.names=levels(datp.graph$iso)) # Standarized mean matrixexp.datp <- reshape(datxy5[,c(1,2,3,4,7)],
timevar = "var",
idvar = c("exp", "esp", "iso"),
direction = "wide")
exp.datp$id= as.factor(paste(exp.datp$iso,exp.datp$exp,sep="."))
names(exp.datp)[4:8]<-c( "inc","incub","lat","les","vitro_gr")
exp.datp= exp.datp[,c(9,4,5,6,7,8)]
expdap=cbind(exp.datp[,2:6], row.names=levels(exp.datp$id))Graf 1
xyplot(val ~ iso|var_id, group=esp, type=c("p","a"), pch=19, # Full graphic
data = datxy, scales = "free",ylab="",xlab="",
auto.key = list(cex=1,x= .78, y = .7, corner = c(0, 0), points = TRUE, pch=19,lines = TRUE))
xyplot(val ~ iso|var_id*exp, group=esp, type=c("p","a"), pch=19, # Exp Split
data = datxy, scales = "free", xlab="", ylab="",
#auto.key = list(cex=1,x= .78, y = .7, corner = c(0, 0), points = TRUE, pch=19,lines = TRUE)
)
Graf 2!
Plot.0 <- ggplot(datxy5, aes( x = var_id, y = val , fill = var_id ) ) +
guides(fill=FALSE)+
theme(legend.position="none")+
geom_boxplot( alpha = 0.6, outlier.colour = c("grey40") , outlier.size=3) +
facet_wrap(~exp) +
theme_bw() +
labs(x="", y="") +
theme(axis.text.x = element_text(colour="black", size = 11, ,angle=45,vjust=1,hjust=1))
Plot.0 
Padronizacao das variaveis
mu<- colMeans(dat) ; mu## inc incub lat les vitro_gr
## 63.8383333 5.5416667 6.3250000 5.3925000 0.7741667s<- diag(var(dat)) ; s## inc incub lat les vitro_gr
## 164.22580667 1.13541667 0.49875000 0.54407750 0.09588417dat$Incidence <- (dat$inc - mu[1])/sqrt(s[1])
dat$Incubation_period <- (dat$incub - mu[2])/sqrt(s[2])
dat$Latency_period <- (dat$lat - mu[3])/sqrt(s[3])
dat$Lesion_size <- (dat$les - mu[4])/sqrt(s[4])
dat$In_vitro_Growth_rate <- (dat$vitro_gr - mu[5])/sqrt(s[5])Graf 2!
Plot.2 <- ggplot(datxy5, aes( x = var_id, y = val.p , fill = var_id ) ) +
guides(fill=FALSE)+
theme(legend.position="none")+
geom_boxplot( alpha = 0.6, outlier.colour = c("grey40") , outlier.size=3) +
theme_bw() +
labs(x="", y="") +
theme(axis.text.x = element_text(colour="black", size = 11, ,angle=45,vjust=1,hjust=1))
Plot.2 
CORRGRAMA
panel.corrgram <-
function(x, y, z, subscripts, at,
level = 0.9, label = FALSE, ...)
{
require("ellipse", quietly = TRUE)
x <- as.numeric(x)[subscripts]
y <- as.numeric(y)[subscripts]
z <- as.numeric(z)[subscripts]
zcol <- level.colors(z, at = at, ...)
for (i in seq(along = z)) {
ell <- ellipse(z[i], level = level, npoints = 50,
scale = c(.2, .2), centre = c(x[i], y[i]))
panel.polygon(ell, col = zcol[i], border = zcol[i], ...)
}
if (label)
panel.text(x = x, y = y, lab = 100 * round(z, 2), cex = 0.8,
col = ifelse(z < 0, "white", "black"))
}cor.dat <- cor(dap[, !sapply(dap, is.factor)], use = "pair")
ord <- order.dendrogram(as.dendrogram(hclust(dist(cor.dat))))
levelplot(cor.dat[ord, ord], at = do.breaks(c(-1.01, 1.01), 20),
xlab = NULL, ylab = NULL, main="",
colorkey = list(space = "top"),
scales = list(x = list(rot = 90)),
panel = panel.corrgram, label = TRUE)
Componentes principais
row.names(dat.graph)=dat.graph$iso
pca<-princomp(dat.graph[,3:7] ,cor=TRUE);
#summary(pca.dap,loadings=TRUE) # Mean data
biplot(pca, col=c("black","red"), ylim=c(-1,1),xlim=c(-1,1), cex=0.9,main="")
abline(h=0,v=0, lty=3)
pca.expdap<-princomp(expdap,cor=TRUE);summary(pca.expdap,loadings=TRUE) # Full data## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 1.8223712 1.0144679 0.65004873 0.44314222
## Proportion of Variance 0.6642073 0.2058290 0.08451267 0.03927501
## Cumulative Proportion 0.6642073 0.8700364 0.95454903 0.99382404
## Comp.5
## Standard deviation 0.175726545
## Proportion of Variance 0.006175964
## Cumulative Proportion 1.000000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## inc 0.383 -0.659 -0.240 0.437 -0.414
## incub -0.457 0.447 -0.415 0.307 -0.570
## lat -0.488 -0.181 0.448 0.669 0.283
## les 0.440 0.397 0.660 0.168 -0.431
## vitro_gr 0.462 0.419 -0.366 0.489 0.488biplot(pca.expdap, col=c("black","red"), ylim=c(-1,1),xlim=c(-1,1), cex=0.9,main="")
Analise de agrupamento
row.names(datp)=datp$iso
d1<-dist(datp[,3:7], method='euclidean')# Matriz de Distancias euclidianas
round(d1,2)## 1 2 3 4 5
## 2 1.63
## 3 1.36 1.21
## 4 2.13 3.08 2.14
## 5 2.99 3.12 2.01 2.51
## 6 4.50 4.53 3.55 3.01 2.32hc <- hclust(d1, "ward") # Metodo de Ward## The "ward" method has been renamed to "ward.D"; note new "ward.D2"hc <- hclust(d1, "ave")
plot(as.dendrogram(hc),horiz=TRUE, xlab="Euclidean distances") # Dendrograma
#title(main="Dendogram")Analise exploratorio de fatores
(scores<-factanal(dap,factors=2, method="mle", rotation="varimax",scores="regression"))##
## Call:
## factanal(x = dap, factors = 2, scores = "regression", rotation = "varimax", method = "mle")
##
## Uniquenesses:
## Incidence Incubation period In vitro Growth rate
## 0.005 0.200 0.005
## Latency period Lesion size
## 0.013 0.352
##
## Loadings:
## Factor1 Factor2
## Incidence 0.178 0.982
## Incubation period -0.370 -0.815
## In vitro Growth rate 0.963 0.261
## Latency period -0.898 -0.426
## Lesion size 0.787 0.170
##
## Factor1 Factor2
## SS loadings 2.520 1.906
## Proportion Var 0.504 0.381
## Cumulative Var 0.504 0.885
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 10.33 on 1 degree of freedom.
## The p-value is 0.00131(cargas=scores$loadings)##
## Loadings:
## Factor1 Factor2
## Incidence 0.178 0.982
## Incubation period -0.370 -0.815
## In vitro Growth rate 0.963 0.261
## Latency period -0.898 -0.426
## Lesion size 0.787 0.170
##
## Factor1 Factor2
## SS loadings 2.520 1.906
## Proportion Var 0.504 0.381
## Cumulative Var 0.504 0.885(dap=cbind(dap,round(scores$scores,2)))## Incidence Incubation period In vitro Growth rate Latency period
## Mf1 1.12814515 -0.6121836 0.72697783 -0.82673738
## Mf2 -0.03382362 -0.5808435 1.01377649 -0.66196693
## Mf3 -0.08557795 -0.2207613 0.55697466 -0.40591803
## Ml1 0.96420147 -0.5527032 -0.70391314 0.37018660
## Ml2 -0.78636461 1.0443261 0.01268569 -0.04267482
## Ml3 -1.18658044 0.9221655 -1.60650153 1.56711057
## Lesion size Factor1 Factor2
## Mf1 0.3675801 0.45 1.16
## Mf2 1.4639318 1.08 -0.23
## Mf3 0.4317712 0.63 -0.20
## Ml1 -0.6510400 -1.07 1.25
## Ml2 -0.8090105 0.30 -0.92
## Ml3 -0.8032326 -1.40 -1.06Analise discriminante
variables = exp.dat[,4:8]
Yp<-as.matrix(variables)
row.names(Yp)=exp.dat$iso
grupo = exp.dat$esp
mandat=manova(Yp~grupo)
test_Wilks<-summary(mandat,test="Wilks")
#Si Pr<0.05, há diferencas entre as medias -> analise discriminante tem sentido- Que variaveis tem maior poder discriminatorio: aquela com menor Pr (p.valor)
(poder_discrim<-summary.aov(mandat))## Response val.inc :
## Df Sum Sq Mean Sq F value Pr(>F)
## grupo 1 253 253.0 0.9598 0.3503
## Residuals 10 2636 263.6
##
## Response val.incub :
## Df Sum Sq Mean Sq F value Pr(>F)
## grupo 1 6.0208 6.0208 3.1345 0.1071
## Residuals 10 19.2083 1.9208
##
## Response val.lat :
## Df Sum Sq Mean Sq F value Pr(>F)
## grupo 1 3.1008 3.10083 5.0488 0.04843 *
## Residuals 10 6.1417 0.61417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response val.les :
## Df Sum Sq Mean Sq F value Pr(>F)
## grupo 1 3.3391 3.3391 3.3162 0.09861 .
## Residuals 10 10.0690 1.0069
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response val.vitro_gr :
## Df Sum Sq Mean Sq F value Pr(>F)
## grupo 1 0.68641 0.68641 17.974 0.001717 **
## Residuals 10 0.38188 0.03819
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Como clasifica la funcion discriminante los datos clasificados inicialmente?
z <- lda(grupo~.,data=variables)
P <- predict(z,variables)
(datos_clasificados<-table(grupo,P$class)) ##
## grupo M. frut M. laxa
## M. frut 6 0
## M. laxa 0 6- Coeficientes de las funciones discriminantes canonicas:
(fun_canonica<-z$scaling)## LD1
## val.inc 0.09497875
## val.incub 1.84420110
## val.lat -2.18238504
## val.les 0.55971924
## val.vitro_gr -9.58172652- Incluir novo isolado
iso_ameixa<-structure(list(inc = 40, incu = 5, lat = 6, les = 8,vitro_gr= 0.6, Grupo = 3),
.Names = c("inc", "incub", "lat", "les","vitro_gr", "Grupo"),
row.names = c("ISO_ameixa"),
class = "data.frame")
attach(iso_ameixa)Classificate new case
# p1<-predict(z, iso_ameixa)
# (nuevos_casos<-p1$class)- Todas as variaveis
names(expdat_graph)## [1] "Experiment" "Specie" "Isolate"
## [4] "Incidence" "Incubation_period" "Latency_period"
## [7] "Lesion_size" "In_vitro_Growth_rate"attach(expdat_graph)
da=cbind(expdat_graph[,4:8], row.names=levels(exp.datp$id))
par(mfrow=c(2,2))
plot(Incidence, Incubation_period, data=expdat_graph, pch=19,
xlab="Incidence",ylab="Incubation period",
col=(1:2)[expdat_graph$Specie])## Warning in plot.window(...): "data" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "data" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "data" is not
## a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "data" is not
## a graphical parameter## Warning in box(...): "data" is not a graphical parameter## Warning in title(...): "data" is not a graphical parameterplot(Latency_period, Lesion_size, data=expdat_graph, pch=19,
xlab="Latency period", ylab="Lesion size",
col=(1:2)[expdat_graph$Specie])## Warning in plot.window(...): "data" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "data" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "data" is not
## a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "data" is not
## a graphical parameter## Warning in box(...): "data" is not a graphical parameter## Warning in title(...): "data" is not a graphical parameterplot(In_vitro_Growth_rate, data=expdat_graph, pch=19,
ylab="In vitro Growth rate",
col=(1:2)[expdat_graph$Specie])## Warning in plot.window(...): "data" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "data" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "data" is not
## a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "data" is not
## a graphical parameter## Warning in box(...): "data" is not a graphical parameter## Warning in title(...): "data" is not a graphical parameterplot(P$x,col=grupo,main='',pch=19, xlim=c(0,16), ylab="Discrimant function",axes=FALSE,xlab="")
abline(h=0,lty=2); box(col="black")
text(P$x, row.names(Yp), cex=0.8, pos=4, col="blue")
axis(2, col="black", col.ticks="black", col.axis="black", cex.axis=0.8)
# plot(c(P$x,p1$x),col=c(grupo,"green"),main='',pch=19, xlim=c(0,16), ylab="Discrimant # function",axes=FALSE,xlab="")
#abline(h=0,lty=2); box(col="black")
#text(c(P$x,p1$x), c(row.names(Yp), "Iso_Ameixa"), cex=0.8, pos=4, col="blue")
#axis(2, col="black", col.ticks="black", col.axis="black", cex.axis=0.8)- Carinhas e estrelinhas —-
layout(1)
stars(dat)
aplpack::faces(dap[,1:5],face.type=1)
## effect of variables:
## modified item Var
## "height of face " "Incidence"
## "width of face " "Incubation period"
## "structure of face" "In vitro Growth rate "
## "height of mouth " "Latency period"
## "width of mouth " "Lesion size"
## "smiling " "Incidence"
## "height of eyes " "Incubation period"
## "width of eyes " "In vitro Growth rate "
## "height of hair " "Latency period"
## "width of hair " "Lesion size"
## "style of hair " "Incidence"
## "height of nose " "Incubation period"
## "width of nose " "In vitro Growth rate "
## "width of ear " "Latency period"
## "height of ear " "Lesion size"