#==========================================================================================
# Aula 16 da disciplina ce223 (13/05/2011)
# Estatística básica, construção de funções com opções de controle
#                                                               Professor Walmes M. Zeviani
#                                                                     www.leg.ufpr.br/ce223
#==========================================================================================

#------------------------------------------------------------------------------------------
# aprimorando a função para teste de hipótese sobre a média incluindo o tipo de teste
# bilateral, unilateral à esquerda, unilateral à direita

my.test <- function(x, y, alpha, tipo){
  ## x, é o vetor com dados
  ## y, é o vetor ou escalar
  ## alpha, é o nível de significância do teste
  ## é o tipo de teste: bilateral, à esquerda, à direita
  if(length(y)==1){ # vai fazer o teste para uma média
    h0 <- y
    m <- mean(x); v <- var(x); n <- length(x); e <- sqrt(v/n)
    tc <- (m-h0)/e
    if(tipo=="bilateral"){ # vai aplicar o teste bilateral, P(xbar!=h0)
      tt <- -pt(alpha/2, df=n-1)
      dec <- ifelse(abs(tc)>tt, "rejeito", "aceito")
    }
    if(tipo=="esquerda"){ # vai aplicar o teste unilateral à esquerda, P(xbar<h0)
      tt <- pt(alpha, df=n-1)
      dec <- ifelse(tc>tt, "rejeito", "aceito")
    }
    if(tipo=="direita"){ # vai aplicar o teste unilateral à direita, P(xbar>h0)
      tt <- -pt(alpha, df=n-1)
      dec <- ifelse(tc>tt, "rejeito", "aceito")
    }
    return(list(stat=c(media=m, variancia=v, erropad=e, tcalc=tc, ttab=tt),
                n=n, decisao=dec, tipo=tipo))
  }
  if(length(y)>1){ # vai fazer o teste para diferença de duas médias
    a1 <- x; a2 <- y; nivel <- alpha
    m1 <- mean(a1); m2 <- mean(a2)
    n1 <- length(a1); n2 <- length(a2)
    v1 <- var(a1); v2 <- var(a2)
    v <- ((n1-1)*v1+(n2-1)*v2)/(n1+n2-2)
    dif <- abs(m2-m1)
    tt <- -qt((1-nivel)/2, df=n1+n2-2)
    ep <- sqrt(v*(1/n1+1/n2))
    tc <- (dif-0)/ep
    return(list(medias=c(m1=m1, m2=m2, dif=dif),
                variancias=c(v1=v1, v2=v2, v=v),
                n=c(n1=n1, n2=n2, gl=n1+n2-2),
                statistica=c(tc=tc, tt=tt)))
  }
}

#------------------------------------------------------------------------------------------

y <- rbeta(100, 2, 3)
z <- rbeta(80, 2.5, 3)
my.test(y, h0=0.5, alpha=0.05, tipo="bilateral")
my.test(y, h0=0.5, alpha=0.05, tipo="esquerda")
my.test(y, h0=0.5, alpha=0.05, tipo="direita")

my.test(y, z, alpha=0.05, tipo="direita")
my.test(y, 0.5, alpha=0.05, tipo="direita")

#------------------------------------------------------------------------------------------

help(t.test, help_type="html")

#------------------------------------------------------------------------------------------

x1 <- rgamma(20, 1, 0.9)
x2 <- rgamma(22, 1.1, 0.9)

t.test(x1, x2, var.equal=TRUE)

dados <- data.frame(amostra=rep(c(1,2), c(length(x1), length(x2))),
                    valor=c(x1, x2))
dados

dados <- dados[sample(1:nrow(dados)),]
dados

subset(dados, amostra==1)

with(dados, t.test(valor~amostra, var=TRUE))

#------------------------------------------------------------------------------------------

dados <- data.frame(peso=rnorm(100, 80, 4),
                    altura=rnorm(100, 1.78, 0.2),
                    idade=rpois(100, lambda=22),
                    renda=rexp(100, rate=0.01))
str(dados)

aux <- apply(dados, 2, t.test)
str(aux)

aux$peso$conf

lapply(aux,
       function(slot){
         c(slot$conf)
       })

sapply(aux,
       function(slot){
         c(slot$conf)
       })

#------------------------------------------------------------------------------------------