questao2.R
viccj
2025-09-10
###################################
## Q2 — Modelos Mistos ##
## Autora: Victória Souza ##
#################################### Pacotes: modelos mistos (lme/lmer) e gráficos.
library (lme4)## Carregando pacotes exigidos: Matrixlibrary (nlme)##
## Anexando pacote: 'nlme'## O seguinte objeto é mascarado por 'package:lme4':
##
## lmListlibrary (ggplot2)
# Ler a base de dados (delimitado por tab)
dados <- read.delim("produtividade.txt")
# Inspeção rápida do cabeçalho
head (dados) #olhar o cabeçalho da base de dados## producao bloco irrigacao densidade fertilizante
## 1 90 A controle baixo N
## 2 95 A controle baixo P
## 3 107 A controle baixo NP
## 4 92 A controle medio N
## 5 89 A controle medio P
## 6 92 A controle medio NP# Visualização geral
dados ## producao bloco irrigacao densidade fertilizante
## 1 90 A controle baixo N
## 2 95 A controle baixo P
## 3 107 A controle baixo NP
## 4 92 A controle medio N
## 5 89 A controle medio P
## 6 92 A controle medio NP
## 7 81 A controle alto N
## 8 92 A controle alto P
## 9 93 A controle alto NP
## 10 80 A irrigado baixo N
## 11 87 A irrigado baixo P
## 12 100 A irrigado baixo NP
## 13 121 A irrigado medio N
## 14 110 A irrigado medio P
## 15 119 A irrigado medio NP
## 16 78 A irrigado alto N
## 17 98 A irrigado alto P
## 18 122 A irrigado alto NP
## 19 83 B controle baixo N
## 20 80 B controle baixo P
## 21 95 B controle baixo NP
## 22 98 B controle medio N
## 23 98 B controle medio P
## 24 106 B controle medio NP
## 25 74 B controle alto N
## 26 81 B controle alto P
## 27 74 B controle alto NP
## 28 102 B irrigado baixo N
## 29 109 B irrigado baixo P
## 30 105 B irrigado baixo NP
## 31 99 B irrigado medio N
## 32 94 B irrigado medio P
## 33 123 B irrigado medio NP
## 34 136 B irrigado alto N
## 35 133 B irrigado alto P
## 36 132 B irrigado alto NP
## 37 85 C controle baixo N
## 38 88 C controle baixo P
## 39 88 C controle baixo NP
## 40 112 C controle medio N
## 41 104 C controle medio P
## 42 91 C controle medio NP
## 43 82 C controle alto N
## 44 78 C controle alto P
## 45 94 C controle alto NP
## 46 60 C irrigado baixo N
## 47 104 C irrigado baixo P
## 48 114 C irrigado baixo NP
## 49 90 C irrigado medio N
## 50 118 C irrigado medio P
## 51 113 C irrigado medio NP
## 52 119 C irrigado alto N
## 53 122 C irrigado alto P
## 54 136 C irrigado alto NP
## 55 86 D controle baixo N
## 56 78 D controle baixo P
## 57 89 D controle baixo NP
## 58 79 D controle medio N
## 59 86 D controle medio P
## 60 87 D controle medio NP
## 61 85 D controle alto N
## 62 89 D controle alto P
## 63 83 D controle alto NP
## 64 73 D irrigado baixo N
## 65 114 D irrigado baixo P
## 66 114 D irrigado baixo NP
## 67 109 D irrigado medio N
## 68 131 D irrigado medio P
## 69 126 D irrigado medio NP
## 70 116 D irrigado alto N
## 71 136 D irrigado alto P
## 72 133 D irrigado alto NP# Tornar colunas acessíveis por nome, sem a necessidade de especificar fatores (dados$)
attach (dados)
#------------------------------------------------
# SELEÇÃO DE FIXOS COM ML (comparação de modelos)
#------------------------------------------------
# Modelo Nulo (apenas intercepto nos fixos)
# Ajustando o modelo nulo, mantendo apenas produção como variável resposta e nenhuma variável explicativa.
# Aleátorios aninhados (~1|bloco/irrigacao/densidade) = (1|bloco)+(1|bloco:irrigacao)+(1|bloco:irrigacao:densidade)
# method = "ML" -> Maximum Likelihood, permite que a comparação entre modelos que diferem nos fatores fixos seja feita usando ANOVA e AIC/BIC.
m.nulo = lme (fixed = producao ~ 1,
random = ~ 1|bloco/irrigacao/densidade,
method = "ML") # retira o efeito de blocos (e níveis aninhados) ao tratá-los como ALEATÓRIOS.
# Modelo 1: fatores isolados + todas as interações entre irrigação, densidade e fertilizante.
m1 = lme (fixed = producao ~ irrigacao*densidade*fertilizante,
random = ~ 1|bloco/irrigacao/densidade,
method = "ML") #ajustando o modelo com produção em função de irrigação, densidade e fertilidade, com interação entre eles (*)
# LRT (via anova), para comparar os modelos (m.nulo vs m1) por verossimilhança e investigar se os fatores fixos explicam algo além do modelo nulo
anova (m.nulo, m1) ## Model df AIC BIC logLik Test L.Ratio p-value
## m.nulo 1 5 605.4183 616.8016 -297.7091
## m1 2 22 573.5108 623.5974 -264.7554 1 vs 2 65.9075 <.0001Interpretação: LRT χ²=65,91 (gl=17), p<0,0001 (significativo): o modelo com fixos (m1) explica a produção melhor que o nulo. Prosseguir com m1.
# ANOVA por termo do Modelo 1, para ver a contribuição de cada termo e simplificar o modelo, se preciso.
anova (m1) ## numDF denDF F-value p-value
## (Intercept) 1 36 2674.6636 <.0001
## irrigacao 1 3 30.9211 0.0115
## densidade 2 12 3.7842 0.0532
## fertilizante 2 36 11.4493 0.0001
## irrigacao:densidade 2 12 5.9119 0.0163
## irrigacao:fertilizante 2 36 5.5204 0.0081
## densidade:fertilizante 4 36 0.8826 0.4841
## irrigacao:densidade:fertilizante 4 36 0.6795 0.6107Interpretação: A interação tripla irrigação×densidade×fertilizante não é significativa (p=0,61): manter apenas interações duplas.
# Modelo 2: mantém principais + todas as duplas
m2 = lme (fixed = producao ~ (irrigacao+densidade+fertilizante)^2,
random = ~ 1|bloco/irrigacao/densidade,
method = "ML") #^2 mantém apenas até as interações duplas
# LRT m1 vs m2
anova (m1, m2) ## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 22 573.5108 623.5974 -264.7554
## m2 2 18 569.0046 609.9845 -266.5023 1 vs 2 3.493788 0.4788Interpretação: m2 (só duplas) é tão bom quanto m1 (p=0,48; não significativo) e mais simples: optar pelo m2.
# ANOVA (m2): identificar duplas não significativas
anova (m2) ## numDF denDF F-value p-value
## (Intercept) 1 40 2872.7400 <.0001
## irrigacao 1 3 33.2110 0.0104
## densidade 2 12 4.0645 0.0449
## fertilizante 2 40 11.4341 0.0001
## irrigacao:densidade 2 12 6.3499 0.0132
## irrigacao:fertilizante 2 40 5.5131 0.0077
## densidade:fertilizante 4 40 0.8815 0.4837Interpretação: A interação dupla densidade×fertilizante não é significativa (p=0,48): remover esse termo (gerar m3).
# Modelo 3: (sem densidade:fertilizante): principais + (irrigação×densidade) + (irrigação×fertilizante).
m3 = lme (fixed = producao ~ irrigacao+densidade+fertilizante+
irrigacao:densidade+irrigacao:fertilizante,
random = ~ 1|bloco/irrigacao/densidade,
method = "ML")
# LRT m2 vs m3
anova (m2, m3) ## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 18 569.0046 609.9845 -266.5023
## m3 2 14 565.1933 597.0667 -268.5967 1 vs 2 4.188774 0.3811Interpretação: m3 é tão bom quanto m2 (p=0,38; não significativo), é mais simples e tem AIC menor: optar pelo m3.
# ANOVA (m3)
anova (m3) ## numDF denDF F-value p-value
## (Intercept) 1 44 3070.8771 <.0001
## irrigacao 1 3 35.5016 0.0095
## densidade 2 12 4.3448 0.0381
## fertilizante 2 44 11.2013 0.0001
## irrigacao:densidade 2 12 6.7878 0.0107
## irrigacao:fertilizante 2 44 5.4008 0.0080Interpretação: todas as variaveis são significativas agora. A produção depende dos três fatores e os efeitos de densidade e de fertilizante mudam com a irrigação.
# Sumário do m3: componentes de variância (aleatórios) e coeficientes (fixos, relativos aos níveis de referência).
summary (m3) ## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 565.1933 597.0667 -268.5967
##
## Random effects:
## Formula: ~1 | bloco
## (Intercept)
## StdDev: 0.0005261003
##
## Formula: ~1 | irrigacao %in% bloco
## (Intercept)
## StdDev: 1.716888
##
## Formula: ~1 | densidade %in% irrigacao %in% bloco
## (Intercept) Residual
## StdDev: 5.722413 8.718327
##
## Fixed effects: producao ~ irrigacao + densidade + fertilizante + irrigacao:densidade + irrigacao:fertilizante
## Value Std.Error DF t-value p-value
## (Intercept) 82.08333 4.756285 44 17.257867 0.0000
## irrigacaoirrigado 27.80556 6.726403 3 4.133793 0.0257
## densidadebaixo 4.83333 5.807347 12 0.832279 0.4215
## densidademedio 10.66667 5.807347 12 1.836754 0.0911
## fertilizanteNP 4.33333 3.835552 44 1.129781 0.2647
## fertilizanteP 0.91667 3.835552 44 0.238992 0.8122
## irrigacaoirrigado:densidadebaixo -29.75000 8.212829 12 -3.622382 0.0035
## irrigacaoirrigado:densidademedio -19.66667 8.212829 12 -2.394628 0.0338
## irrigacaoirrigado:fertilizanteNP 16.83333 5.424290 44 3.103325 0.0033
## irrigacaoirrigado:fertilizanteP 13.50000 5.424290 44 2.488805 0.0167
## Correlation:
## (Intr) irrgcr dnsddb dnsddm frtlNP frtlzP
## irrigacaoirrigado -0.707
## densidadebaixo -0.610 0.432
## densidademedio -0.610 0.432 0.500
## fertilizanteNP -0.403 0.285 0.000 0.000
## fertilizanteP -0.403 0.285 0.000 0.000 0.500
## irrigacaoirrigado:densidadebaixo 0.432 -0.610 -0.707 -0.354 0.000 0.000
## irrigacaoirrigado:densidademedio 0.432 -0.610 -0.354 -0.707 0.000 0.000
## irrigacaoirrigado:fertilizanteNP 0.285 -0.403 0.000 0.000 -0.707 -0.354
## irrigacaoirrigado:fertilizanteP 0.285 -0.403 0.000 0.000 -0.354 -0.707
## irrgcrrgd:dnsddb irrgcrrgd:dnsddm irr:NP
## irrigacaoirrigado
## densidadebaixo
## densidademedio
## fertilizanteNP
## fertilizanteP
## irrigacaoirrigado:densidadebaixo
## irrigacaoirrigado:densidademedio 0.500
## irrigacaoirrigado:fertilizanteNP 0.000 0.000
## irrigacaoirrigado:fertilizanteP 0.000 0.000 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.58166960 -0.51480884 0.07893407 0.60157076 2.19570825
##
## Number of Observations: 72
## Number of Groups:
## bloco irrigacao %in% bloco
## 4 8
## densidade %in% irrigacao %in% bloco
## 24Interpretação: Ser irrigado é o que mais aumenta a produção (coef. de irrigaçãoirrigado = +27,8)
# Modelo 4, apenas para confirmar se esse é o MMA (modelo mínimo adequado): retirando a interação (irrigação:fertilizante), o modelo piora?
m4 = lme (fixed = producao ~ irrigacao+densidade+fertilizante+irrigacao:densidade,
random = ~ 1|bloco/irrigacao/densidade,
method = "ML")
# ANOVA(m4) apenas para referência dos termos em m4
anova (m4)## numDF denDF F-value p-value
## (Intercept) 1 46 3169.893 <.0001
## irrigacao 1 3 36.646 0.0090
## densidade 2 12 4.485 0.0351
## fertilizante 2 46 9.167 0.0004
## irrigacao:densidade 2 12 7.007 0.0096# LRT m3 vs m4
anova (m3, m4) ## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 14 565.1933 597.0667 -268.5967
## m4 2 12 572.3373 599.6573 -274.1687 1 vs 2 11.14397 0.0038Interpretação: Remover irrigação×fertilizante piora o ajuste (LRT p=0,0038): manter a interação; m3 segue como modelo final.
#---------------------------------------------
# O MODELO MISTO É MESMO NECESSÁRIO?
#---------------------------------------------
# Modelo linear simples (sem aleatório de bloco) com os principais e uma das duplas.
ml = lm (producao ~ irrigacao+densidade+fertilizante+irrigacao:densidade) #Modelo sem fator aleatório
# LRT misto (m3) vs (ml)
anova.lme (m3, ml) ## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 14 565.1933 597.0667 -268.5967
## ml 2 9 569.9342 590.4242 -275.9671 1 vs 2 14.74085 0.0115Interpretação: Significativo (LRT p=0,0115): o fator aleatório (bloco) melhora o ajuste, ou seja, o modelo misto é necessário. É coerente retirar o efeito dos blocos.
#---------------------------------------------
# DIAGNÓSTICOS DO MODELO (PRESSUPOSTOS)
#---------------------------------------------
# Resíduos padronizados vs valores ajustados: avaliar homocedasticidade e tendência.
plot (m3) 
Interpretação: Dispersão aleatória em torno de zero, sem “funil” visível: homocedasticidade ok.
# Observado vs ajustado: verificar alinhamento (qualidade de ajuste).
plot (m3, producao ~ fitted(.)) 
Interpretação: a variavel resposta é razoavelmente uma função linear dos valores ajustados: ajuste adequado.
# QQ-plot dos resíduos, por bloco: checar normalidade dos erros em cada bloco.
qqnorm (m3, ~ resid(.)| bloco) 
Interpretação: Em cada bloco, os pontos formam um padrão aproximadamente retilíneo, com pequenas curvaturas apenas nas caudas. Isso é compatível com normalidade aceitável dos resíduos.
# Ajuste final em REML (inferência): reportar F/p dos fixos e variâncias (bloco, residual).
m3_REML <- update(m3, method="REML")
# Tabela final de F/p
anova(m3_REML) ## numDF denDF F-value p-value
## (Intercept) 1 44 2674.6034 <.0001
## irrigacao 1 3 30.9204 0.0115
## densidade 2 12 3.7842 0.0532
## fertilizante 2 44 11.9239 0.0001
## irrigacao:densidade 2 12 5.9119 0.0163
## irrigacao:fertilizante 2 44 5.7492 0.0061Interpretação: Em REML, irrigação (p=0,0115) e fertilizante (p=0,0001) são significativos; densidade é marginal (p=0,053). As interações duplas com irrigação permanecem significativas. Ou seja, a produção depende de irrigação e fertilizante, e os efeitos de densidade e de fertilizante mudam com a irrigação
# Componentes de variância
VarCorr(m3_REML) ## Variance StdDev
## bloco = pdLogChol(1)
## (Intercept) 3.097876e-07 0.0005565857
## irrigacao = pdLogChol(1)
## (Intercept) 3.930759e+00 1.9826142192
## densidade = pdLogChol(1)
## (Intercept) 4.980311e+01 7.0571318966
## Residual 8.291915e+01 9.1059950596Interpretação: variância de bloco praticamente nula, a variabilidade concentra-se no nível residual (≈61% (82,91/136,65)). Confirma que o efeito de blocos foi adequadamente controlado.Obs: percentuais: 3.1e-7, 3.93, 49.80, 82.92 ⇒ total ≈136.65.
# VISUALIZAÇÃO DAS INTERAÇÕES - Figuras Finais# 1) Interação Fertilizante × Irrigação
if (interactive() && !isTRUE(getOption("knitr.in.progress"))) windows(10, 8)
interaction.plot(fertilizante, irrigacao, producao)
Interpretação: Irrigação eleva a produção; sob irrigação, NP tende a superar N e P; as linhas não paralelas indicam interação entre fertilizante e irrigação.
# 2) Interação Densidade × Irrigação
if (interactive() && !isTRUE(getOption("knitr.in.progress"))) windows(10, 8)
interaction.plot(densidade, irrigacao, producao)
Interpretação: Sob irrigação, densidades maiores tendem a apresentar médias mais altas; as linhas não paralelas indicam interação entre densidade e irrigação.
Conclusão final: A produção depende de irrigação (F(1,3)=30,92; p=0,0115) e de fertilizante (F(2,44)=11,92; p=0,0001). O efeito de densidade é marginal no geral (F(2,12)=3,78; p≈0,053), mas há interações significativas com irrigação: irrigação×densidade (p=0,016) e irrigação×fertilizante (p=0,006). Isso indica que os efeitos de densidade e de fertilizante dependem do regime de irrigação. A variância de bloco é praticamente nula, confirmando que seu efeito foi adequadamente controlado.