Mestranda em Ecologia, Instituto de Biologia, USP.

Março de 2009

Uma função para análise exploratória de dados que sejam dois vetores através de apresentações gráficas, resumos e opções de simulação de uma distribuição nula para teste unicaudal, bicaudal e teste t. Apresenta gráficos para exploração das variáveis, como histograma, boxplot, plot, qqnorm, e também retorna um gráfico de simulação de uma distribuição nula para teste unicaudal, bicaudal e teste t. É uma função integrando a análise exploratória de dados através de apresentações gráficas e opções de simulação usadas na função “simula”.

Idéia promissora, mas algumas coisas não ficaram claras para mim, veja se estão para você:

1. boxplot e qqnorm são gráficos para uma variável e sua função receberá duas. Um gráfico para cada uma?
2. Que simulação será feita e que gráfico será retornado? Há vários pontos importantes aqui, a começar pela pergunta que se pretende responder com a simulação, o que vai determinar o tipo de simulação feita, o melhor modo de mostrar os resultados, etc.

Uma função para planilha de dados (data frame) organizados com cada indivíduo da amostra nas linhas e com uma coluna com o nome das espécies correspondente a cada indivíduo.

A função retorna uma matriz e gráfico do rank de abundância das espécies, e terá argumento de opção para aplicar a função para apenas algum ou vários subconjuntos da coleta dos dados (como locais, armadilhas, meses ou parcelas) os quais você queira saber essas informações (abundância de cada espécie).

E quem sabe também apresente cálculos de medidas ecológicas (como índices de diversidade e equitabilidade).

Bem geral e factível. Se entendi sua planilha tem uma linha para cada indivíduo, obrigatoriamente um fator que identifica a espécie, e outro fatores opcionais que podem ser usados para separar os dados em subconjuntos, é isto? Uma dica: resista à tentação de usar a função sort, e opte pela order, que parece mais difícil, mas é bem mais poderosa. Só ela pode resolver problemas de ordenação mais complicados como este. Outra sugestão: você pode oferecer a opção da planilha ter uma linha por espécie, e um campo de abundância de cada espécie.

analix                package:nenhum                R Documentation



Análise exploratória de duas variáveis e opção para simulação nula com teste 
unicaudal, bicaudal e teste t.


Description:

Apresenta sumário das variáveis e os principais gráficos para uma análise exploratória
de dados, índice de correlação entre duas variáveis, e retorna uma simulação de uma 
distribuição nula com testes estatísticos.  


Usage:

     analix <- function(dados1,dados2,teste="N")


Arguments:

 dados1: Vetor numérico. Valores de uma amostra ou variável. 

 dados2: numérico. Valores de uma amostra ou variável.
 
 teste: lógico. Qual o tipo de teste na simulação nula?
 

Details: 

Os valores das amostras (dados1,dados2) são plotados em gráficos de 
histograma, boxplot, qqnorm e "dados2" em função de "dados1", com linha 
da regressão linear simples.
É apresentado um sumário de cada variável (ou amostra) com o índice de 
correlação entre elas.
Também são simulados uma hipótese nula para os dados com opção de teste 
unicaudal (teste="uni"), teste bicaudal (teste="bi"), teste t de Student 
(teste="t") ou sem simulação e teste estatístico (teste="N"). 
As simulações são feitas com repetição de 1000 de uma distribuição
uniforme dos dados. 

 

Value:

 São gerados um conjunto de seis gráficos para análise exploratória dos 
dados, e um gráfico com o resultado da simulação de hipótese nula com o 
teste estatístico escolhido.
 

Author(s):

Paula Alves Condé

paula.conde@usp.br


References:
 
Exploratory data analysis. http://en.wikipedia.org/wiki/Exploratory_data_analysis

Null hypothesis. http://en.wikipedia.org/wiki/Null_hypothesis

Statistical hypothesis testing. http://en.wikipedia.org/wiki/Hypothesis_testing


See Also:

'qqnorm' para analisar a normalidade dos dados.
'simula' simulação da hipótese nula com testes estatísticos. 


Examples:

  dados1=round(rnorm(10,mean=6,sd=3),2)
  dados2=round(rnorm(10,mean=7.5,sd=3.2),2)

## Sem simulações e teste, apenas exploração gráfica das variáveis:	      
  analix(dados1,dados2,teste="N")  

