#Codigo da funcao "edm" com intervalos de confianca e graficos

edm <- function (x, y, data=data, IC=FALSE, estima=FALSE, N=10)
{
  #chama os vetores contendo as variaveis de interesse "x" e "y" e produz um "data frame".
  edm.call <- Call <- match.call()
  x <- as.character(edm.call[[2L]])
  y <- as.character(edm.call[[3L]])
  dados <- data.frame(data[x], data[y])
  #numero de reamostragens em caso de estimativa de erros (default=10).
  N <- N
  #testa se a variavel independente pertence a classe fator, e faz a transformacao
  #em caso negativo
  if (!is.factor(dados[, 1]))
  {
    dados[, 1] <- as.factor(dados[, 1])
  }
  #estima o numero de amostras e de replicas.
  n <- length(unique(dados[, 1]))
  k <- length(dados[, 1])/n
  #realiza uma anova com as variaveis de interesse, das quais sao retiradas os valores
  #de desvios quadrados medios, utilizados no calculo de erro e repetibilidade.
  funcao <- formula(dados[, 2] ~ dados[, 1])
  varis <- anova(aov(funcao, data = dados))
  num.df <- n-1
  denom.df <- n*(k-1)
  MSa <- varis [1, 3]
  Sw <- MSw <- varis [2, 3]
  Sa <- (MSa - MSw)/k
  #calcula o valor de repetibilidade.
  r <- Sa/(Sw + Sa)
  #calcula o valor de EM.
  EM <- Sw/(Sw + Sa)
  #argumento que chama a estimativa de multiplas repetibilidades.
  if (estima == TRUE)
  {
    result3 <- data.frame(rep(2:k, each=N), rep(NA, ((k-1)*N)), rep(NA, ((k-1)*N)))
    #produz N estimativas de repetibilidade para subamostras das 2:k replicas das
    #medicoes.
    for(l in 2:k){
      result2 <- data.frame(rep(NA, N), rep(NA, N))
      #refaz a estimativa de repetibilidade N vezes.
      for(j in 1:N){
        result <- data.frame(rep(1:n, each=l) , rep(NA, n*l))               
        #loop que calcula a repetibilidade
        for(i in 1:n){
          result[, 2][which(result[, 1]==i)] = sample((dados[which(dados[, 1]==i), ][, 2])
          , size=l)
          result[, 1] <- as.factor(result[, 1])
        }
        funcao <- anova(aov(result[, 2] ~ result[, 1], data=dados))
        n <- length(unique(result[, 1]))
        rep <- length(result[, 1])/n
        MSa <- funcao [1, 3]
        Sw <- MSw <- funcao [2, 3]
        Sa <- (MSa - MSw)/k
        r <- Sa/(Sw + Sa)
        result2[, 1][j] = r
        result2[, 2][j] = l
      }  
      result3[, 2][which(result3[, 1]==l)] = result2[, 1][1:N]
      result3[, 3][which(result3[, 1]==l)] = result2[, 2][1:N]
    }
    #abre uma janela e plota grafico de caixa com as estimativas de repetibilidade.
    x11()
    boxplot(result3[, 2] ~ result3[, 1])
    mtext("Repetibilidade", side=2, cex=1.5, line=2.5, las=0)
    mtext("Número de réplicas", side=1, cex=1.5, line=2.5)
  }
  #argumento que ativa o calculo de intervalos de confianca
  if (IC == TRUE) 
  {
    #calculo dos intervalos de confianca segundo Bonett(2002).
    F.stat <- varis [1, 4]
    low.F <- qf(0.05/2, num.df, denom.df, lower.tail = FALSE)
    up.F <- qf(0.05/2, denom.df, num.df, lower.tail = FALSE)
    FL <- F.stat/low.F
    FU <- F.stat * up.F
    ICLow <- (FL - 1)/(FL + k - 1)
    ICUp <- (FU - 1)/(FU + k - 1)
    ICdif <- ICUp-ICLow
    #retorna a lista com valores calculados pela funcao
    return(list(edm = EM, repetibilidade = r, ICsup = ICUp, ICinf = ICLow, ICdif = ICdif, 
    individuos = n, replicas = k))
  }
  #retorna resultados simplificados, sem calculo de ICs.
  else
    return(list(edm = EM, repetibilidade = r, individuos = n, replicas = k))
}