################################################
## Funcao para simular o teste t por permutacao
## Alexandre Adalardo de Oliveira
## 08 jun 2022 adaptado da antiga simula.r
################################################
simulaT = function(dados1, dados2, nsim = 1000, teste = c("menor", "maior", "difere"), anima = TRUE)
{
teste <- match.arg(teste)
dados1 <- dados1[!is.na(dados1)]
dados2 <- dados2[!is.na(dados2)]
dif <- mean(dados1) - mean(dados2)
difAbs= abs(dif)
#cat("\n Diferenca absoluta observada entre as medias das variaveis = ", dif.abs, "\nPara salvar os valores gerados, atribua a funcao a um objeto/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)
#################################################
med1sim <- apply(matrix(rnorm((n1* nsim) , mean = med, sd = des), nrow = nsim, ncol = n1), 1, mean)
med2sim <- apply(matrix(rnorm((n2* nsim) , mean = med, sd = des), nrow = nsim, ncol = n2), 1, mean)
res <- c( med1sim - med2sim)
res[1] <- dif
if (teste=="maior")
{
    pvalor<- sum(res>=dif)/nsim
}
if (teste == "menor")
{
    pvalor <- sum(res<=dif)/nsim
}
if (teste=="difere")
{
    pvalor =  sum(abs(res) >= abs(dif))/nsim
}
if(anima)
{
    par(mar = c(5, 5, 2, 2), cex.lab = 1.5, cex.axis = 1.2, las = 1) 
    histSimT <- hist(res, plot = FALSE)
    if(nsim >= 5000)
    {
        seqHist <-  seq(10, nsim, by = 5)
    }else{
        seqHist <- 1:nsim
    }
    cores <- rev(topo.colors(length(seqHist), alpha = .5)) 
    for(i in 1:length(seqHist))
    {
         hist(res[1:seqHist[i]], col = cores[i], border=NA, breaks = histSimT$breaks, ylim= c(0, max(histSimT$counts)), main="", ylab="Frequência", xlab="Diferença entre as médias") 
    }
    hist(res, col = rgb(0, 0, 1, 0.4), border=NA, breaks = histSimT$breaks, ylim= c(0, max(histSimT$counts)), main=paste("Diferenças entre médias (", teste,", pvalor = ", round(pvalor,4), ")", sep = ""), ylab="Frequência", xlab="Diferença entre as médias", bg = rgb(0,0,0,,15) )
    rect(par()$usr[1], par()$usr[3], par()$usr[2], par()$usr[4], col = rgb(0, 0, 0, 0.15))
    if (teste=="maior")
    {
        abline(v = dif, lty= 2 , lwd = 3, col = "white")
        vfMaior <- histSimT$breaks >= dif
        poly <- c(rep(histSimT$counts[which(vfMaior) - 1], each = 2), 0,0, histSimT$counts[ (which(vfMaior) -1)[1] ])
        polx <- c(dif, rep(histSimT$breaks[vfMaior], each = 2), dif, dif)
        polygon(x = polx, y = poly, col = rgb(1,0,0, 0.7), border = NA)
    }
    if(teste == "menor")
    {
        abline(v = dif, lty= 2 , lwd = 3, col = "white")
        vfMenor <- histSimT$breaks <= dif
        poly <- c(0, histSimT$counts[rep(which(vfMenor), each = 2)], 0, 0)
        polx <- c(rep(histSimT$breaks[vfMenor], each = 2), dif, dif, histSimT$breaks[1])
        polygon(x = polx, y = poly, col = rgb(1,0,0, 0.7), border = NA)
    }
    if(teste == "difere")
    {
        abline(v = c(difAbs, -1* difAbs), lty= 2 , lwd = 3, col = "white")
        vfMaiorU <- histSimT$breaks >= difAbs
        vfMenorL <- histSimT$breaks <= (-1 * difAbs)
        polyU <- c(rep(histSimT$counts[which(vfMaiorU) - 1], each = 2), 0,0, histSimT$counts[ (which(vfMaiorU) -1)[1] ])
        polxU <- c(dif, rep(histSimT$breaks[vfMaiorU], each = 2), difAbs, difAbs)
        polygon(x = polxU, y = polyU, col = rgb(1,0,0, 0.7), border = NA)
        polyL <- c(0, histSimT$counts[rep(which(vfMenorL), each = 2)], 0, 0)
        polxL <- c(rep(histSimT$breaks[vfMenorL], each = 2),-1* difAbs, -1* difAbs, histSimT$breaks[1])
        polygon(x = polxL, y = polyL, col = rgb(1,0,0, 0.7), border = NA)
    }
}    
invisible(list(difSimulado = res, tipo = teste, n_simula = nsim,  p_valor = pvalor))
}
## ## testes
## macho <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112)
## femea <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111)
## x11()
## simulaT(macho, femea, teste = "difere", anima = TRUE)
## set.seed(3)
## m1 = 7.5; m2 = 6; sd1 = 2.2; sd2 = 1.9; ns = 20
## dados1 <- rnorm(ns, m1, sd1)
## dados2 <- rnorm(ns, m2, sd2)
## simDif <- simulaT(dados1, dados2, teste = "difere", anima = TRUE)
## simDif <- simulaT(dados1, dados2, teste = "maior", anima = TRUE)