## Exploração gráfica + Simulação de uma distribuição nula com teste bicaudal:   
  analix(dados1,dados2,teste="bi")  

## Exploração gráfica + Simulação de uma distribuição nula com teste unicaudal:
  analix(dados1,dados2,teste="uni")  

## Exploração gráfica + Simulação de uma distribuição nula com teste t:
  analix(dados1,dados2,teste="t")  

analix <- function(dados1,dados2,teste="N")
{       
    ### Função utilizada para analise exploratória de duas variáveis.
    ### A entrada de dados deve ser feita por dois objetos vetores de mesmo tamanho.
    ### Também faz simulação de uma distribuição nula para teste unicaudal, bicaudal 
## e teste t.

cat("Para fazer simulação de uma distribuição nula, no caso de teste unicaudal o 
+ primeiro vetor deve ser o dos dados que seriam maiores.\n As opçoes de teste são 
+ (entre aspas):\n\t uni (normal unicaudal)\n\t bi (normal bicaudal) \n\t 
+ t (distribuição t)\n")
       

 ### Autora: Paula Alves Condé
 ### Data: 01/04/2009
 ### Trabalho final 
 ### Disciplina "Uso da Linguagem R para Analise de Dados Ecologicos" - BIE5782.
 
        x11()
        par(mfrow=c(3,2),bty="l")
        
          ## Produz Histograma dos "dados1"
       hist(dados1,col="red",main=paste(c("Histograma",substitute(dados1))))
          ## Produz Histograma dos "dados2"
       hist(dados2,col="blue",main=paste(c("Histograma",substitute(dados2))))
         
          ## Produz Boxplot das duas variáveis
       boxplot(dados1,dados2,main="Boxplot",names=c(substitute(dados1),
 + substitute(dados2)))
         
          ## Plota "dados2" em função dos "dados1" com a reta de regressão 
 ## linear simples para ver a relação entre as variáveis
         plot(dados2~dados1,main=paste(c(substitute(dados2),"~",substitute(dados1)))
+ ,pch=16,xlab=deparse(substitute(dados1)),ylab=deparse(substitute(dados2)))
         abline(lm(dados2~dados1),col="purple",lwd=2)
         
            ## Gráfico QQNorm para "Dados1" 
         qqnorm(dados1)
         qqline(dados1,col="red")
         mtext(deparse(substitute(dados1)),side=3)  
            
            ## Gráfico QQNorm para "Dados2"
         qqnorm(dados2)
         qqline(dados2,col="blue")
         mtext(deparse(substitute(dados2)),side=3)
                
                 ## Sumário das variáveis e índice de correlação
          name.xy <- paste(deparse(substitute(dados1)), "e", deparse(substitute
+ (dados2)))
	    cat("Dados:", name.xy, "\n")
          n <- length(dados1)
          sumx <- sum(dados1^2) - sum(dados1)^2/n
	    sumy <- sum(dados2^2) - sum(dados2)^2/n
	    sumxy <- sum(dados1 * dados2) - sum(dados1) * sum(dados2)/n
	     correl <- sumxy/sqrt(sumx * sumy)
	       concluindo<- list(summary(dados1),summary(dados2),correl)
	       names(concluindo)<- c(substitute(dados1), substitute(dados2),
+ "correlação")
 

  ### Fazendo simulação de uma distribuição nula para teste unicaudal, bicaudal
## ou teste t.  
    
      nsim=1000         
         resultado<-rep(NA,nsim)
         dif = mean(dados1) - mean(dados2)
         dif.abs = round(abs(dif), 1) 
    cat("\n Diferença absoluta observada entre as médias das variáveis = ",dif.abs,
+ "\n")
    v1 = var(dados1)
    v2 = var(dados2)
    n1 = length(dados1)
    n2 = length(dados2)
    s12 = sqrt((v1/n1) + (v2/n2))
    tvalor = abs(dif/s12)
    med = mean(c(dados1, dados2))
    des = sd(c(dados1, dados2))
    res = rep(NA,nsim)
    arco = rainbow(nsim)
x11()

  if (teste == "bi") 
        {
        meu.cex = 900/nsim
        plot(runif(50, 0, (dif.abs * 1.5)), 0:49, type = "n", 
            xlim = c(0, dif.abs * 1.5), ylim = c(0, (0.08 * nsim)), 
            xlab = "Diferença absoluta", ylab = "Frequência", 
            main = "Simulação")
        for (i in 1:nsim) {
            res[i] = abs(round(mean(rnorm(n1, mean = med, sd = des)) - 
                mean(rnorm(n2, mean = med, sd = des)), 1))
            stripchart(res, method = "stack", add = T, cex = meu.cex, 
                pch = 15, col = arco[i])
         }
        res.menor <- res[res < dif.abs]
        stripchart(res.menor, method = "stack", add = T, cex = meu.cex, 
            pch = 15, col = "dark green")
        legend(5, (6 * nsim), legend = paste(sum(res >= dif.abs), 
            "valores >= ", round(dif.abs, 1)), bty = "n", text.col = "red")
        abline(v = dif.abs, col = "red")
         }
    if (teste == "uni") 
        {
        meu.cex = 1700/nsim
        plot(runif(50, (dif.abs * -1.5), (dif.abs * 1.5)), 0:49, 
            type = "n", xlim = c((dif * -1.5), (dif * 1.5)), 
            ylim = c(0, (0.04 * nsim)), xlab = " Diferença Absoluta", 
            ylab = "Frequência", main = "Simulação")
        for (i in 1:nsim) {
            res[i] = round(mean(rnorm(n1, mean = med, sd = des)) - 
                mean(rnorm(n2, mean = med, sd = des)), 1)
            stripchart(res, method = "stack", add = T, cex = meu.cex, 
                pch = 15, col = arco[i])
        }
        res.menor <- res[res < dif.abs & res > (-1 * dif.abs)]
        stripchart(res.menor, method = "stack", add = T, cex = meu.cex, 
            pch = 15, col = "dark green")
        legend((dif.abs * 0.5), (0.038 * nsim), legend = paste(sum(res >= 
            dif.abs), "valores >= ", round(dif, 1)), bty = "n", 
            text.col = "red")
        legend((dif.abs * -1.4), (0.038 * nsim), legend = paste(sum(res <= 
            (dif.abs * -1)), "valores <= -", round(dif, 1)), 
            bty = "n", text.col = "red")
        abline(v = dif.abs, col = "red")
        abline(v = (-1 * dif.abs), col = "red")
        }
    
   if (teste == "t") 
        {
        meu.cex = 1300/nsim
        cat("\t\nValor t observado = ", tvalor, "\n")
        plot(runif(50, (-0.7 * dif.abs), (0.7 * dif.abs)), runif(50, 
            0, (nsim/20)), type = "n", xlim = c((-0.8 * dif.abs), 
            (0.8 * dif.abs)), ylim = c(0, (nsim/20)), xlab = " valor t ", 
            ylab = "Frequência", main = "Simulação")
        for (i in 1:nsim) {
            simula1 = rnorm(n1, mean = med, sd = des)
            simula2 = rnorm(n2, mean = med, sd = des)
            difs = mean(simula1) - mean(simula2)
            vs1 = var(simula1)
            vs2 = var(simula2)
            ss12 = sqrt((vs1/n1) + (vs2/n2))
            tsimula = difs/ss12
            res[i] = round(tsimula, 1)
            stripchart(res, method = "stack", add = T, cex = meu.cex, 
                pch = 15, col = arco[i])
        }
        tvalor.vf = res < round(tvalor, 1) & res > round((-1 * 
            tvalor), 1)
        ntvalor = sum(tvalor.vf == F)
        res.menor = res[tvalor.vf]
        stripchart(res.menor, method = "stack", add = T, cex = meu.cex, 
            pch = 15, col = "dark green")
        legend((tvalor * 0.2), (nsim/22), legend = paste(sum(res >= 
            round(tvalor, 1)), "valores >= ", round(tvalor, 1)), 
            bty = "n", text.col = "red")
        legend((tvalor * -0.9), (nsim/22), legend = paste(sum(res <= 
            round((tvalor * -1), 1)), "valores <= -", round(tvalor, 
            1)), bty = "n", text.col = "red")
        abline(v = tvalor, col = "red")
        abline(v = (-1 * tvalor), col = "red")
        cat("\n\tProbabilidade erro I = ", ntvalor/nsim, "\n")
       }
       return(concluindo)
}

analix.r