n.trat <- 3 # número de tratamentos
n.rep <- 2 # número de repetições
n.time <- 4 # número de tempos
n.ind <- n.trat * n.rep # número de indivíduos
n <- n.trat * n.rep * n.time # número total de registros
id <- rep(1:n.ind, each = n.time) # identificação do indivíduo
time <- rep(1:n.time, times = n.ind) # identificação do tempo
trat <- rep(factor(LETTERS[1:n.trat]), each = n.rep*n.time) # tratamentos
dados <- data.frame(id, time, trat) # estrutura do delineamento
y <- function(x) ifelse(dados$time < 2, rbinom(1, 1, .78), rbinom(1, 1, .35)) ## funcao que simula um ensaio
nsmc <- 1000 # numero de simulacoes de monte carlo
varios <- replicate(nsmc, cbind(dados, y = y(x)), simplify = FALSE) # repete o ensaio nsmc vezes
modelo <- function( plan ) glm(y ~ trat, data = plan, family = binomial)$coef # funcao que faz o ajuste
ajustes <- do.call(rbind, lapply(varios, modelo)) # varios ajustes
medias <- apply(ajustes, 2, mean) # médias dos coeficientes
Primeiro criei uma função que atribui os resultados ao vetor y, segundo suas premissas. O replicate() repete esse procedimento por nsmc vezes e armazena cada uma das simulações em uma posição da lista.
Depois disso é só usar a sua função de ajuste (usei a glm normal) em uma função e aplica com o lapply(). Fiz um artifício de já organizar os coeficientes do modelo numa matriz, uma linha para cada simulação e uma coluna por parâmetro.
Espero ter ajudado.
abraço,
FH