## Carregando pacote para usar os dados library(faraway) data(pulp) ## Estimando o modelo de efeitos fixos lmod <- aov(bright ~ operator, pulp) summary(lmod) ## os coeficientes são interpretados de acordo com a opção de contrastes: coef(lmod) ## Mudando a interpretação dos parâmetros de média op <- options(contrasts=c("contr.sum", "contr.poly")) options(op) lmod <- aov(bright ~ operator, pulp) ## No ajuste de modelo considerano efeitos aleatórios... ## a estimativa da variãncia dos tratamentos pelo método da anova é: (0.447-0.106)/5 ## Ajustando novamente o modelo de efeitos aleatórios com outro **metodo de estimação** library(lme4) ## Método de estimação: REML mmod <- lmer(bright ~ 1+(1|operator), pulp) summary(mmod) ## Método de estimação: ML smod <- lmer(bright ~ 1+(1|operator), pulp, REML=FALSE) summary(smod) ## Teste do efeito dos "operadores" (tratamento) pelo TRV ## Para o TRV precisamos ajustar o modelo nulo nullmod <- lm(bright ~ 1, pulp) ## TRV supondo problema "bem comportado" (TRVest <- as.numeric(2*(logLik(smod)-logLik(nullmod)))) #pchisq(2.5684,1,lower=FALSE) pchisq(TRVest,1,lower=FALSE) ## ... porém o teste que queremos de $\sigma^2_T = 0" não é "bem comportado" ## bootrap (paramétrico) ##y <- simulate(nullmod) ## objeto para guardar e estística de teste em cada simulação lrstat <- numeric(1000) ## construção de dist. amostral da estatística do TRV for(i in 1:1000){ y <- unlist(simulate(nullmod)) bnull <- lm(y ~ 1) balt <- lmer(y ~ 1 + (1|operator),pulp,method="ML") lrstat[i] <- as.numeric(2*(logLik(balt)-logLik(bnull))) } ## distribuição amostral hist(lrstat, prob=T); lines(density(lrstat)) curve(dchisq(x, df=1), add=T) ## falta de ajuste da distribuição chi-quadrado para o TRVneste caso prop0 <- mean(lrstat < 0.00001) ## cálculo do p-valor "empírico" (bootstrap) mean(lrstat > 2.5684) mean(lrstat > TRVest) ## "margem de erro" parea o meu p-valor sqrt(prop0*(1-prop0)/1000) ## mesmo considerando a margem de erro para o p0-valor, rejeitaríamos H0 ## estimativas dos efeitos pelo modelo aleatório (EfAl <- ranef(mmod)$operator) ## estimativas dos efeitos pelo modelo fixo (cc <- model.tables(lmod)) ## notamos que as do efeito aleatório sao menos variáveis ## (shrinkage) cc[[1]]$operator/ranef(mmod)$operator ## Médias estimadas (pradições) fixef(mmod) fixef(mmod)+EfAl ## Diagnśticos ## notar que há "diferentes" possíveis resíduos qqnorm(resid(mmod),main="") plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals") abline(0,0) ## data(penicillin) summary(penicillin) op <- options(contrasts=c("contr.sum", "contr.poly")) lmod <- aov(yield ~ blend + treat, penicillin) summary(lmod) coef(lmod) mmod <- lmer(yield ~ treat + (1|blend), penicillin) summary(mmod) options(op) ranef(mmod)$blend anova(mmod) amod <- lmer(yield ~ treat + (1|blend), penicillin, method="ML") nmod <- lmer(yield ~ 1 + (1|blend), penicillin, method="ML") anova(amod,nmod) lrstat <- numeric(1000) for(i in 1:1000){ ryield <- unlist(simulate(nmod)) nmodr <- lmer(ryield ~ 1 + (1|blend), penicillin, method="ML") amodr <- lmer(ryield ~ treat + (1|blend), penicillin, method="ML") lrstat[i] <- 2*(logLik(amodr)-logLik(nmodr)) } plot(qchisq((1:1000)/1001,3),sort(lrstat),xlab=expression(chi[3]^2),ylab="Simulated LRT") abline(0,1) mean(lrstat > 4.05) rmod <- lmer(yield ~ treat + (1|blend), penicillin) nlmod <- lm(yield ~ treat, penicillin) 2*(logLik(rmod)-logLik(nlmod,REML=TRUE)) lrstatf <- numeric(1000) for(i in 1:1000){ ryield <- unlist(simulate(nlmod)) nlmodr <- lm(ryield ~ treat, penicillin) rmodr <- lmer(ryield ~ treat + (1|blend), penicillin) lrstatf[i] <- 2*(logLik(rmodr)-logLik(nlmodr,REML=TRUE)) } mean(lrstatf < 0.00001) cs <- lrstatf[lrstatf > 0.00001] ncs <- length(cs) plot(qchisq((1:ncs)/(ncs+1),1),sort(cs),xlab=expression(chi[1]^2),ylab="Simulated LRT") abline(0,1) mean(lrstatf > 2.7629) data(irrigation) summary(irrigation) lmod <- lmer(yield ~ irrigation * variety + (1|field) +(1|field:variety),data=irrigation) logLik(lmod) lmodr <- lmer(yield ~ irrigation * variety + (1|field),data=irrigation) logLik(lmodr) summary(lmodr) anova(lmodr) plot(fitted(lmodr),resid(lmodr),xlab="Fitted",ylab="Residuals") qqnorm(resid(lmodr),main="") mod <- lm(yield ~ irrigation * variety, data=irrigation) anova(mod) data(eggs) summary(eggs) cmod <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs) summary(cmod) cmodr <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs) anova(cmod,cmodr) VarCorr(cmodr) lrstat <- numeric(1000) for(i in 1:1000){ rFat <- unlist(simulate(cmodr)) nmod <- lmer(rFat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs) amod <- lmer(rFat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs) lrstat[i] <- 2*(logLik(amod)-logLik(nmod)) } mean(lrstat < 0.00001) 2*(logLik(cmod)-logLik(cmodr)) mean(lrstat > 1.6034) data(abrasion) matrix(abrasion$material,4,4) lmod <- aov(wear ~ material + run + position, abrasion) summary(lmod) mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion) anova(mmod) summary(mmod) data(jsp) jspr <- jsp[jsp$year==2,] plot(jitter(math)~jitter(raven),data=jspr,xlab="Raven score",ylab="Math score") boxplot(math ~ social, data=jspr,xlab="Social class",ylab="Math score") glin <- lm(math ~ raven*gender*social,jspr) anova(glin) glin <- lm(math ~ raven*social,jspr) anova(glin) glin <- lm(math ~ raven+social,jspr) summary(glin) table(jspr$school) mmod <- lmer(math ~ raven*social*gender+(1|school)+(1|school:class),data=jspr) anova(mmod) jspr$craven <- jspr$raven-mean(jspr$raven) mmod <- lmer(math ~ craven*social+(1|school)+(1|school:class),jspr) summary(mmod) qqnorm(resid(mmod),main="") plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals") qqnorm(ranef(mmod)$school[,1],main="School effects") qqnorm(ranef(mmod)$"school:class"[,1],main="Class effects") adjscores <- ranef(mmod)$school[,1] rawscores <- coef(lm(math ~ school-1,jspr)) rawscores <- rawscores-mean(rawscores) plot(rawscores,adjscores) sint <- c(9,14,29) text(rawscores[sint],adjscores[sint]+0.2,c("9","15","30")) schraven <- lm(raven ~ school, jspr)$fit mmodc <- lmer(math ~ craven*social+schraven*social+(1|school)+(1|school:class),jspr) anova(mmodc)