[R-br] Sistema é computacionalmente singular - Matriz/Mahalanobis
Tiago Souza Marçal
tiagosouzamarcal em hotmail.com
Sexta Junho 26 14:01:06 BRT 2015
Boa tarde Andre,
como o Walmes afirmou seu exemplo ainda não é reproduzível. Entretanto, segue um exemplo para servir de inspiração para o seu problema.
Uma dica é você buscar entendimento da lógica de obtenção da matriz de variância e covariância Σ nas documentações da função mahalanobis.dist()
do pacote StatMatch.
O exemplo que trago a seguir é para uma condição experimental então utilizei a matriz de variância e covariância S
(matriz de soma de quadrados e produtos residuais/graus de liberdade dos resíduos). Entretanto, você pode tentar adaptar para sua condição
definindo a matriz de variância e covariância de interesse.
# Dados para o exemplo.
dat <- structure(list(trat = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"), x1 = c(50.039225756996, 51.7825110686358, 48.1384396336807, 50.4489824448722, 48.4220046908771, 44.5221770682025, 52.2588301216718, 54.1412387912868, 48.8071610383588, 53.8184419408507, 53.5005661405145, 48.8892313217333), x2 = c(40.0577191986925, 36.6154386398641, 45.1013819536244, 37.9422687414707, 41.1111204733925, 33.3243892012254, 33.0139305834796, 35.7751641361826, 47.7243822772713, 39.8005702375829, 38.3100567194852, 44.0757780378149), x3 = c(10.4591588838242, 10.6373435622737, 10.3189328544074, 10.0386650378322, 9.96842284038367, 10.6003051066032, 10.8695539960311, 10.143593400873, 9.45082912208145, 9.28742328152705, 9.60271226420169, 10.1566398296013)), .Names = c("trat", "x1", "x2", "x3"), row.names = c(NA, -12L), class = "data.frame")
attach(dat)
# Matriz contendo as variáveis que serão analizadas.
resp = as.matrix( dat[,2:ncol(dat)] )
# Realização da Manova para obtenção da matriz de covariância residual (S).
mav <- manova(resp ~ trat)
quadro = summary(mav,test="Wilks")
glr <- df.residual(mav)
S <- quadro$SS$Residuals/glr
# Obtenção das médias dos tratamentos
medtrat <- aggregate(resp, by = list(trat), FUN = mean)[,2:ncol(dat)]
medtrat <- as.matrix(medtrat)
# Obtenção da matriz D2 de Mahalanobis
nla = numeric(0)
nlb = numeric(0)
acd = (nrow(medtrat)-1):1
acc=1:(nrow(medtrat)-1)
acc1=2:nrow(medtrat)
D2 <- mat.or.vec(nrow(medtrat),nrow(medtrat))
for(i in 1:length(acd)){ nla <- c(nla, rep( acc[i],each=acd[i] )) nlb <- c(nlb, acc1[i]:nrow(medtrat) ) }
for(i in 1:length(nla)){ D2[nla[i], nlb[i]] <- t(medtrat[nla[i],] - medtrat[nlb[i],])%*%solve(S)%*%(medtrat[nla[i],] - medtrat[nlb[i],]) D2[nlb[i], nla[i]] <- D2[nla[i], nlb[i]] }
# Matriz D2 de Mahalanobis
D2
################################################################# Tiago de Souza Marçal - Graduando em Agronomia pelo CCA-UFES Bolsista de Iniciação Científica da área de Genética e Melhoramento de Plantas #################################################################
Date: Fri, 26 Jun 2015 00:34:46 +0000
From: andreolsouza em yahoo.com.br
To: walmeszeviani em gmail.com; r-br em listas.c3sl.ufpr.br
Subject: Re: [R-br] Sistema é computacionalmente singular - Matriz/Mahalanobis
Desculpas pela falta de clareza Walmes, é que não citei que a função faz é D_{ij}^2 = (x_i - x_j)' (Σ^-1) (x_i - x_j) e a matriz Σ, neste caso não tem inversa comum. Já tentei várias funções e de fato todas funções para tal método usam a função solve(Σ) como padrão. Desisti de usar a matriz de Mahalanobis, pensei em distancia Euclidiana, mas seria importante eu conseguir usar um método que envolvesse a matriz de covariâncias para obter D2. Minha pergunta é se teria uma forma de usar, neste caso, um outro tipo de inversa (outro algorítimo) e sem comprometer a matriz D2 ou se vou mesmo ter ficar com distancia Euclidiana e variantes.
setwd("/home/andre")
dados <- read.csv("dados.txt", sep="")
str(dados)
'data.frame': 6 obs. of 10 variables: $ ch : int 146 26 44 56 73 101 $ lam : int 126 14 20 104 54 156 $ esfr : int 208 1 1 36 73 22 $ cut : int 43 10 9 20 22 40 $ fleh : int 16 1 1 2 9 5 $ tmont : int 45 5 12 15 11 27 $ mont : int 18 8 16 25 21 32 $ montco : int 6 0 3 2 2 3 $ vocal : int 106 14 21 105 97 151 $ respaud: int 123 44 64 113 56 79
dados ch lam esfr cut fleh tmont mont montco vocal respaudA 146 126 208 43 16 45 18 6 106 123B 26 14 1 10 1 5 8 0 14 44C 44 20 1 9 1 12 16 3 21 64D 56 104 36 20 2 15 25 2 105 113E 73 54 73 22 9 11 21 2 97 56F 101 156 22 40 5 27 32 3 151 79
require(StatMatch)D2=mahalanobis.dist(dados)mahalanobis.dist(dados)Error in solve.default(cov, ...) : sistema é computacionalmente singular: condição recíproca número = 1.53112e-19
require(biotools)Cov=cov(dados)
D2=D2.dist(dados, Cov) D2=D2.dist(dados, Cov) Error in solve.default(cov, (x[i, ] - x[j, ])) : sistema é computacionalmente singular: condição recíproca número = 1.53112e-19
André Oliveira Souza. Graduação em Matemática, mestrado em estatística aplicada.Instituto Federal de Educação, Ciência e Tecnologia do Espirito Santo. IFES
Em Quinta-feira, 25 de Junho de 2015 19:41, walmes . <walmeszeviani em gmail.com> escreveu: Seu exemplo é 100% não reproduzível.À disposição.Walmes.
_______________________________________________
R-br mailing list
R-br em listas.c3sl.ufpr.br
https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20150626/f568646c/attachment.html>
Mais detalhes sobre a lista de discussão R-br