Análise Multivariada de um Banco de Dados.
Natacha N.
26 de outubro de 2016
R Markdown
Introdução
Inicialmente, para este relatório, estabeleceremos algumas questões, nas quais, as respostas surgirão no decorrer dos “Resultados”. Primeiramente, levemos em conta a seguinte pergunta: As variáveis (X1, X2, X3, X4, X5, X6, X7, X8, X9) contidas no banco de dados selecionado tem alguma relação? Seja ela positiva, negativa? Ou não há relação alguma?. Em segundo lugar: caso haja relação, essas variáveis seguem uma normalidade ? Para respondermos essas questões e, assim chegarmos a uma conclusão, seja ela satisfatória ou não, foram necessários a aplicabilidade de alguns testes. São eles:
Henze-Zirkler’s Multivariate Normality Test: O teste Henze-Zirkler baseia-se numa distância funcional não negativa, que mede a distância entre as duas funções de distribuição.
Shapiro-Wilk’s Normality Test: O teste de Shapiro-Wilk utiliza a hipótese nula, a princípio para verificar se uma amostra \(x(1), ..., x(n)\) veio a partir de uma distribuição normal populacional.
Mardia’s Multivariate Normality Test: Esta função calcula os coeficientes de assimetria e curtose multivariadas, bem como a sua significância estatística correspondente. Também pode calcular a versão corrigida do coeficiente de assimetria para uma amostra de tamanho pequeno (n <20). Para normalidade multivariada, ambos os valores de p das estatísticas de assimetria e curtose deve ser maior do que 0,05.
Royston’s Multivariate Normality Test: O teste de Royston usa a estatística de Shapiro-Wilk/Shapiro-Francia para testar a normalidade multivariada. Se a curtose dos dados for superior a 3, ele utiliza o teste de Shapiro-Francia para as distribuições leptocúrticas, caso contrário, ele usa o teste de Shapiro-Wilk para as distribuições platicúrticas.
Após essas avaliações, foi calculada a MANOVA (Análise de Variância Multivariada), onde essa análise é a generalização da ANOVA (Análise de Variância). Na MANOVA destacam - se os seguintes testes:
- Teste de Pillai:
- Teste de Wilk:
- Teste de Roy:
Teste de Hotelling-Lawley:
Onde, o teste de Wilks é o mais utilizado para testar a hipótese \(Ho\) da MANOVA. Além disso, os outros testes citados acima também são utilizados, porém podem apresentar resultados diferentes para a mesma análise.
Material
O banco de dados utilizado neste trabalho está relacionado com as medidas de mandíbulas para 5 grupos de cães (cães Modernos da Tailândia (TM), Chacais Dourados (CD), Cuons (CU), Lobos Indianos (LI) e Cães Pré-históricos Tailandeses (PH)). Nos “Resultados” veremos o banco de dados propriamente dito, com suas respectivas variáveis e observações. OBS: Na variável Sexo, por exemplo, temos os seguintes significados para os valores mostrados: 1 (um) para macho, 2 (dois) para fêmea e, 0 (zero) para desconhecido.
Resultados
Carregando os pacotes necessários para as análises
require("MVA")## Loading required package: MVA## Loading required package: HSAUR2## Loading required package: toolsrequire("ellipse")## Loading required package: ellipserequire(MVN) ## Loading required package: MVN## sROC 0.1-2 loadedrequire(ggplot2)## Loading required package: ggplot2require(GGally)## Loading required package: GGallyrequire(CCA)## Loading required package: CCA## Loading required package: fda## Loading required package: splines## Loading required package: Matrix##
## Attaching package: 'fda'## The following object is masked from 'package:graphics':
##
## matplot## Loading required package: fields## Loading required package: spam## Loading required package: grid## Spam version 1.4-0 (2016-08-29) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: maps#Lendo os dados para análises
setwd("D:/Unidade D/Documentos/Analise Multivariada I/Trabalho_seg_unid/dados")
dados=read.table("dados.txt",head=TRUE)
#Apresentação do Banco de dados
dados## Sistema X1 X2 X3 X4 X5 X6 X7 X8 X9
## 1 TM 123 10.1 23 23 19 7.8 32 33 5.6
## 2 TM 137 9.6 19 22 19 7.8 32 40 5.8
## 3 TM 121 10.2 18 21 21 7.9 35 38 6.2
## 4 TM 130 10.7 24 22 20 7.9 32 37 5.9
## 5 TM 149 12.0 25 25 21 8.4 35 43 6.6
## 6 TM 125 9.5 23 20 20 7.8 33 37 6.3
## 7 TM 126 9.1 20 22 19 7.5 32 35 5.5
## 8 TM 125 9.7 19 19 19 7.5 32 37 6.2
## 9 TM 121 9.6 22 20 18 7.6 31 35 5.3
## 10 TM 122 8.9 20 20 19 7.6 31 35 5.7
## 11 TM 115 9.3 19 19 20 7.8 33 34 6.5
## 12 TM 112 9.1 19 20 19 6.6 30 33 5.1
## 13 TM 124 9.3 21 21 18 7.1 30 36 5.5
## 14 TM 128 9.6 22 21 19 7.5 32 38 5.8
## 15 TM 130 8.4 23 20 19 7.3 31 40 5.8
## 16 TM 127 10.5 25 23 20 8.7 32 35 6.1
## 17 CD 120 8.2 18 17 18 7.0 32 35 5.2
## 18 CD 107 7.9 17 17 20 7.0 32 34 5.3
## 19 CD 110 8.1 18 16 19 7.1 31 32 4.7
## 20 CD 116 8.5 20 18 18 7.1 32 33 4.7
## 21 CD 114 8.2 19 18 19 7.9 32 33 5.1
## 22 CD 111 8.5 19 16 18 7.1 30 33 5.0
## 23 CD 113 8.5 17 18 19 7.1 30 34 4.6
## 24 CD 117 8.7 20 17 18 7.0 30 34 5.2
## 25 CD 114 9.4 21 19 19 7.5 31 35 5.3
## 26 CD 112 8.2 19 17 19 6.8 30 34 5.1
## 27 CD 110 8.5 18 17 19 7.0 31 33 4.9
## 28 CD 111 7.7 20 18 18 6.7 30 32 4.5
## 29 CD 107 7.2 17 16 17 6.0 28 35 4.7
## 30 CD 108 8.2 18 16 17 6.5 29 33 4.8
## 31 CD 110 7.3 19 15 17 6.1 30 33 4.5
## 32 CD 105 8.3 19 17 17 6.5 29 32 4.5
## 33 CD 107 8.4 18 17 18 6.2 29 31 4.3
## 34 CD 106 7.8 19 18 18 6.2 31 32 4.4
## 35 CD 111 8.4 17 16 18 7.0 30 34 4.7
## 36 CD 111 7.6 19 17 18 6.5 30 35 4.6
## 37 CU 123 9.7 22 21 20 7.8 27 36 6.1
## 38 CU 135 11.8 25 21 23 8.9 31 38 7.1
## 39 CU 138 11.4 25 25 22 9.0 30 38 7.3
## 40 CU 141 10.8 26 25 21 8.1 29 39 6.6
## 41 CU 135 11.2 25 25 21 8.5 29 39 6.7
## 42 CU 136 11.0 22 24 22 8.1 31 39 6.8
## 43 CU 131 10.4 23 23 23 8.7 30 36 6.8
## 44 CU 137 10.6 25 24 21 8.3 28 38 6.5
## 45 CU 135 10.5 25 25 21 8.4 29 39 6.9
## 46 CU 131 10.9 25 24 21 8.5 29 35 6.2
## 47 CU 130 11.3 22 23 21 8.7 29 37 7.0
## 48 CU 144 10.8 24 26 22 8.9 30 42 7.1
## 49 CU 139 10.9 26 23 22 8.7 30 39 6.9
## 50 CU 123 9.8 23 22 20 8.1 26 34 5.6
## 51 CU 137 11.3 27 26 23 8.7 30 39 6.5
## 52 CU 128 10.0 22 23 22 8.7 29 37 6.6
## 53 CU 122 9.9 22 22 20 8.2 26 36 5.7
## 54 LI 167 11.5 29 28 25 9.5 41 45 7.2
## 55 LI 164 12.3 27 26 25 10.0 42 47 7.9
## 56 LI 150 11.5 21 24 25 9.3 41 46 8.5
## 57 LI 145 11.3 28 24 24 9.2 36 41 7.2
## 58 LI 177 12.4 31 27 27 10.5 43 50 7.9
## 59 LI 166 13.4 32 27 26 9.5 40 47 7.3
## 60 LI 164 12.1 27 24 25 9.9 42 45 8.3
## 61 LI 165 12.6 30 26 25 7.7 40 43 7.9
## 62 LI 131 11.8 20 24 23 8.8 38 40 6.5
## 63 LI 163 10.8 27 24 24 9.2 39 48 7.0
## 64 LI 164 10.7 24 23 26 9.5 43 47 7.6
## 65 LI 141 10.4 20 23 23 8.9 38 43 6.0
## 66 LI 148 10.6 26 21 24 8.9 39 40 7.0
## 67 LI 158 10.7 25 25 24 9.8 41 45 7.4
## 68 PH 112 10.1 17 18 19 7.7 31 33 5.8
## 69 PH 115 10.0 18 23 20 7.8 33 36 6.0
## 70 PH 136 11.9 22 25 21 8.5 36 39 7.0
## 71 PH 111 9.9 19 20 18 7.3 29 34 5.3
## 72 PH 130 11.2 23 27 20 9.1 35 35 6.6
## 73 PH 125 10.7 19 26 20 8.4 33 37 6.3
## 74 PH 132 9.6 19 20 19 9.7 35 38 6.6
## 75 PH 121 10.7 21 23 19 7.9 32 35 6.0
## 76 PH 122 9.8 22 23 18 7.9 32 35 6.1
## 77 PH 124 9.5 20 24 19 7.6 32 37 6.0head(dados)## Sistema X1 X2 X3 X4 X5 X6 X7 X8 X9
## 1 TM 123 10.1 23 23 19 7.8 32 33 5.6
## 2 TM 137 9.6 19 22 19 7.8 32 40 5.8
## 3 TM 121 10.2 18 21 21 7.9 35 38 6.2
## 4 TM 130 10.7 24 22 20 7.9 32 37 5.9
## 5 TM 149 12.0 25 25 21 8.4 35 43 6.6
## 6 TM 125 9.5 23 20 20 7.8 33 37 6.3attach(dados)
dados1=data.frame(cbind(X4,X5,X6))
head(dados1)## X4 X5 X6
## 1 23 19 7.8
## 2 22 19 7.8
## 3 21 21 7.9
## 4 22 20 7.9
## 5 25 21 8.4
## 6 20 20 7.8mean(X4)## [1] 21.49351attach(dados1)## The following objects are masked from dados:
##
## X4, X5, X6mean(X5)## [1] 20.49351mean(X6)## [1] 8#-------------------------------RELAÇÃO DAS VARIÁVEIS GRAFICAMENTE------------------------------------------
#Visualizando graficamente
#Retirando a primeira coluna
plot(dados[,-1]) #Verificando a correlaçao
pairs(dados[,-1])
#MATRIZ DE VARIÂNCIAS E COVARIÂNCIAS
S=cov(dados[,-1])
S## X1 X2 X3 X4 X5 X6 X7
## X1 306.31511 20.281869 53.630212 47.157724 39.512987 15.2565789 55.421565
## X2 20.28187 1.968462 3.953213 4.216849 2.869481 1.2147368 3.277085
## X3 53.63021 3.953213 12.839371 9.302290 6.947027 2.6407895 7.145762
## X4 47.15772 4.216849 9.302290 11.411141 6.226931 2.7960526 6.634997
## X5 39.51299 2.869481 6.947027 6.226931 6.200615 2.1763158 7.713944
## X6 15.25658 1.214737 2.640789 2.796053 2.176316 1.0478947 2.757895
## X7 55.42157 3.277085 7.145762 6.634997 7.713944 2.7578947 17.410800
## X8 73.19481 4.610629 11.468558 10.640807 9.627649 3.6000000 14.459159
## X9 15.76514 1.269684 2.731596 2.834706 2.241285 0.9334211 2.756408
## X8 X9
## X1 73.194805 15.7651401
## X2 4.610629 1.2696839
## X3 11.468558 2.7315960
## X4 10.640807 2.8347061
## X5 9.627649 2.2412850
## X6 3.600000 0.9334211
## X7 14.459159 2.7564081
## X8 19.401572 3.7640123
## X9 3.764012 1.0397779R=cor(S)
R## X1 X2 X3 X4 X5 X6 X7
## X1 1.0000000 0.9973665 0.9947229 0.9931756 0.9996113 0.9994007 0.9853891
## X2 0.9973665 1.0000000 0.9971192 0.9987356 0.9963850 0.9986676 0.9733676
## X3 0.9947229 0.9971192 1.0000000 0.9955204 0.9931959 0.9947020 0.9651695
## X4 0.9931756 0.9987356 0.9955204 1.0000000 0.9915450 0.9956434 0.9635248
## X5 0.9996113 0.9963850 0.9931959 0.9915450 1.0000000 0.9990162 0.9871172
## X6 0.9994007 0.9986676 0.9947020 0.9956434 0.9990162 1.0000000 0.9824499
## X7 0.9853891 0.9733676 0.9651695 0.9635248 0.9871172 0.9824499 1.0000000
## X8 0.9990410 0.9940102 0.9897529 0.9888114 0.9990336 0.9978383 0.9898187
## X9 0.9994683 0.9987788 0.9950473 0.9957280 0.9990497 0.9998881 0.9817481
## X8 X9
## X1 0.9990410 0.9994683
## X2 0.9940102 0.9987788
## X3 0.9897529 0.9950473
## X4 0.9888114 0.9957280
## X5 0.9990336 0.9990497
## X6 0.9978383 0.9998881
## X7 0.9898187 0.9817481
## X8 1.0000000 0.9979308
## X9 0.9979308 1.0000000#Construindo as elipses de confiança - caso bivariado
Test1<-matrix(c(dados$X4,dados$X5),77,2) #77 é o numero de linhas e 2 dois o número de colunas
Test2<-matrix(c(dados$X4,dados$X6),77,2)
Test3<-matrix(c(dados$X5,dados$X6),77,2)
X11()
bvbox(Test1,xlab="X4",ylab="X5")
X11()
bvbox(Test2,xlab="X4",ylab="X6")
X11()
bvbox(Test3,xlab="X5",ylab="X6")
#---------------------DADOS BIVARIADOS---------------
# Estrutura multivariada dos dados - Bivariada
Test1## [,1] [,2]
## [1,] 23 19
## [2,] 22 19
## [3,] 21 21
## [4,] 22 20
## [5,] 25 21
## [6,] 20 20
## [7,] 22 19
## [8,] 19 19
## [9,] 20 18
## [10,] 20 19
## [11,] 19 20
## [12,] 20 19
## [13,] 21 18
## [14,] 21 19
## [15,] 20 19
## [16,] 23 20
## [17,] 17 18
## [18,] 17 20
## [19,] 16 19
## [20,] 18 18
## [21,] 18 19
## [22,] 16 18
## [23,] 18 19
## [24,] 17 18
## [25,] 19 19
## [26,] 17 19
## [27,] 17 19
## [28,] 18 18
## [29,] 16 17
## [30,] 16 17
## [31,] 15 17
## [32,] 17 17
## [33,] 17 18
## [34,] 18 18
## [35,] 16 18
## [36,] 17 18
## [37,] 21 20
## [38,] 21 23
## [39,] 25 22
## [40,] 25 21
## [41,] 25 21
## [42,] 24 22
## [43,] 23 23
## [44,] 24 21
## [45,] 25 21
## [46,] 24 21
## [47,] 23 21
## [48,] 26 22
## [49,] 23 22
## [50,] 22 20
## [51,] 26 23
## [52,] 23 22
## [53,] 22 20
## [54,] 28 25
## [55,] 26 25
## [56,] 24 25
## [57,] 24 24
## [58,] 27 27
## [59,] 27 26
## [60,] 24 25
## [61,] 26 25
## [62,] 24 23
## [63,] 24 24
## [64,] 23 26
## [65,] 23 23
## [66,] 21 24
## [67,] 25 24
## [68,] 18 19
## [69,] 23 20
## [70,] 25 21
## [71,] 20 18
## [72,] 27 20
## [73,] 26 20
## [74,] 20 19
## [75,] 23 19
## [76,] 23 18
## [77,] 24 19result1 <- hzTest(Test1, qqplot=TRUE)
result1## Henze-Zirkler's Multivariate Normality Test
## ---------------------------------------------
## data : Test1
##
## HZ : 2.194819
## p-value : 2.688081e-05
##
## Result : Data are not multivariate normal.
## ---------------------------------------------#Densidade bivariada tridimencional
mvnPlot(result1, type = "persp", default = TRUE) 
Test2## [,1] [,2]
## [1,] 23 7.8
## [2,] 22 7.8
## [3,] 21 7.9
## [4,] 22 7.9
## [5,] 25 8.4
## [6,] 20 7.8
## [7,] 22 7.5
## [8,] 19 7.5
## [9,] 20 7.6
## [10,] 20 7.6
## [11,] 19 7.8
## [12,] 20 6.6
## [13,] 21 7.1
## [14,] 21 7.5
## [15,] 20 7.3
## [16,] 23 8.7
## [17,] 17 7.0
## [18,] 17 7.0
## [19,] 16 7.1
## [20,] 18 7.1
## [21,] 18 7.9
## [22,] 16 7.1
## [23,] 18 7.1
## [24,] 17 7.0
## [25,] 19 7.5
## [26,] 17 6.8
## [27,] 17 7.0
## [28,] 18 6.7
## [29,] 16 6.0
## [30,] 16 6.5
## [31,] 15 6.1
## [32,] 17 6.5
## [33,] 17 6.2
## [34,] 18 6.2
## [35,] 16 7.0
## [36,] 17 6.5
## [37,] 21 7.8
## [38,] 21 8.9
## [39,] 25 9.0
## [40,] 25 8.1
## [41,] 25 8.5
## [42,] 24 8.1
## [43,] 23 8.7
## [44,] 24 8.3
## [45,] 25 8.4
## [46,] 24 8.5
## [47,] 23 8.7
## [48,] 26 8.9
## [49,] 23 8.7
## [50,] 22 8.1
## [51,] 26 8.7
## [52,] 23 8.7
## [53,] 22 8.2
## [54,] 28 9.5
## [55,] 26 10.0
## [56,] 24 9.3
## [57,] 24 9.2
## [58,] 27 10.5
## [59,] 27 9.5
## [60,] 24 9.9
## [61,] 26 7.7
## [62,] 24 8.8
## [63,] 24 9.2
## [64,] 23 9.5
## [65,] 23 8.9
## [66,] 21 8.9
## [67,] 25 9.8
## [68,] 18 7.7
## [69,] 23 7.8
## [70,] 25 8.5
## [71,] 20 7.3
## [72,] 27 9.1
## [73,] 26 8.4
## [74,] 20 9.7
## [75,] 23 7.9
## [76,] 23 7.9
## [77,] 24 7.6result2 <- hzTest(Test2, qqplot=TRUE)
result2## Henze-Zirkler's Multivariate Normality Test
## ---------------------------------------------
## data : Test2
##
## HZ : 1.017776
## p-value : 0.02962667
##
## Result : Data are not multivariate normal.
## ---------------------------------------------#Densidade bivariada tridimencional
mvnPlot(result2, type = "persp", default = TRUE) 
Test3## [,1] [,2]
## [1,] 19 7.8
## [2,] 19 7.8
## [3,] 21 7.9
## [4,] 20 7.9
## [5,] 21 8.4
## [6,] 20 7.8
## [7,] 19 7.5
## [8,] 19 7.5
## [9,] 18 7.6
## [10,] 19 7.6
## [11,] 20 7.8
## [12,] 19 6.6
## [13,] 18 7.1
## [14,] 19 7.5
## [15,] 19 7.3
## [16,] 20 8.7
## [17,] 18 7.0
## [18,] 20 7.0
## [19,] 19 7.1
## [20,] 18 7.1
## [21,] 19 7.9
## [22,] 18 7.1
## [23,] 19 7.1
## [24,] 18 7.0
## [25,] 19 7.5
## [26,] 19 6.8
## [27,] 19 7.0
## [28,] 18 6.7
## [29,] 17 6.0
## [30,] 17 6.5
## [31,] 17 6.1
## [32,] 17 6.5
## [33,] 18 6.2
## [34,] 18 6.2
## [35,] 18 7.0
## [36,] 18 6.5
## [37,] 20 7.8
## [38,] 23 8.9
## [39,] 22 9.0
## [40,] 21 8.1
## [41,] 21 8.5
## [42,] 22 8.1
## [43,] 23 8.7
## [44,] 21 8.3
## [45,] 21 8.4
## [46,] 21 8.5
## [47,] 21 8.7
## [48,] 22 8.9
## [49,] 22 8.7
## [50,] 20 8.1
## [51,] 23 8.7
## [52,] 22 8.7
## [53,] 20 8.2
## [54,] 25 9.5
## [55,] 25 10.0
## [56,] 25 9.3
## [57,] 24 9.2
## [58,] 27 10.5
## [59,] 26 9.5
## [60,] 25 9.9
## [61,] 25 7.7
## [62,] 23 8.8
## [63,] 24 9.2
## [64,] 26 9.5
## [65,] 23 8.9
## [66,] 24 8.9
## [67,] 24 9.8
## [68,] 19 7.7
## [69,] 20 7.8
## [70,] 21 8.5
## [71,] 18 7.3
## [72,] 20 9.1
## [73,] 20 8.4
## [74,] 19 9.7
## [75,] 19 7.9
## [76,] 18 7.9
## [77,] 19 7.6result3 <- hzTest(Test3, qqplot=TRUE)
result3## Henze-Zirkler's Multivariate Normality Test
## ---------------------------------------------
## data : Test3
##
## HZ : 2.115799
## p-value : 4.143491e-05
##
## Result : Data are not multivariate normal.
## ---------------------------------------------#Densidade bivariada tridimencional
mvnPlot(result3, type = "persp", default = TRUE) 
####------- Testes para análise univariada-----
#Construindo os QQ-plots Univariados
uniPlot(dados[,-1], type = "qqplot") 
#Construindo os histogramas Univariados
uniPlot(dados[,-1], type = "histogram") 
#Construindo os scatter-plots Univariados
uniPlot(dados[,-1], type = "scatter")
#Construindo os Box-Plots Univariados
uniPlot(dados[,-1], type = "box")## Warning in uniPlot(dados[, -1], type = "box"): Box-Plots are based on
## standardized values (centered and scaled).#Testando a normalidade univariada
uniNorm(dados[,-1], type = "SW", desc = TRUE)## $`Descriptive Statistics`
## n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
## X1 77 128.974 17.502 125.0 105.0 177.0 114.0 137.0 0.824 -0.065
## X2 77 9.961 1.403 10.0 7.2 13.4 8.7 10.9 0.034 -0.754
## X3 77 21.948 3.583 22.0 17.0 32.0 19.0 25.0 0.677 -0.237
## X4 77 21.494 3.378 22.0 15.0 28.0 18.0 24.0 -0.148 -1.133
## X5 77 20.494 2.490 20.0 17.0 27.0 19.0 22.0 0.760 -0.396
## X6 77 8.000 1.024 7.9 6.0 10.5 7.1 8.7 0.145 -0.641
## X7 77 32.519 4.173 31.0 26.0 43.0 30.0 33.0 1.086 0.231
## X8 77 37.403 4.405 36.0 31.0 50.0 34.0 39.0 0.973 0.190
## X9 77 6.075 1.020 6.1 4.3 8.5 5.3 6.8 0.169 -0.748
##
## $`Shapiro-Wilk's Normality Test`
## Variable Statistic p-value Normality
## 1 X1 0.9208 0.0001 NO
## 2 X2 0.9824 0.3672 YES
## 3 X3 0.9359 0.0008 NO
## 4 X4 0.9538 0.0070 NO
## 5 X5 0.9116 0.0001 NO
## 6 X6 0.9861 0.5705 YES
## 7 X7 0.8588 0.0000 NO
## 8 X8 0.9047 0.0000 NO
## 9 X9 0.9766 0.1675 YESrequire(MASS)## Loading required package: MASSsolo1<-dados$X4
solo1## [1] 23 22 21 22 25 20 22 19 20 20 19 20 21 21 20 23 17 17 16 18 18 16 18
## [24] 17 19 17 17 18 16 16 15 17 17 18 16 17 21 21 25 25 25 24 23 24 25 24
## [47] 23 26 23 22 26 23 22 28 26 24 24 27 27 24 26 24 24 23 23 21 25 18 23
## [70] 25 20 27 26 20 23 23 24boxcox(solo1~1)
solo2<-dados$X5
boxcox(solo2~1)
boxcox(solo2~1,lam=c(-2,2,0.5))
abline(v=0.15,col="red")
solo3<-dados$X6
boxcox(solo3~1)
nova<-(dados$X6)**0.5
boxcox(nova~1)
# Testando a normalidade multivariada:
# Mardia's Multivariate Normality Test
r1 <- mardiaTest(dados[,-1], qqplot = TRUE)
r1## Mardia's Multivariate Normality Test
## ---------------------------------------
## data : dados[, -1]
##
## g1p : 24.93864
## chi.skew : 320.0458
## p.value.skew : 5.346495e-12
##
## g2p : 111.7058
## z.kurtosis : 3.961722
## p.value.kurt : 7.441108e-05
##
## chi.small.skew : 335.0928
## p.value.small : 1.169286e-13
##
## Result : Data are not multivariate normal.
## ---------------------------------------# Henze-Zirkler's Multivariate Normality Test
r2 <- hzTest(dados[,-1], qqplot =TRUE)
r2## Henze-Zirkler's Multivariate Normality Test
## ---------------------------------------------
## data : dados[, -1]
##
## HZ : 1.142696
## p-value : 0
##
## Result : Data are not multivariate normal.
## ---------------------------------------------# Royston's Multivariate Normality Test
r3<-roystonTest(dados[,-1], qqplot = TRUE)
r3## Royston's Multivariate Normality Test
## ---------------------------------------------
## data : dados[, -1]
##
## H : 42.32811
## p-value : 1.582013e-08
##
## Result : Data are not multivariate normal.
## ---------------------------------------------#-------------verificando se há outiler----------------------
# Verificando os outlier's no conjunto de dados--- caso bivariado
test1 <- dados[,5:6]
# Mahalanobis distance
result <- mvOutlier(test1, qqplot = TRUE, method = "quan")
# Adjusted Mahalanobis distance
result <- mvOutlier(test1, qqplot = TRUE, method = "adj.quan")
# Verificando outlier no conjunto de dados----multivariado
test2 <- dados[,5:7] #considerando as três variáveis
# Mahalanobis distance
result <- mvOutlier(test2, qqplot = TRUE, method = "quan")
# Adjusted Mahalanobis distance
result <- mvOutlier(test2, qqplot = TRUE, method = "adj.quan")
####Verificando, através da função "ellipse", se há relações entre as variáveis. Caso bivariado.
library(mvtnorm)
library(ellipse)
###Primeiro com as variáveis X4 e X5
####Construir a elipse
dados<-data.frame(cbind(X4,X5))
S<-cov(dados)
S## X4 X5
## X4 11.411141 6.226931
## X5 6.226931 6.200615media<-apply(dados,2,mean)
media## X4 X5
## 21.49351 20.49351corr<-cov2cor(S)
corr## X4 X5
## X4 1.0000000 0.7402734
## X5 0.7402734 1.0000000# Plotar os dados
max(X4)## [1] 28min(X4)## [1] 15max(X5)## [1] 27min(X5)## [1] 17plot(dados,xlim=c(-5,40),ylim=c(-5,40))
# Plotar a ellipse
lines( ellipse(corr, centre = media) , col='red')
abline(h=0)
abline(v=0)
###########Agora, com as variáveis X4 e X6:
library(mvtnorm)
library(ellipse)
#Construir a elipse
dados<-data.frame(cbind(X4,X6))
S<-cov(dados)
S## X4 X6
## X4 11.411141 2.796053
## X6 2.796053 1.047895media<-apply(dados,2,mean)
media## X4 X6
## 21.49351 8.00000corr<-cov2cor(S)
corr## X4 X6
## X4 1.0000000 0.8085781
## X6 0.8085781 1.0000000# Plotar os dados
max(X4)## [1] 28min(X4)## [1] 15max(X6)## [1] 10.5min(X6)## [1] 6plot(dados,xlim=c(-5,40),ylim=c(-5,40))
# Plotar a ellipse
lines( ellipse(corr, centre = media) , col='red')
abline(h=0)
abline(v=0)
########E, por fim, com as variáveis X5 e X6:
library(mvtnorm)
library(ellipse)
#Construir a elipse
dados<-data.frame(cbind(X5,X6))
S<-cov(dados)
S## X5 X6
## X5 6.200615 2.176316
## X6 2.176316 1.047895media<-apply(dados,2,mean)
media## X5 X6
## 20.49351 8.00000corr<-cov2cor(S)
corr## X5 X6
## X5 1.0000000 0.8537794
## X6 0.8537794 1.0000000# Plotar os dados
max(X5)## [1] 27min(X5)## [1] 17max(X6)## [1] 10.5min(X6)## [1] 6plot(dados,xlim=c(-5,40),ylim=c(-5,40))
# Plotar a ellipse
lines( ellipse(corr, centre = media) , col='red')
abline(h=0)
abline(v=0)
Outra alternativa para testar a normalidade—>Teste Shapiro-Wilks
shapiro.test(X4)##
## Shapiro-Wilk normality test
##
## data: X4
## W = 0.95377, p-value = 0.007023shapiro.test(X5)##
## Shapiro-Wilk normality test
##
## data: X5
## W = 0.91161, p-value = 5.368e-05shapiro.test(X6)##
## Shapiro-Wilk normality test
##
## data: X6
## W = 0.98614, p-value = 0.5705Multivariate Analysis of Variance (MANOVA)
# Lendo os dados
setwd("D:/Unidade D/Documentos/Analise Multivariada I/Trabalho_seg_unid/dados")
sistema<-read.table("dados.txt",header=T)
head(sistema)## Sistema X1 X2 X3 X4 X5 X6 X7 X8 X9
## 1 TM 123 10.1 23 23 19 7.8 32 33 5.6
## 2 TM 137 9.6 19 22 19 7.8 32 40 5.8
## 3 TM 121 10.2 18 21 21 7.9 35 38 6.2
## 4 TM 130 10.7 24 22 20 7.9 32 37 5.9
## 5 TM 149 12.0 25 25 21 8.4 35 43 6.6
## 6 TM 125 9.5 23 20 20 7.8 33 37 6.3# Attaching the data frame:
attach(sistema)## The following objects are masked from dados1:
##
## X4, X5, X6## The following objects are masked from dados:
##
## Sistema, X1, X2, X3, X4, X5, X6, X7, X8, X9names(sistema)## [1] "Sistema" "X1" "X2" "X3" "X4" "X5" "X6"
## [8] "X7" "X8" "X9"sistema.manova<-manova(cbind(X4, X5, X6)~Sistema)
summary(sistema.manova)## Df Pillai approx F num Df den Df Pr(>F)
## Sistema 4 1.4727 17.357 12 216 < 2.2e-16 ***
## Residuals 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Análise de variância Multivariada utilizando a estatística Lambda de Wilks
summary(sistema.manova, test="Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## Sistema 4 0.054688 30.906 12 185.49 < 2.2e-16 ***
## Residuals 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#outras estatísticas:
summary(sistema.manova, test="Roy")## Df Roy approx F num Df den Df Pr(>F)
## Sistema 4 6.3897 115.02 4 72 < 2.2e-16 ***
## Residuals 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(sistema.manova, test="Hotelling-Lawley")## Df Hotelling-Lawley approx F num Df den Df Pr(>F)
## Sistema 4 7.8341 44.828 12 206 < 2.2e-16 ***
## Residuals 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(sistema.manova, test="Pillai")## Df Pillai approx F num Df den Df Pr(>F)
## Sistema 4 1.4727 17.357 12 216 < 2.2e-16 ***
## Residuals 72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-------------------Checando graficamente a suposição de normalidade da multivariada-----------------------------------------
chisplot <- function(x) {
if (!is.matrix(x)) stop("x is not a matrix")
### determinando as dimensões
n <- nrow(x)
p <- ncol(x)
xbar <- apply(x, 2, mean)
S <- var(x)
S <- solve(S)
index <- (1:n)/(n+1)
xcent <- t(t(x) - xbar)
di <- apply(xcent, 1, function(x,S) x %*% S %*% x,S)
quant <- qchisq(index,p)
plot(quant, sort(di), ylab = "Distâncias Ordenadas",
xlab = "Quantis da Qui-quadrado", lwd=2,pch=1)
}
###
###############
chisplot(residuals(sistema.manova))
#----sobre a homogeneidade de variâncias e covariâncias de cada tratamento(sistema)
# Examinando as matrizes de variâncias e covariâncias de cada grupo:
by(sistema[,-1], Sistema, var) #matriz de variâncias E COVARIâNCIAS DE CADA SISTEMA## Sistema: CD
## X1 X2 X3 X4 X5 X6
## X1 15.0526316 0.80000000 1.52631579 1.10526316 0.68421053 1.11578947
## X2 0.8000000 0.25326316 0.19684211 0.23684211 0.15684211 0.15663158
## X3 1.5263158 0.19684211 1.30526316 0.52631579 -0.07368421 0.12210526
## X4 1.1052632 0.23684211 0.52631579 0.94736842 0.36842105 0.21578947
## X5 0.6842105 0.15684211 -0.07368421 0.36842105 0.69473684 0.26526316
## X6 1.1157895 0.15663158 0.12210526 0.21578947 0.26526316 0.23081579
## X7 2.2631579 0.14947368 0.25263158 0.47368421 0.61052632 0.36289474
## X8 2.1578947 0.02842105 -0.11578947 0.05263158 0.13684211 0.12078947
## X9 0.6473684 0.06905263 0.04947368 0.05263158 0.13578947 0.09939474
## X7 X8 X9
## X1 2.2631579 2.15789474 0.64736842
## X2 0.1494737 0.02842105 0.06905263
## X3 0.2526316 -0.11578947 0.04947368
## X4 0.4736842 0.05263158 0.05263158
## X5 0.6105263 0.13684211 0.13578947
## X6 0.3628947 0.12078947 0.09939474
## X7 1.2921053 0.13421053 0.17184211
## X8 0.1342105 1.39736842 0.21921053
## X9 0.1718421 0.21921053 0.09734211
## --------------------------------------------------------
## Sistema: CU
## X1 X2 X3 X4 X5 X6 X7
## X1 41.316176 2.6191176 7.2977941 7.4007353 3.5698529 1.0404412 7.06250
## X2 2.619118 0.3706618 0.5610294 0.3713235 0.3694853 0.1309191 0.70625
## X3 7.297794 0.5610294 2.8088235 1.4595588 0.5955882 0.1882353 0.93750
## X4 7.400735 0.3713235 1.4595588 2.4926471 0.4264706 0.1643382 0.93750
## X5 3.569853 0.3694853 0.5955882 0.4264706 1.0147059 0.2496324 1.31250
## X6 1.040441 0.1309191 0.1882353 0.1643382 0.2496324 0.1173529 0.32500
## X7 7.062500 0.7062500 0.9375000 0.9375000 1.3125000 0.3250000 2.25000
## X8 10.636029 0.5823529 1.3933824 2.0147059 0.8345588 0.2400735 1.75000
## X9 2.278309 0.2103309 0.2367647 0.3044118 0.3253676 0.1151471 0.61250
## X8 X9
## X1 10.6360294 2.2783088
## X2 0.5823529 0.2103309
## X3 1.3933824 0.2367647
## X4 2.0147059 0.3044118
## X5 0.8345588 0.3253676
## X6 0.2400735 0.1151471
## X7 1.7500000 0.6125000
## X8 3.7205882 0.6286765
## X9 0.6286765 0.2286029
## --------------------------------------------------------
## Sistema: LI
## X1 X2 X3 X4 X5 X6
## X1 156.401099 4.84670330 37.1483516 14.6483516 11.9560440 3.80164835
## X2 4.846703 0.80489011 2.1203297 1.1703297 0.5703297 0.05851648
## X3 37.148352 2.12032967 14.9505495 4.6043956 2.8351648 0.52252747
## X4 14.648352 1.17032967 4.6043956 3.6043956 1.1428571 0.36483516
## X5 11.956044 0.57032967 2.8351648 1.1428571 1.2967033 0.37252747
## X6 3.801648 0.05851648 0.5225275 0.3648352 0.3725275 0.44554945
## X7 18.994505 0.46648352 1.8736264 1.3736264 1.7582418 0.77637363
## X8 30.543956 0.78736264 4.5879121 2.9340659 2.5494505 1.30054945
## X9 4.920330 0.27170330 0.9060440 0.3637363 0.5175824 0.14587912
## X7 X8 X9
## X1 18.9945055 30.5439560 4.9203297
## X2 0.4664835 0.7873626 0.2717033
## X3 1.8736264 4.5879121 0.9060440
## X4 1.3736264 2.9340659 0.3637363
## X5 1.7582418 2.5494505 0.5175824
## X6 0.7763736 1.3005495 0.1458791
## X7 4.1813187 4.5109890 0.9214286
## X8 4.5109890 9.2582418 0.9785714
## X9 0.9214286 0.9785714 0.4607143
## --------------------------------------------------------
## Sistema: PH
## X1 X2 X3 X4 X5 X6
## X1 70.844444 3.3311111 10.3333333 13.8666667 4.5111111 4.8755556
## X2 3.331111 0.5937778 0.8000000 1.2822222 0.5533333 0.1571111
## X3 10.333333 0.8000000 3.7777778 3.7777778 0.3333333 0.4333333
## X4 13.866667 1.2822222 3.7777778 8.1000000 1.5888889 0.5433333
## X5 4.511111 0.5533333 0.3333333 1.5888889 0.9000000 0.2811111
## X6 4.875556 0.1571111 0.4333333 0.5433333 0.2811111 0.5498889
## X7 15.844444 0.9088889 1.8888889 3.2000000 1.5111111 1.2755556
## X8 12.533333 0.4044444 0.8888889 2.3222222 1.0333333 0.7211111
## X9 3.782222 0.2268889 0.4888889 0.7522222 0.3322222 0.2796667
## X7 X8 X9
## X1 15.8444444 12.5333333 3.7822222
## X2 0.9088889 0.4044444 0.2268889
## X3 1.8888889 0.8888889 0.4888889
## X4 3.2000000 2.3222222 0.7522222
## X5 1.5111111 1.0333333 0.3322222
## X6 1.2755556 0.7211111 0.2796667
## X7 4.4000000 2.8666667 0.9822222
## X8 2.8666667 3.4333333 0.6633333
## X9 0.9822222 0.6633333 0.2290000
## --------------------------------------------------------
## Sistema: TM
## X1 X2 X3 X4 X5 X6 X7
## X1 72.329167 4.3350000 10.2916667 9.8083333 2.4250000 2.1450000 5.2708333
## X2 4.335000 0.7180000 0.9966667 1.0433333 0.4766667 0.2893333 0.7983333
## X3 10.291667 0.9966667 5.3166667 2.2166667 0.3833333 0.5833333 0.3083333
## X4 9.808333 1.0433333 2.2166667 2.6500000 0.5500000 0.4633333 0.9250000
## X5 2.425000 0.4766667 0.3833333 0.5500000 0.7833333 0.2700000 1.1083333
## X6 2.145000 0.2893333 0.5833333 0.4633333 0.2700000 0.2313333 0.4616667
## X7 5.270833 0.7983333 0.3083333 0.9250000 1.1083333 0.4616667 2.0625000
## X8 19.908333 0.9300000 1.6166667 1.6500000 0.9500000 0.4166667 1.9583333
## X9 1.497917 0.1828333 0.1858333 0.1108333 0.2925000 0.1385000 0.4887500
## X8 X9
## X1 19.9083333 1.4979167
## X2 0.9300000 0.1828333
## X3 1.6166667 0.1858333
## X4 1.6500000 0.1108333
## X5 0.9500000 0.2925000
## X6 0.4166667 0.1385000
## X7 1.9583333 0.4887500
## X8 7.4500000 0.5741667
## X9 0.5741667 0.1796250##### As diferenças de médias para fazer comparações múltiplas de Bonferroni:
medias_grupo <- cbind(tapply(X4,Sistema,mean),
tapply(X5,Sistema,mean),
tapply(X6,Sistema,mean))
#------------------------------------------------------------------------------------------------
#Calculando o vetor de médias
SX4<-mean(sistema$X4)
SX5<-mean(sistema$X5)
SX6<-mean(sistema$X6)
medias<-matrix(c(SX4,SX5,SX6),3,1);medias## [,1]
## [1,] 21.49351
## [2,] 20.49351
## [3,] 8.00000#Cálculo da matriz de variâncias e covariâncias
Sp<-cov(sistema[,-1]);Sp## X1 X2 X3 X4 X5 X6 X7
## X1 306.31511 20.281869 53.630212 47.157724 39.512987 15.2565789 55.421565
## X2 20.28187 1.968462 3.953213 4.216849 2.869481 1.2147368 3.277085
## X3 53.63021 3.953213 12.839371 9.302290 6.947027 2.6407895 7.145762
## X4 47.15772 4.216849 9.302290 11.411141 6.226931 2.7960526 6.634997
## X5 39.51299 2.869481 6.947027 6.226931 6.200615 2.1763158 7.713944
## X6 15.25658 1.214737 2.640789 2.796053 2.176316 1.0478947 2.757895
## X7 55.42157 3.277085 7.145762 6.634997 7.713944 2.7578947 17.410800
## X8 73.19481 4.610629 11.468558 10.640807 9.627649 3.6000000 14.459159
## X9 15.76514 1.269684 2.731596 2.834706 2.241285 0.9334211 2.756408
## X8 X9
## X1 73.194805 15.7651401
## X2 4.610629 1.2696839
## X3 11.468558 2.7315960
## X4 10.640807 2.8347061
## X5 9.627649 2.2412850
## X6 3.600000 0.9334211
## X7 14.459159 2.7564081
## X8 19.401572 3.7640123
## X9 3.764012 1.0397779#Cálculo da matriz de correlação
Rp<-cor(sistema[,-1]);Rp## X1 X2 X3 X4 X5 X6 X7
## X1 1.0000000 0.8259623 0.8551719 0.7976348 0.9066471 0.8515578 0.7589012
## X2 0.8259623 1.0000000 0.7863480 0.8897336 0.8213389 0.8457847 0.5597767
## X3 0.8551719 0.7863480 1.0000000 0.7685169 0.7785916 0.7199512 0.4779333
## X4 0.7976348 0.8897336 0.7685169 1.0000000 0.7402734 0.8085781 0.4707245
## X5 0.9066471 0.8213389 0.7785916 0.7402734 1.0000000 0.8537794 0.7424201
## X6 0.8515578 0.8457847 0.7199512 0.8085781 0.8537794 1.0000000 0.6456683
## X7 0.7589012 0.5597767 0.4779333 0.4707245 0.7424201 0.6456683 1.0000000
## X8 0.9494620 0.7460676 0.7266386 0.7151408 0.8777774 0.7984086 0.7867110
## X9 0.8833714 0.8874866 0.7476086 0.8229495 0.8826925 0.8942284 0.6478342
## X8 X9
## X1 0.9494620 0.8833714
## X2 0.7460676 0.8874866
## X3 0.7266386 0.7476086
## X4 0.7151408 0.8229495
## X5 0.8777774 0.8826925
## X6 0.7984086 0.8942284
## X7 0.7867110 0.6478342
## X8 1.0000000 0.8380353
## X9 0.8380353 1.0000000#Calculando a inversa da matriz Sp
Sdinv <- solve(Sp);Sdinv## X1 X2 X3 X4 X5 X6
## X1 0.11575237 0.0609763 -0.18064853 -0.036576225 0.00646519 -0.11462084
## X2 0.06097630 4.2736021 -0.30472698 -0.798923418 -0.32320260 -0.35886445
## X3 -0.18064853 -0.3047270 0.55939635 -0.017785521 -0.20577090 0.24019469
## X4 -0.03657622 -0.7989234 -0.01778552 0.523492320 0.14815877 -0.34085189
## X5 0.00646519 -0.3232026 -0.20577090 0.148158770 1.34005053 -0.57355564
## X6 -0.11462084 -0.3588645 0.24019469 -0.340851895 -0.57355564 5.94925880
## X7 -0.05977183 -0.1129111 0.14880024 0.075700835 -0.13191531 -0.04732921
## X8 -0.22242193 0.1410147 0.28252965 -0.009926485 -0.19769248 0.23699109
## X9 -0.20261792 -2.3567639 0.50070058 -0.028430302 -0.87499465 -2.36250689
## X7 X8 X9
## X1 -0.059771831 -0.222421926 -0.20261792
## X2 -0.112911068 0.141014740 -2.35676394
## X3 0.148800244 0.282529647 0.50070058
## X4 0.075700835 -0.009926485 -0.02843030
## X5 -0.131915306 -0.197692481 -0.87499465
## X6 -0.047329215 0.236991087 -2.36250689
## X7 0.215055636 -0.008976164 0.23607356
## X8 -0.008976164 0.766819347 -0.05371557
## X9 0.236073563 -0.053715569 9.24940057Conclusões
A partir dos resultados obtidos, podemos concluir que, de acordo com plot em pares, há uma relação positiva entre as variáveis. Após isso, temos as matrizes de variâncias e covariâncias, respectivamente, comprovando a relação positiva. Após isto, realizamos as elipses de confiança para os casos bivariados. E, como a autora deste trabalho ficou com três variáveis(X4, X5, X6), ocorreram as seguintes combinações: X4 e X5; X4 e X6 e, X5 e X6. A partir dos boxplots para os casos bivariados, notamos que as observações, com 95% de confiança, não se encontram de maneira uniforme, logo há presença de erros sistemáticos para ambas as combinações. Além disso, há pontos fora da elipse, representando resultados não compatíveis com os demais, presentes nos gráficos 2 (X4 e X6) e 3 (X5 E X6). Iniciando agora, o uso dos testes -caso bivariados- para, de fato, confirmar se há ou não a presença da normalidade. E, segundo o teste de Henze-Zirkler os dados não seguem uma normalidade multivariada, pois o p-valor foi bem menor que o nível de significância estabelecido (5% ou 0,05) para ambas as combinações. E, o que nos faz comprovar ainda mais essa conclusão são o fato de, o gráfico de densidade não ter formato de um “chapéu perfeito”, os qqplot’s se comportarem na forma de escada, principalmente nas variáveis X1, X3, X4, X5, X7 e X8,e, os pontos ou observações não se ajustarem a curva de densidade dos histogramas, com exceção da variável X2, pois foi a que mais se aproximou da curva. Logo após, temos as estatísticas descritivas (tamanhos amostrais, respectivas médias, desvios padrões, mediana, valores mínimos e máximos e etc); teste de Shapiro-Wilk (para os casos univariados) indicando que as variáveis X2, X6 e X9 possuem normalidade; a função boxcox responsável pelos gráficos da verossimilhança, onde os valores máximos da curva devem constar dentro do intervalo para, então, seguir uma normalidade. E, por fim, analisamos a MANOVA.Temos ainda que, O teste de normalidade multivariada de Royston indica que os dados não seguem uma distribuição normal. De acordo com as verficações em relação aos possíveis outliers e considerando a distância de Mahalanobis, podemos observar que, à medida que os quantis aumentam, os outliers diminuem. Depois, utlizou-se o Teste de Shapiro para o caso univariado e, para esse teste, as variáveis indicaram não possuir normalidade. E, por fim, analisamos a MANOVA. De acordo com o p-valor de Wilks (p-valor< 0.05 ou 5%), a MANOVA é estatisticamente significativa. Em outras palavras, há uma diferença estatisticamente significativa entre medidas das mandíbulas dos 5 (cinco) grupos de cães.