EDMARY SILVEIRA BARRETO ARAÚJO1
NORMA SOUZA DE OLIVEIRA2
UNIVERSIDADE FEDERAL DE LAVRAS - UFLA
DOUTORADO INTERINSTITUCIONAL (DINTER) EM ESTATÍSTICA
E EXPERIMENTAÇÃO AGROPECUÁRIA
5ª LISTA DE EXERCÍCIOS
SALVADOR 2010
1
Edmary Silveira Barreto Araújo: Mestre em Matemática, Doutoranda em Estatística e Experimentação
Agropecuária pela Universidade Federal de Lavras em convênio com o Instituto Federal da Bahia – IFBA,
Professora Assistente do IFBA. www.professoraedmary.com
2
Norma Souza de Oliveira: Licenciada em Matemática, Especialista em Educação, Mestre em Pedagogia
Profissional, aluna do Doutorado em Estatística e Experimentação Agropecuária pela Universidade Federal de
Lavras em convênio com o Instituto Federal da Bahia – IFBA, Professora de Matemática do IFBA.
1. Um experimento para avaliar a produção de cana-de-açúcar foi realizado no
delineamento em blocos casualizados com 9 tratamentos e 3 repetições. Os tratamentos
constituem um fatorial 3x3, sendo 3 níveis de fósforo (X1 = 0,50 e 100 kg/ha de P2O5) e
3 níveis de potássio (X2 = 0,50 e 100 kg/ha de K2O). Os resultados da produção de
cana-de-açucar, em t/ha, estão apresentados na tabela 1.
Tabela 1: Dados de produção de cana-de-açucar, em t/ha
Fósforo (X1)
0
0
0
50
50
50
100
100
100
Potássio (X2)
0
50
100
0
50
100
0
50
100
1
57
60
64
66
71
80
72
80
87
Blocos
2
49
52
48
70
76
72
78
91
94
3
54
62
65
69
77
84
80
85
96
Considerando o modelo
Yi = b 0 + b1 X 1i + b 2 X 2i + b 3 X 12i + b 4 X 22i + b 5 X 1i X 2i + ei ,
em que ei ~ N (0,s 2 ) , pede-se:
a) Obter a equação de regressão ajustada;
^
Yi = 53.5741 + 0.3544 X 1i + 0.1067 X 2i - 0.0012 X 12i - 0.0005 X 22i + 0.0010 X 1i X 2i
Usando o software R.
library(fBasics)
Y<matrix(c(57,60,64,66,71,80,72,80,87,49,52,48,70,76,72,78,91,94,54,62,65,69,77,84,80,85,96),27)
X<-matrix(c(rep(1,27),rep(0,3),rep(50,3),rep(100,3),rep(0,3),rep(50,3),rep(100,3),rep(0,3),rep(50,3),
rep(100,3),0,50,100,0,50,100,0,50,100,0,50,100,0,50,100,0,50,100,
0,50,100,0,50,100,0,50,100,rep(0,3),rep(50^2,3),rep(100^2,3),rep(0,3),rep(50^2,3),rep(100^2,3),
rep(0,3),rep(50^2,3),rep(100^2,3),0,(50^2),(100^2),0,(50^2),(100^2),0,(50^2),(100^2),
0,(50^2),(100^2),0,(50^2),(100^2),0,(50^2),(100^2),0,(50^2),(100^2),0,(50^2),(100^2),0,(50^2),
(100^2),rep(0,4),(50^2),50*100,0,100*50,(100^2),rep(0,4),(50^2),50*100,0,100*50,(100^2),
rep(0,4),(50^2),50*100,0,100*50,(100^2)),27)
rk(X) # 6
A<-t(X)%*%X
B<-t(X)%*%Y
rk(A) # 6
round(solve(A),2)
Betas<- solve(A)%*%t(X)%*%Y
round(Betas,4)
[1,] 53.5741
[2,] 0.3544
[3,] 0.1067
[4,] -0.0012
[5,] -0.0005
[6,] 0.0010
#
#
#
#
#
#
Beta0
Beta1
Beta2
Beta3
Beta4
Beta5
b) Efetuar a análise de variância da regressão testando a hipótese
H 0 : b1 = b 2 = b 3 = b 4 = b 5 = 0;
OBS.: Observe que têm-se 8 graus de liberdade para Tratamentos e o modelo de
regressão proposto tem p + 1 = 6 parâmetros (5 graus de liberdade para
Regressão). Portanto, têm-se 3 graus de liberdade para a falta de ajustamento.
reg<-lm(Y~X)
anova(reg)
Analysis of Variance Table
Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
X
5 4172.3 834.5 33.972 2.228e-09 ***
Residuals 21 515.8
24.6
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Rejeita H 0 : b1 = b 2 = b 3 = b 4 = b 5 = 0 , ao nível de 1%; pelo menos um bi ¹ 0.
Observe:
Analysis of Variance Table
Response: Y
Df
Trat
8
Bloco
2
Residuals 16
--Signif. codes:
Sum Sq Mean Sq
4173.4
521.7
112.5
56.3
402.1
25.1
F value Pr(>F)
20.7556 4.604e-07 ***
2.2383
0.1390
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Análise de variância da regressão para a falta de ajustamento
F.V
G.L
S.Q
Q.M
Regressão
5
4172.3
834.46
Falta de ajuste
3
1.1
0.367
Tratamentos
(8)
(4173.4)
Blocos
2
112.5
Resíduo
16
402.1
25.1
Total
26
4688
F
33.24
0.014
c) Calcular o coeficiente de determinação múltipla, interpretando-o corretamente.
R2 = SQReg/SQTotal = 0.8899767
89% da variabilidade da produção de cana-de-açúcar, em t/ha, é explicada pelo
modelo considerado.
^
d) Apresentar as estimativas dos coeficientes de regressão ( b i , com i = 1,2,3,4,5 ) e
os respectivos intervalos de confiança, considerando o nível de confiança de
95%;
^
b 0 = 53.5741
^
^
^
^
^
b1 = 0.3544 b 2 = 0.1067 b 3 = -0.0012 b 4 = -0.0005 b 5 = 0.0010
Usando o software R.
X1<-X[,2]
X2<-X[,3]
X1X1<-X[,4]
X2X2<-X[,5]
X1X2<-X[,6]
rega<-lm(Y~X1+X2+X1X1+X2X2+X1X2)
confint(rega)
(Intercept)
X1
X2
X1X1
X2X2
X1X2
2.5 %
48.2333
0.1694
-0.0783
-0.0029
-0.0022
-0.0002
97.5 %
58.9149
0.5395
0.2917
0.0004
0.0012
0.0022
e) Verificar se os coeficientes de regressão são significativos ao nível de 5% de
probabilidade pelo teste t, isto é, testar as hipóteses H 0 : b i = 0 vs H a : b i ¹ 0 ;
summary(rega)
Call:
lm(formula = Y ~ X1 + X2 + X1X1 + X2X2 + X1X2)
Residuals:
Min
1Q Median
3Q Max
-11.130 -4.157 1.259 3.426 5.870
Coefficients:
Estimat Std. Error t value Pr(>|t|)
(Intercept) 53.5740741 2.5681736 20.861 1.61e-15 ***
X1
0.3544444 0.0889641 3.984 0.000675 ***
X2
0.1066667 0.0889641 1.199 0.243889
X1X1
-0.0012444 0.0008093 -1.538 0.139069
X2X2
-0.0005111 0.0008093 -0.632 0.534509
X1X2
0.0010000 0.0005723 1.747 0.095176 .
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.956 on 21 degrees of freedom
Multiple R-squared: 0.89,
Adjusted R-squared: 0.8638
F-statistic: 33.97 on 5 and 21 DF, p-value: 2.228e-09
^
^
b 0 = 53.5741 e b1 = 0.3544 são significativos ao nível de 5%. Rejeita H0.
^
^
^
^
b 2 = 0.1067 , b 3 = -0.0012 , b 4 = -0.0005 e b 5 = 0.0010 não são significativos ao
nível de 5%. Não rejeita (Aceita) H0.
f) Determinar o ponto crítico, estudar sua natureza e calcular a resposta estimada
neste ponto.
Yi = 53.5741 + 0.3544 X 1i + 0.1067 X 2i - 0.0012 X 12i - 0.0005 X 22i + 0.0010 X 1i X 2i
dYi = 0.3544 dX 1i - 0.0024 X 1i dX 1i + 0.0010 X 2i dX 1i + 0.1067dX 2i - 0.001X 2i dX 2i + 0.0010 X 1i dX 2i = 0
dYi
= 0.3544 - 0.0024 X 1i + 0.0010 X 2i = 0
dX 1i
dYi
= 0.1067 + 0.001X 1i - 0.0010 X 2i = 0
dX 2i
^
^
P = ( X 1i , X 2i ) = ( 329.36, 436.06)
d 2Yi = -0.0024dX 12i + 0.0010dX 1i dX 2i + 0.0010dX 2 i dX 1i - 0.001dX 22i
d 2Yi = [dX 1i
é- 0.0024 0.001 ù é dX 1i ù
dX 2i ]ê
- 0.001úû êëdX 2i úû
ë 0.001
Forma Quadrática.
d<-matrix(c(-0.0024,0.001,0.0010,-0.001),2)
eigen(d) # l1 = -0.0005 l2 = -0.003 (negativos)
P é ponto máximo.
Resposta estimada neste ponto:
X1=329.36
X2=436.06
Yhat=53.57+(0.36*X1)+(0.1066666667*X2)-(0.001)*(X1^2)-(0.0005)*(X2^2)+(0.001)*X1*X2
Yhat = 158.7212
2. Um experimento foi conduzido utilizando-se o delineamento inteiramente
casualizado, obtendo-se os dados apresentados na tabela 2.
Tabela 2: Resultados experimentais
Tratamentos
1
2
5
6
4
8
3
7
4
9
4
8
3
9
8
10
8
11
a) Representar o modelo matricialmente, na forma Y = Xb + e ;
é e11 ù
é 5 ù é1 1 0 0ù
êe ú
ê 4 ú ê1 1 0 0ú
ê 12 ú
ê ú ê
ú
ê e13 ú
ê 3 ú ê1 1 0 0ú
ê ú
ê ú ê
ú
ê e14 ú
ê 4 ú ê1 1 0 0ú
ê e15 ú
ê 4 ú ê1 1 0 0ú
ê ú
ê ú ê
ú
êe21 ú
ê 6 ú ê1 0 1 0ú
m
é
ù
êe ú
ê 8 ú ê1 0 1 0ú
ê ú ê
ú êt 1 ú ê 22 ú
ê 7 ú = ê1 0 1 0ú.ê ú + êe23 ú
ê 9 ú ê1 0 1 0ú êt 2 ú êe ú
ê ú ê
ú êët 3 úû ê 24 ú
ê 8 ú ê1 0 1 0ú { êe25 ú
ê ú
ê ú ê
ú b
9
1
0
0
1
ê e31 ú
ê ú ê
ú
êe32 ú
ê 8 ú ê1 0 0 1 ú
ê ú
ê ú ê
ú
êe33 ú
ê10ú ê1 0 0 1 ú
êe ú
ê 8 ú ê1 0 0 1 ú
ê 34 ú
ê ú ê
ú
êë11úû êë1 0 0 1 úû
êëe35 úû
{ 14
{
4244
3
Y
X
e
b) Apresentar o sistema de equações normais (S.E.N.): X ´ Xb 0 = X ´Y ;
Usando o software R.
Y<-matrix(c(5,4,3,4,4,6,8,7,9,8,9,8,10,8,11),15)
X<-matrix(c(rep(1,20),rep(0,15),rep(1,5),rep(0,15),rep(1,5)),15)
A<-t(X)%*%X
B<-t(X)%*%Y
ì15m 0 + 5t 10 + 5t 02 + 5t 03 = 104
ï 0
0
ï 5m + 5t 1 = 20
S.E.N : A.b = B. í 0
0
ï 5m + 5t 2 = 38
ï 5m 0 + 5t 0 = 46
3
î
c) Obter duas soluções diferentes para o S.E.N. ;
Usando o software R.
library(fBasics)
rk(A)
inv<-ginv(A) # A+
Am<-matrix(c(rep(0,5),1/5,rep(0,4),1/5,rep(0,4),1/5),4) A-
é 5,2 ù
ê- 1,2ú
0
ú
Betas1<- inv%*%t(X)%*%Y # b1 = ê
ê 2,4 ú
ê
ú
ë 4,0 û
é0,0ù
ê4,0ú
0
Betas2<-Am%*%t(X)%*%Y # b 2 = ê ú
ê 7 ,6 ú
ê ú
ë9,2û
d) Efetuar a análise de variância, testando a hipótese
considerando a = 5%.
H 0 : m1 = m 2 = m3 ,
Usando o software R.
library(fBasics)
dados<-data.frame(Y=c(5,4,3,4,4,6,8,7,9,8,9,8,10,8,11), Trat=as.factor(c(rep(1:3,each=5))))
n<-length(Y)
X1<-X[,1]
I<- diag(n)
P<- X %*% ginv(t(X)%*%X) %*% t(X)
rk(P)
P1<- X1 %*% solve(t(X1)%*%X1) %*% t(X1)
P.to<- I - P1
P.trat<- P - P1
P.erro<- I – P
SQTotal<- t(Y) %*% P.to %*% Y # 84.93333
SQTrat<- t(Y) %*% P.trat %*% Y # 70.93333
SQErro<- t(Y) %*% P.erro %*% YSQTotal # 14
GLTo<- rk(P.to) # 14
GLTrat<- rk(P.trat) # 2
GLErro<- rk(P.erro) # 12
F.calc<- (SQTrat/GLTrat)/(SQErro/GLErro) # 30.4
1-pf(F.calc,GLTrat,GLErro) # 2.005853e-05
Rejeita H 0 : m1 = m 2 = m3 ao nível de 5%.
Outra maneira:
anova(resp)
Analysis of Variance Table
Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
Trat
2 70.933 35.467
30.4
2.006e-05 ***
Residuals 12 14.000 1.167
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Portanto, rejeita H 0 : m1 = m + t 1 = m 2 = m + t 2 = m3 = m + t 3 Þ H 0 : t 1 = t 2 = t 3 ao
nível de 5%. Há evidências de que ao menos duas médias de tratamentos diferem
entre si.
e) Quais dentre as funções paramétricas dadas a seguir são estimáveis no modelo
considerado? Utilize pelo menos dois critérios, justificando sua resposta.
Modelo: yij = m + t i + eij
básica estimável.
i = 1,2,3
j = 1,2,3,4,5. E ( yij ) = m + t i é a função
1º Critério: Teorema: Uma condição necessária e suficiente para que a função
l´b seja estimável no modelo Y = Xb + e , Gauss Markov, é que l Î C ( X ´).
~
a) Posto [X´: l ] = posto [X´].
b) Posto [X´X : l ] = posto [X´X].
2º Critério: l´1 b , l´2 b ,..., l´p b
~
~
p funções lineares paramétricas estimáveis
~
p
Þ l´b = å l´i b é estimável.
~
i =1
~
I
I
i =1
i =1
3º Critério: Seja l ´b = c0 m + å cit i . c0 = å ci Þ l´b é estimável.
~
~
i) m + t 1 ;
3º Critério:
m + t 1 é uma função básica estimável com c0 = 1 = c1 = 1.
1º Critério:
l´ = [1,1,0,0]
l<-matrix(c(1,1,0,0),4) # lambda
rk(t(X)) # 3
rk(cbind(t(X),l)) # 3
rk(A) # posto( A = X´X) # 3
rk(cbind(A,l)) # posto [A = X´X: lambda] # 3
l Î C ( X ´).
Portanto, m + t 1 é estimável
ii) m + 2t 1 ;
3º Critério:
c0 = 1 ¹ c1 = 1
1º Critério:
l´ = [1,2,0,0]
l<-matrix(c(1,2,0,0),4) # lambda
rk(t(X)) # 3
rk(cbind(t(X),l)) # 4
rk(A) # posto( A = X´X) # 3
rk(cbind(A,l)) # posto [A = X´X: lambda] # 4
l Ï C (X ´).
Portanto, m + 2t 1 não é estimável.
iii) 2m + t 1 + t 2 ;
2º Critério:
2 m + t 1 + t 2 = (m + t 1 ) + (m + t 2 ) (combinação linear de funções estimáveis).
123 1
424
3
é estimável
é estimável
1º Critério:
l´ = [2,1,1,0]
l<-matrix(c(2,1,1,0),4) # lambda
rk(t(X)) # 3
rk(cbind(t(X),l)) # 3
rk(A) # posto( A = X´X) # 3
rk(cbind(A,l)) # posto [A = X´X: lambda] # 3
l Î C ( X ´).
Portanto, 2m + t 1 + t 2 é estimável.
iv) t 1 + t 2 - 2t 3 ;
1º Critério:
l´ = [0,1,1,-2]
l<-matrix(c(0,1,1,-2),4) # lambda
rk(t(X)) # 3
rk(cbind(t(X),l)) # 3
rk(A) # posto( A = X´X) # 3
rk(cbind(A,l)) # posto [A = X´X: lambda] # 3
l Î C (X ´).
3º Critério:
c0 = 0 = 1 + 1 - 2 = c1 + c2 + c3
Portanto, t 1 + t 2 - 2t 3 é estimável.
v) t 1 - t 2 ;
2º
t1 - t 2 = m + t1 - ( m + t 2 )
123
123
Critério:
é estimável
(combinação
linear
de
funções
e
para
é estimável
estimáveis).
3º Critério:
c0 = 0 = 1 - 1 = c1 + c2.
Portanto, t 1 - t 2 é estimável.
vi) t 1 + t 2 - t 3 .
1º Critério:
l´ = [0,1,1,-1]
l<-matrix(c(0,1,1,-1),4) # lambda
rk(t(X)) # 3
rk(cbind(t(X),l)) # 4
rk(A) # posto( A = X´X) # 3
rk(cbind(A,l)) # posto [A = X´X: lambda] # 4
l Ï C (X ´).
3º Critério:
c0 = 0 ¹ 1 + 1 – 1 = 1 = c1 + c2 + c3.
Portanto, t 1 + t 2 - t 3 não é estimável.
ém ù
ê ú
é0 1 - 1 0 ù êt 1 ú é0ù
f) Considerando
H 0 : B´b = f Þ ê
=ê ú
ú.
ë0 1 1 - 2û êt 2 ú ë0û
ê ú
ët 3 û
( X ´ X ) simétrica, calcule:
^
^
SQH 0 = ( B´b - B´b )´[B´( X ´ X ) - B ]-1 ( B´b - B´b )
Sabendo-se que
SQH 0
posto[ B]
~ F ( posto[ B]; n - posto[ X ])
SQErro
n - posto[ X ]
Aplique o teste F, concluindo corretamente.
No software R, temos:
B<- matrix(c(0,1,-1,0,0,1,1,-2),4)
BlBhat<-t(B)%*%Betas2
SQH<-t(BlBhat)%*%solve(t(B)%*%Am%*%B)%*%BlBhat ## 70.9333
^
^
SQH 0 = ( B´b - B´b )´[B´( X ´ X ) - B ]-1 ( B´b - B´b ) = 70.93333
Fc<-(SQH/rk(B))/(SQErro/n-rk(X)) ## 30.40
pvalue<- 1-pf(Fc,rk(B),n-rk(X)) # 2.006326e-05
ì t1 = t 2
H0 : í
ît 1 + t 2 = 2t 3
Rejeita H0 ao nível de 5%.
g) Sejam as duas funções paramétricas (contrastes);
l´1b = t 1 - t 2 = m1 - m2 ;
l´2 b = t 1 + t 2 - 2t 3 = m1 + m2 - 2m3
Teste as hipóteses a seguir pelo teste F, concluindo corretamente. Utilize
a = 5%.
i. H 0 : l´1b = 0; ( Tratamento1) vs (Tratamento2)
No software R.
Cr1<-matrix(c(0,1,-1,0),1, byrow=T)
Theta1=0
SQCr1=t(Cr1%*%Betas2-Theta1) %*% solve(Cr1%*%Am%*%t(Cr1)) %*% (Cr1%*%Betas2-Theta1)
SQCr1 # 32.4
FcalCr1=SQCr1/(rk(Cr1)*(SQErro/(n-rk(X)))) # 27.77143
1-pf(FcalCr1,rk(Cr1),n-rk(X)) # 0.0001978544
Rejeita H0 ao nível de 5%. Em média existe diferença entre os Trat 1 e 2.
ii. H 0 : l´2 b = 0. (Tratamento 1 e Tratamento 2) vs ( Tratamento 3)
Cr2<-matrix(c(0,1,1,-2),1, byrow=T)
Theta2=0
SQCr2=t(Cr2%*%Betas2-Theta2) %*% solve(Cr2%*%Am%*%t(Cr2)) %*% (Cr2%*%Betas2Theta2) # 38.5333
Fcalc2=SQCr2/(rk(Cr2)*(SQErro/(n-rk(X)))) # 33.02857
1-pf(Fcalc2,rk(Cr2),n-rk(X)) # 9.20592e-05
Rejeita H0 ao nível de 5%. Existe diferença entre (Trat 1 e 2) vs (Trat 3).