########################################################################################## ### Modelos Mistos ####################################################################### ########################################################################################## ## Modelo com apenas um intercepto aleatório Gaussiano simula.modelo <- function(n.individuo, n.rep, beta0, sigma, tau){ b <- rnorm(n.individuo, 0, sd=sigma) mu <- beta0 + b y <- matrix(NA, ncol=n.individuo, nrow=n.rep) for(i in 1:n.individuo){ y[,i] <- rnorm(n.rep, mu, sd=tau) } ID <- rep(1:n.individuo,each=n.rep) saida <- data.frame(y = c(y), ID = ID) return(saida) } dados = simula.modelo(n.individuo = 10, n.rep = 10, beta0= 5, sigma = 1, tau=0.5) ## mostrando os dados boxplot(y ~ as.factor(ID),data=dados) abline(h = 5) ## desvio da média tem média 0 mean(tapply(dados$y, dados$ID, mean) - mean(dados$y)) ## Verossimilhança de um individuo Z <- rep(1, n.rep) sigma*diag(n.rep) ########################################################################################## ## Detalhado ############################################################################# ########################################################################################## ## Simulando o b b <- rnorm(2, 0,sd=tau) Z <- rep(1,n.rep) ## Simulando o y y1 <- rnorm(5,beta0 + Z*b[1], sd= sigma) y2 <- rnorm(5,beta0 + Z*b[2], sd= sigma) ## Matriz de covariancia de um individuo Sigma <- tau*(Z%*%t(Z)) + sigma*diag(5) ## Quem é a função de verossimilhança verossimilhanca <- function(beta0,sigma, tau, dados){ dados.id <- split(dados,dados$ID) M <- length(dados.id) Z <- rep(1,M) Sigma <- tau*(Z%*%t(Z)) + sigma*diag(M) ll <-c() for(i in 1:M){ ll[i] <- dmvnorm(dados.id[[i]]$y, mean=rep(beta0,M), sigma=Sigma, log=TRUE) } return(-sum(ll)) } ########################################################################################## fit <- mle2(verossimilhanca, start=list(beta0 =5, sigma = 1, tau = 0.5), 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) ##########################################################################################