[R-br] obter media dos coeficientes de modelos

[FHRB Toledo]

Maurício,

Então, se bem entendi o seu problema acho que aí vaí uma solução:

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


2013/4/10 Maurício Lordêlo <mslordelo em gmail.com>

> Caros,
> No CMR abaixo crio uma matriz de delineamento que estah associada ao nro
> de tratamentos, nro de repeticoes de cada tratamento e o tempo que cada
> resposta (neste caso binaria) foi observada.
> Suponha que o vetor y seja gerado levando em consideracao que no tempo=1,
> a probabilidade eh 0.78; nos tempos>1, esta probabilidade eh 0.35.
>
> 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
>
> y = matrix(data = NA, nrow = n, ncol = 1, byrow = FALSE,dimnames = NULL)
> data1=data.frame(id,time,trat,y)
> for(i in 1:n)
> {
> y[i]=ifelse(data1[i,2]<2,rbinom(1,1,.78),rbinom(1,1,.35))
> }
> data2=data.frame(id,time,trat,y)
> data2
>
> Neste caso tenho formado meu conjunto de dados (data2) já com os valores
> de y gerados.
> Agora vem a minha dificuldade. Preciso rodar um mesmo modelo repetidas
> vezes, como por exemplo:
>
> modelo1=geeglm(y~trat, id=id, data=data2, family=binomial)
>
> e depois obter uma media para cada um dos coeficientes (com seus
> erros-padrao) obtidos.
> Ou seja, terei outros conjuntos de dados gerados (com o mesmo nro de
> tratamentos, repeticoes e tempos)
> sendo que a unica diferenca estarah nos valores de y. Preciso obter a
> média dos coeficientes de cada um
> dos modelos rodados a partir destes diferentes y´s.
> Agradeco quem puder auxiliar.
> Mauricio
>
>
> _______________________________________________
> R-br mailing list
> R-br em listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
> código mínimo reproduzível.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20130411/d8ecd410/attachment.html>