[R-br] Teste de t^2 multivariado de Hotellings

[Alexandre dos Santos]

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

 

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20110415/94a1b955/attachment.html>