####################################
### Cálculo do poder da anova.MC ###
####################################
# ARGUMENTOS
# N - número de amostras (testes anova.MC a serem feitos)
# n - número de permutações a serem feitos pelo anova.MC
# p - valor crítico de alfa
# dados - data.frame com 4 colunas e b linhas (b = número de blocos)
anova.power=function(dados,N,n,pc=.05){
############################################################################################################################################################
anova.MC=function(dados,n){
aleat=function(dados){ # função que faz as permutações
mc=dados # crio o objeto que receberá as permutações, idêntico aos dados originais
mc$resp=rep(NA,length(dados$resp)) # retiro os valores originais
for(j in 1:length(unique(dados$bloco))){ # número de blocos
mc$resp[mc$bloco==j]=sample(dados$resp[dados$bloco==j]) # permutação dos valores encontrados
}
return(summary(aov(resp~A*B+Error(bloco/(A*B)),mc))[5][[1]][[1]][[4]][1:3]) # extração e cálculo de F
}
result=data.frame(A=rep(NA,n),B=rep(NA,n),A.B=rep(NA,n)) # crio o objeto que receberá os valores de F
for(i in 1:n){ # número de permutações
result[i,]=aleat(dados) # gera n valores de F
}
real=summary(aov(resp~A*B+Error(bloco/(A*B)),data=dados))[5][[1]][[1]][[4]][1:3] # extração e cálculo de F dos dados reais
p=data.frame(fatores="p",
A=(length(result[result[,1]>=real[1],1])+1)/(n+1), # probabilidade de se encontrar o efeito do fator A ao acaso
B=(length(result[result[,2]>=real[2],2])+1)/(n+1), # o mesmo para o fator B
A.B=(length(result[result[,3]>=real[3],3])+1)/(n+1) # e idem para a interação A:B
)
return(p)
}
############################################################################################################################################################
data.gen=function(nb,int,varb,vare,efA,efB,efAB){
d=data.frame(
bloco=rep(1:nb,each=4), ######### crio o fator BLOCO
A=rep(0:1,each=2), ####### idem para o fator A
B=0:1, ##### idem para o fator B
resp=NA ### crio coluna vazia para a variável resposta
)
efb=rnorm(nb,0,varb) # crio o efeito randômico BLOCO com média zero; a variância diz quão diferentes serão os blocos
for(i in 1:nb){ # ciclo número de blocos
for(j in 0:1){ ######### ciclo fator A
for(k in 0:1){ ####### ciclo fator B
d$resp[d$bloco==i&d$A==j&d$B==k]=round(( #
int+ # crio cada valor resposta a partir de uma média
efA*j+ #################### adiciono o efeito do fator A
efB*k+ ############################## idem para o fator B
efAB*j*k+ ##########################idem para a interação A:B
efb[i]+ ########################## idem para o fator BLOCO
rnorm(1,0,vare) ##### adiciono o erro associado a cada observação
),2) }}} ############################# arredonda os dados
#
return(d)} # retorna a tabela de dados gerados
############################################################################################################################################################
medias=aggregate(dados$resp,list(dados$A,dados$B),mean) # calculo as médias por tratamento
desvios=aggregate(dados$resp,list(dados$A,dados$B),sd) # idem para os desvios
medias.bloco=aggregate(dados$resp,list(dados$bloco),mean) # calculo as médias por bloco
nb=length(unique(dados$bloco)) # conto o número de blocos do desenho
int=medias[1,3] # intercepto ou média dos controles
varb=sd(medias.bloco[,2]) # desvios entre os blocos
efA=medias[2,3]-medias[1,3] # calculo o efeito do tratamento para o fator A
efB=medias[3,3]-medias[1,3] # idem para o fator B
efAB=medias[4,3]-medias[2,3]-medias[3,3]+medias[1,3] # idem para a interação entre A e B
efbloco=numeric()
for(i in 1:nb){
efbloco[i]=medias.bloco[i,2]-mean(dados$resp) # calculo o efeito do bloco
}
residuos=numeric()
for(i in 1:length(dados$resp)){
residuos[i]=dados[i,4]-int-efA*dados[i,2]-efB*dados[i,3]-efAB*dados[i,2]*dados[i,3]-efbloco[dados[i,1]] # cálculo para encontrar os resíduos
}
vare=sd(residuos) # desvio padrão dos resíduos
d=list()
for(i in 1:N){
d[[i]]=data.gen(nb,int,varb,vare,efA,efB,efAB)
}
result=data.frame(pA=rep(NA,N),pB=rep(NA,N),pAB=rep(NA,N))
for(i in 1:N){
result[i,]=anova.MC(d[[i]],n)[2:4]
}
power=data.frame(row.names="poder(%)",
A=length(result[result$pA<=pc,1])*100/length(result$pA),
B=length(result[result$pB<=pc,2])*100/length(result$pB),
AB=length(result[result$pAB<=pc,3])*100/length(result$pAB)
)
return(power)
}
# testando...
anova.power(dados,2,999)