########################################################################################## ## Modelo Gaussiana para medidas repetidas ############################################### ########################################################################################## ## Carregando pacotes adicionais require(bbmle) require(mvtnorm) ## Função para simular do modelo marginal induzido pelo condicional simula.marginal <- function(beta0, sigma, tau, n.ua, n.rep){ Z <- rep(1, n.rep) Sigma <- tau*(Z%*%t(Z)) + sigma*diag(length(Z)) set.seed(123) y.sim = t(rmvnorm(n.ua, mean=rep(beta0,n.rep), sigma=Sigma)) ID <- rep(1:n.ua, each=n.rep) dados <- data.frame("y"=c(y.sim), ID = ID) return(dados) } ## Simulando um conjunto de dados set.seed(123) dados = simula.marginal(beta0= 5, sigma = 1, tau= 1, n.ua = 10, n.rep = 15) ## Função de verossimilhança verossimilhanca <- function(beta0,sigma, tau, dados){ dados.id <- split(dados,dados$ID) M <- length(dados.id) n.rep <- dim(dados.id[[1]])[1] Z <- rep(1,n.rep) Sigma <- tau*(Z%*%t(Z)) + sigma*diag(n.rep) ll <-c() for(i in 1:M){ ll[i] <- dmvnorm(dados.id[[i]]$y, mean=rep(beta0,n.rep), sigma=Sigma, log=TRUE) } return(-sum(ll)) } ## Maximizando via algoritmo quase-Newton fit <- mle2(verossimilhanca, start=list(beta0 =5, sigma = 1, tau = 1), method="L-BFGS-B", lower=list(beta0 = -Inf, sigma = 1e-5, tau = 1e-5), upper=list(beta0=Inf, tau = Inf, sigma = Inf), data=list(dados=dados)) summary(fit) confint(fit,method="quad") ## Perfil de verossimilhança perfil <- profile(fit) plot(perfil) confint(perfil) ########################################################################################## ## Ajustando via lme4 #################################################################### ########################################################################################## require(lme4) fit2 <- lmer(y ~ (1 |ID), REML=FALSE, data=dados) summary(fit2) ## Via REML fit3 <- lmer(y ~ (1 |ID), REML=TRUE, data=dados) summary(fit3)