Re: [R-br] [R_STAT] Teste de esfericidade
Michele, Acho que você estava se referindo a função sphericiy.test que eu implementei e fiz referẽncia em um material sobre ANOVA com medidas repetidas. Deixei o código disponível aqui: http://www.fernandohrosa.com.br/en/P/sphericity-test-for-covariance-matrices... Fernando, 2008/2/25 michele barbosa <mixxele@yahoo.com.br>
**
Como encontrar no help conteúdo para teste de esfericidade? Algo referente ao comando: *sphericity.test* ** obrigada Michele
------------------------------ Abra sua conta no Yahoo! Mail<http://br.rd.yahoo.com/mail/taglines/mail/*http://br.mail.yahoo.com/>, o único sem limite de espaço para armazenamento!
__._,_.___ Mensagens neste tópico <http://br.groups.yahoo.com/group/R_STAT/message/2683;_ylc=X3oDMTM1ZG1uYWw5BF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRtc2dJZAMyNjgzBHNlYwNmdHIEc2xrA3Z0cGMEc3RpbWUDMTIwMzk3ODE1OQR0cGNJZAMyNjgz>( 1) Responder (através da web) <http://br.groups.yahoo.com/group/R_STAT/post;_ylc=X3oDMTJxdDB1bmdlBF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRtc2dJZAMyNjgzBHNlYwNmdHIEc2xrA3JwbHkEc3RpbWUDMTIwMzk3ODE1OQ--?act=reply&messageNum=2683>| Adicionar um novo tópico <http://br.groups.yahoo.com/group/R_STAT/post;_ylc=X3oDMTJmOXFwdmd2BF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNudHBjBHN0aW1lAzEyMDM5NzgxNTk-> Mensagens<http://br.groups.yahoo.com/group/R_STAT/messages;_ylc=X3oDMTJmNW9qdjR2BF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNtc2dzBHN0aW1lAzEyMDM5NzgxNTk->| Arquivos<http://br.groups.yahoo.com/group/R_STAT/files;_ylc=X3oDMTJnaG5zYzJzBF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNmaWxlcwRzdGltZQMxMjAzOTc4MTU5>| Fotos<http://br.groups.yahoo.com/group/R_STAT/photos;_ylc=X3oDMTJmbWRzZ3BjBF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNwaG90BHN0aW1lAzEyMDM5NzgxNTk->| Links<http://br.groups.yahoo.com/group/R_STAT/links;_ylc=X3oDMTJnamVkZWw1BF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNsaW5rcwRzdGltZQMxMjAzOTc4MTU5>| Banco de dados<http://br.groups.yahoo.com/group/R_STAT/database;_ylc=X3oDMTJkamg0ZWVqBF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNkYgRzdGltZQMxMjAzOTc4MTU5>| Enquetes<http://br.groups.yahoo.com/group/R_STAT/polls;_ylc=X3oDMTJncjBzczVxBF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNwb2xscwRzdGltZQMxMjAzOTc4MTU5> [image: Yahoo! Grupos]<http://br.groups.yahoo.com/;_ylc=X3oDMTJlZmFiNHFxBF9TAzk3NDkwNDM1BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNnZnAEc3RpbWUDMTIwMzk3ODE1OQ--> Alterar configurações via web<http://br.groups.yahoo.com/group/R_STAT/join;_ylc=X3oDMTJnaGRiMnBvBF9TAzk3NDkwNDM1BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNzdG5ncwRzdGltZQMxMjAzOTc4MTU5>(Requer Yahoo! ID) Alterar configurações via e-mail: Alterar recebimento para lista diária de mensagens<R_STAT-digest@yahoogrupos.com.br?subject=+Recebimento+de+e-mail:+Lista+de+mensagens>| Alterar formato para o tradicional<R_STAT-traditional@yahoogrupos.com.br?subject=Alterar+formato+de+distribui%C3%A7%C3%A3o:+Tradicional> Visite seu Grupo <http://br.groups.yahoo.com/group/R_STAT;_ylc=X3oDMTJldjFsajNkBF9TAzk3NDkwNDM1BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDZnRyBHNsawNocGYEc3RpbWUDMTIwMzk3ODE1OQ-->| Termos de uso do Yahoo! Grupos <http://br.yahoo.com/info/utos.html> | Sair do grupo <R_STAT-unsubscribe@yahoogrupos.com.br?subject=> Atividade nos últimos dias
Visite seu Grupo <http://br.groups.yahoo.com/group/R_STAT;_ylc=X3oDMTJmanAxcGFzBF9TAzk3NDkwNDM3BGdycElkAzExOTIzMjc1BGdycHNwSWQDMjEzNzExMTYwNQRzZWMDdnRsBHNsawN2Z2hwBHN0aW1lAzEyMDM5NzgxNTk-> Yahoo! Mail
Proteção anti-spam
Muito mais espaço Yahoo! Barra
Buscar sites na web
Checar seus e-mails . Yahoo! Grupos
Crie seu próprio grupo<http://br.groups.yahoo.com/;_ylc=X3oDMTJwcWI3bGZvBF9TAzk3NDkwNDM3BF9wAzMEZ3JwSWQDMTE5MjMyNzUEZ3Jwc3BJZAMyMTM3MTExNjA1BHNlYwNuY21vZARzbGsDZ3JvdXBzMgRzdGltZQMxMjAzOTc4MTU5>
A melhor forma de comunicação
.
__,_._,___
-- "Though this be randomness, yet there is structure in't." Fernando H Rosa - Statistician http://www.fernandohrosa.com.br / http://www.feferraz.net - Estatística, Matemática e Computação BankReview.com.br <http://www.bankreview.com.br/> - Escolha melhor seus serviços financeiros! AprendaAlemao.net <http://aprendaalemao.net/> - Seu ponto de partida para melhorar seu Alemão! @fhrosa
Estava olhando a discussão, tentei rodar a função e não consegui, como faço?
sfericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){
+ #### Performs a hypothesis test that a covariance matrix is of specified
+ #### form. Test is of the form H0: S1=sigma^2*S2. n is the number of
+ #### observations on which the sample covariance matrix is based.
+ #### If the input parameter estsigma is TRUE:
+ #### Perform test of the hypothesis that S1=sigma^2 S2, for unknown sigma.
+ #### If S2 not specified, assumed that S2=I. Reference is Basilevsky,
+ #### Statistical Factor Analysis and Related Methods, page 191.
+ #### If the input parameter estsigma is FALSE:
+ #### Perform test of the hypothesis that S1=S2. If S2 not specified,
+ #### assumed that S2=I. Reference is Seber, Multivariate Observations,
+ #### sec 3.5.4
+ #### Only the lower triangle+diagonal is required at entry, and the upper
+ #### triangle is ignored.
+ #### DAW July 2000
+ dname <- paste(substitute(s1))
+ p<-nrow(s1)
+ for (i in 1:(p-1)){for (j in ((i+1):p)){
+ s1[i,j]<-s1[j,i]
+ s2[i,j]<-s2[j,i] }}
+ if (!is.null(s2)){
+ b<-eigen(s2,symmetric=T,only.values=F)
+ r<-b$vectors %*% diag(1/sqrt(b$values))
+ s<-t(r) %*% s1 %*% r }
+ else { s<-s1 }
+
+ d<-eigen(s,symmetric=T,only.values=T)$values
+ ldet<-sum(log(d))
+ tr<-sum(d)
+
+ if (estsigma==TRUE){
+ sighat<-tr/p
+ cc<--(n-(2*p^2+p+2)/(6*p))*(ldet-p*log(tr/p))
+ statistic <- cc
+ sighat<-sighat
+ names(statistic) <- "L statistic"
+ parameter <- 0.5*(p+2)*(p-1)
+ names(parameter) <- "df"
+ rval<-list(data.name=dname,sighat=sighat,statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Sphericity test") }
+ else {
+ cc<--n*(p+ldet-tr)
+ statistic <- cc
+ names(statistic) <- "L statistic"
+ parameter <- 0.5*(p+1)*p
+ names(parameter) <- "df"
+ rval<-list(data.name="",statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Covariance equality test statistic")
+ }
+ class(rval) <- "htest"
+ return(rval)
+ }
> pw <- function(q,n) {
+ pdf <- function(w) { 1/2 * (n-2) * w^((n-3)/2) }
+ integrate(pdf,0,q)
+ }
>
> varcomp <- function(covmat,n) {
+ if (is.list(covmat)) {
+ if (length(covmat) < 2)
+ stop("covmat must be a list with at least 2 elements")
+ ps <- as.vector(sapply(covmat,dim))
+ if (sum(ps[1] == ps) != length(ps))
+ stop("all covariance matrices must have the same dimension")
+ p <- ps[1]
+ q <- length(covmat)
+ if (length(n) == 1)
+ Ng <- rep(n,q)
+ else if (length(n) == q)
+ Ng <- n
+ else
+ stop("n must be equal length(covmat) or 1")
+
+ DNAME <- deparse(substitute(covmat))
+ }
+
+ else
+ stop("covmat must be a list")
+
+ ng <- Ng - 1
+ Ag <- lapply(1:length(covmat),function(i,mat,n) { n[i] * mat[[i]] },mat=covmat,n=ng)
+ A <- matrix(colSums(matrix(unlist(Ag),ncol=p^2,byrow=T)),ncol=p)
+ detAg <- sapply(Ag,det)
+ detA <- det(A)
+ V1 <- prod(detAg^(ng/2))/(detA^(sum(ng)/2))
+ kg <- ng/sum(ng)
+ l1 <- prod((1/kg)^kg)^(p*sum(ng)/2) * V1
+ rho <- 1 - (sum(1/ng) - 1/sum(ng))*(2*p^2+3*p-1)/(6*(p+1)*(q-1))
+ w2 <- p*(p+1) * ((p-1)*(p+2) * (sum(1/ng^2) - 1/(sum(ng)^2)) - 6*(q-1)*(1-rho)^2) / (48*rho^2)
+ f <- 0.5 * (q-1)*p*(p+1)
+ STATISTIC <- -2*rho*log(l1)
+ PVAL <- 1 - (pchisq(STATISTIC,f) + w2*(pchisq(STATISTIC,f+4) - pchisq(STATISTIC,f)))
+ names(STATISTIC) <- "corrected lambda*"
+ names(f) <- "df"
+ RVAL <- structure(list(statistic = STATISTIC, parameter = f,p.value = PVAL, data.name = DNAME, method = "Equality of Covariances Matrices Test"),class="htest")
+ return(RVAL)
+ }
> sphericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){
+ >
> sphericity.test(10,2)
Erro em 1:(p - 1) : argumento de comprimento zero
Michele,
Acho que você estava se referindo a função sphericiy.test que eu implementei e fiz referẽncia em um material sobre ANOVA com medidas repetidas.
Deixei o código disponível aqui: http://www.fernandohrosa.com.br/en/P/sphericity-test-for-covariance-matrices-in-r-sphericity-test/
Fernando,
.
Você chamou a função com dois escalares:
sphericity.test(10,2)
O segumento argumento tem que ser uma matriz de covariância, e não um
escalar (2, no caso).
2011/10/9 Mauro Sznelwar <sznelwar@uol.com.br>
> **
> *Estava olhando a discussão, tentei rodar a função e não consegui, como
> faço?*
>
> sfericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){
> + #### Performs a hypothesis test that a covariance matrix is of specified
> + #### form. Test is of the form H0: S1=sigma^2*S2. n is the number of
> + #### observations on which the sample covariance matrix is based.
> + #### If the input parameter estsigma is TRUE:
> + #### Perform test of the hypothesis that S1=sigma^2 S2, for unknown
> sigma.
> + #### If S2 not specified, assumed that S2=I. Reference is Basilevsky,
> + #### Statistical Factor Analysis and Related Methods, page 191.
> + #### If the input parameter estsigma is FALSE:
> + #### Perform test of the hypothesis that S1=S2. If S2 not specified,
> + #### assumed that S2=I. Reference is Seber, Multivariate Observations,
> + #### sec 3.5.4
> + #### Only the lower triangle+diagonal is required at entry, and the upper
>
> + #### triangle is ignored.
> + #### DAW July 2000
> + dname <- paste(substitute(s1))
> + p<-nrow(s1)
> + for (i in 1:(p-1)){for (j in ((i+1):p)){
> + s1[i,j]<-s1[j,i]
> + s2[i,j]<-s2[j,i] }}
> + if (!is.null(s2)){
> + b<-eigen(s2,symmetric=T,only.values=F)
> + r<-b$vectors %*% diag(1/sqrt(b$values))
> + s<-t(r) %*% s1 %*% r }
> + else { s<-s1 }
> +
> + d<-eigen(s,symmetric=T,only.values=T)$values
> + ldet<-sum(log(d))
> + tr<-sum(d)
> +
> + if (estsigma==TRUE){
> + sighat<-tr/p
> + cc<--(n-(2*p^2+p+2)/(6*p))*(ldet-p*log(tr/p))
> + statistic <- cc
> + sighat<-sighat
> + names(statistic) <- "L statistic"
> + parameter <- 0.5*(p+2)*(p-1)
> + names(parameter) <- "df"
> + rval<-list(data.name=dname,sighat=sighat,statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Sphericity
> test") }
> + else {
> + cc<--n*(p+ldet-tr)
> + statistic <- cc
> + names(statistic) <- "L statistic"
> + parameter <- 0.5*(p+1)*p
> + names(parameter) <- "df"
> + rval<-list(data.name="",statistic=statistic,parameter=parameter,p.value=1-pchisq(statistic,parameter),method="Covariance
> equality test statistic")
> + }
> + class(rval) <- "htest"
> + return(rval)
> + }
> > pw <- function(q,n) {
> + pdf <- function(w) { 1/2 * (n-2) * w^((n-3)/2) }
> + integrate(pdf,0,q)
> + }
> >
> > varcomp <- function(covmat,n) {
> + if (is.list(covmat)) {
> + if (length(covmat) < 2)
> + stop("covmat must be a list with at least 2 elements")
> + ps <- as.vector(sapply(covmat,dim))
> + if (sum(ps[1] == ps) != length(ps))
> + stop("all covariance matrices must have the same dimension")
> + p <- ps[1]
> + q <- length(covmat)
> + if (length(n) == 1)
> + Ng <- rep(n,q)
> + else if (length(n) == q)
> + Ng <- n
> + else
> + stop("n must be equal length(covmat) or 1")
> +
> + DNAME <- deparse(substitute(covmat))
> + }
> +
> + else
> + stop("covmat must be a list")
> +
> + ng <- Ng - 1
> + Ag <- lapply(1:length(covmat),function(i,mat,n) { n[i] * mat[[i]]
> },mat=covmat,n=ng)
> + A <- matrix(colSums(matrix(unlist(Ag),ncol=p^2,byrow=T)),ncol=p)
> + detAg <- sapply(Ag,det)
> + detA <- det(A)
> + V1 <- prod(detAg^(ng/2))/(detA^(sum(ng)/2))
> + kg <- ng/sum(ng)
> + l1 <- prod((1/kg)^kg)^(p*sum(ng)/2) * V1
> + rho <- 1 - (sum(1/ng) - 1/sum(ng))*(2*p^2+3*p-1)/(6*(p+1)*(q-1))
> + w2 <- p*(p+1) * ((p-1)*(p+2) * (sum(1/ng^2) - 1/(sum(ng)^2)) -
> 6*(q-1)*(1-rho)^2) / (48*rho^2)
> + f <- 0.5 * (q-1)*p*(p+1)
> + STATISTIC <- -2*rho*log(l1)
> + PVAL <- 1 - (pchisq(STATISTIC,f) + w2*(pchisq(STATISTIC,f+4) -
> pchisq(STATISTIC,f)))
> + names(STATISTIC) <- "corrected lambda*"
> + names(f) <- "df"
> + RVAL <- structure(list(statistic = STATISTIC, parameter = f,p.value =
> PVAL, data.name = DNAME, method = "Equality of Covariances Matrices
> Test"),class="htest")
> + return(RVAL)
> + }
> > sphericity.test<-function(n,s1,s2=NULL,estsigma=TRUE){
>
> + >
> > sphericity.test(10,2)
> Erro em 1:(p - 1) : argumento de comprimento zero
>
>
> Michele,
>
> Acho que você estava se referindo a função sphericiy.test que eu
> implementei e fiz referẽncia em um material sobre ANOVA com medidas
> repetidas.
>
> Deixei o código disponível aqui:
> http://www.fernandohrosa.com.br/en/P/sphericity-test-for-covariance-matrices-in-r-sphericity-test/
>
> Fernando,
>
> .
>
>
> _______________________________________________
> R-br mailing list
> R-br@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.
>
--
"Though this be randomness, yet there is structure in't."
Fernando H Rosa - Statistician
http://www.fernandohrosa.com.br / http://www.feferraz.net - Estatística,
Matemática e Computação
BankReview.com.br <http://www.bankreview.com.br/> - Escolha melhor seus
serviços financeiros!
AprendaAlemao.net <http://aprendaalemao.net/> - Seu ponto de partida para
melhorar seu Alemão!
@fhrosa
participantes (2)
-
Fernando Henrique Ferraz Pereira da Rosa -
Mauro Sznelwar