Teste de t^2 multivariado de Hotellings

Boa tarde pessoal, Estava tentando comparar dois grupos de indivíduos com 5 variáveis, usando o teste de t^2 de Hotelling, ate ai tudo bem, mais quando simulei os indivíduos e suas variáveis com pelos menos um NA no vetor, a função não conseguia realizar o cálculo, dando a mensagem Erro em solve.default(v) : sistema é computacionalmente singular: condição recíproca número = 0 Tentei inserir na função na.omit(), sem sucesso, alguem poderia me dar um help, obrigado, segue a rotina: #####################Teste de t multivariado de Hotellings############################################ ##Depois de algum teste de cluster (ex. Mahalanobis) vou observar se os dois grupos possuem diferenças ##Grupo 1 ind.pop.1<-1:30 #######Os indivíduos med.pop.1<-data.frame(mma=Mod(rnorm(30)),mmc=Mod(rnorm(30)),mmd=Mod(rnorm(30 )),mca=Mod(rnorm(30)),mca2=Mod(rnorm(30)))#######As medidas dos indivíduos dados.pop.1<-cbind(ind.pop.1,med.pop.1)#######Planilha com os indivíduos e suas medidas head(dados.pop.1) ##Grupinho Oligo 3 ind.pop.3<-1:6 #######Os indivíduos med.pop.3<-data.frame(mma=Mod(rnorm(6)),mmc=Mod(rnorm(6)),mmd=Mod(rnorm(6)), mca=Mod(rnorm(6)))#######As medidas dos indivíduos mca2<-c(1.22222,1.666666,0.4500000,NA ,2.33333,0.789000) dados.pop.3<-cbind(ind.pop.3,med.pop.3,mca2)#######Planilha com os indivíduos e suas medidas head(dados.pop.3) hotelling = function(y1, y2) { ### Função que faz o cálculo de Hotellings by http://sas-and-r.blogspot.com/2010/05/example-737-calculation-of-hotellings. html # check the appropriate dimension k = ncol(y1) k2 = ncol(y2) if (k2!=k) stop("input matrices must have the same number of columns: y1 has ", k, " while y2 has ", k2) # calculate sample size and observed means n1 = nrow(y1) n2 = nrow(y2) ybar1= apply(y1, 2, mean); ybar2= apply(y2, 2, mean) diffbar = ybar1-ybar2 # calculate the variance of the difference in means v = ((n1-1)*var(y1)+ (n2-1)*var(y2)) /(n1+n2-2) # calculate the test statistic and associated quantities t2 = n1*n2*diffbar%*%solve(v)%*%diffbar/(n1+n2) f = (n1+n2-k-1)*t2/((n1+n2-2)*k) pvalue = 1-pf(f, k, n1+n2-k-1) # return the list of results return(list(pvalue=pvalue, f=f, t2=t2, diff=diffbar)) } ## Fazendo o teste entre o grupo 1 e 3 (result= hotelling(dados.pop.1,dados.pop.3)) Alexandre dos Santos Ingenieur forestier, Msc. INRA- Biostatistique et Processus Spatiaux (BioSP) Domaine Saint-Paul Site Agroparc 84914 - Avignon - France Tél. : +33 (0)6 87 95 16 29
participantes (1)
-
Alexandre dos Santos