#Exercicio 3#
terrenos<-read.table("terrenos.txt", header=T)
head(terrenos)
#ajustamento do modelo com 1 factor de efeitos fixos (variedade) e 1 factor de efeitos aleatórios (terreno), método REML#
terrenoslme1<-lme(rend~variedade, random=~1|terreno, data=terrenos)
summary(terrenoslme1)

#cálculo da estatística do teste aos efeitos fixos#

#matriz de covariancias estimadas dos estimadores dos efeitos fixos#
vcov(terrenoslme1)
c<-vcov(terrenoslme1)
#Construção da matriz L#
l<-matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),4,3)
l
#vector com as estimativas dos efeitos fixos#
b<-c(1.556,-0.0238462, -0.3890769,-0.3778462)
#cálculo do numerador da estatística F do teste aos efeitos fixos#
numerador<-t(b)%*%l%*%solve(t(l)%*%c%*%l)%*%t(l)%*%b
#informação sobre a matriz L, nomeadamente sobre a característica da matriz#
qr(l)
#cálculo da estatística do teste aos efeitos fixos#
F<-numerador/3
F
#teste aos efeitos fixos directamente com o comando anova#
anova(terrenoslme1)

terrenos.lme1<-lme(rend~variedade, random=~1|terreno, data=terrenos)
plot(terrenos.lme1)
qqnorm(terrenos.lme1, ~resid(.))
qqnorm(terrenos.lme1, ~ranef(.))



#Exercicio 4#
#com o lme4, função lmer#
#para um estudo descritivo preliminar#
plot.design(Machines)
attach(Machines)
interaction.plot(Machine, Worker, score)
detach()
#Ajustamento do modelo descrito na alínea a): modelo linear misto, 1 factor de efeitos fixos (máquina, com 3 níveis, 1 factor de efeitos aleatórios (trabalhador, com 6 níveis), com interacção#
machines1r<-lmer(score~Machine+(1|Worker)+(1|Worker:Machine), data=Machines)
#para ver os resultados do ajustamento#
summary(machines1r)
#Para ver as estimativas de máxima verosimlhança#
summary(lmer(score~Machine+(1|Worker)+(1|Worker:Machine), REML=FALSE, data=Machines))
#comandos para efectuar o teste de hipóteses à componente de variância da interacção#
#ajustamento do modelo sem interacção#
machines2r<-lmer(score~Machine+(1|Worker), data=Machines)
summary(machines2r)
#log-verosimilhança restrita do modelo com interacção"
lr1<-logLik(machines1r)
lr1
#log-verosimilhança restrita do modelo sem interacção"
lr0<-logLik(machines2r)
lr0
#cálculo da estatistica do teste de razão de verosimilhanças#
estcal<-2*(lr1-lr0)
estcal
#p-value#
1-pchisq(estcal,1)

#comandos para efectuar o teste de hipóteses à componente de variância associada factor trabalhador#
#ajustamento do modelo sem factor trabalhador#
machines3r<-lmer(score~Machine+(1|Worker:Machine), data=Machines)
#log-verosimilhança restrita do modelo sem factor trabalhador"
lr30<-logLik(machines3r)
lr30
#cálculo da estatistica do teste de razão de verosimilhanças#
estcal2<-2*(lr1-lr30)
estcal2
#p-value#
1-pchisq(estcal2,1)
#p-value corrigido de acordo com a verdadeira distribuição da estatistica do teste#
0.5*(1-pchisq(estcal2,1))

#Para efectuar o teste de hipóteses aos efeitos fixos do modelo, consultar o valor da estatistica do teste calculada#
anova(machines1r)
#valor que numa distribuição F com 2 e 10 graus de liberdade deixa à direita uma região de probabilidade 0.05#
qf(1-0.05,2,10)
#p-value#
1-pf(20.576,2,10)

#gráficos de diagnóstico para validação de alguns dos pressupostos do modelo#
#gráfico dos resíduos vs valores ajustados"
plot(machines1r)
"qqnorm dos residuos#
residuos<-resid(machines1r)
qqnorm(residuos)
#qqnorm dos preditores dos efeitos dos trabalhadores#
eblupsworker<-ranef(machines1r)$Worker
qqnorm(eblupsworker[,1])
#qqnorm dos preditores dos efeitos da interacção#
eblupsinteraccao<-ranef(machines1r)$`Worker:Machine`
qqnorm(eblupsinteraccao[,1])