valores = c(
17,14,13,12,
14,14,13,10,
12,13,12,9,
13,11,11,12,
11,12,10,8
)
trat = factor(rep(1:5,each=4))
bloc = factor(c(
2,3,4,5,
1,2,4,5,
1,3,4,5,
1,2,3,4,
1,2,3,5
))
data.frame(valores, trat, bloc)## valores trat bloc
## 1 17 1 2
## 2 14 1 3
## 3 13 1 4
## 4 12 1 5
## 5 14 2 1
## 6 14 2 2
## 7 13 2 4
## 8 10 2 5
## 9 12 3 1
## 10 13 3 3
## 11 12 3 4
## 12 9 3 5
## 13 13 4 1
## 14 11 4 2
## 15 11 4 3
## 16 12 4 4
## 17 11 5 1
## 18 12 5 2
## 19 10 5 3
## 20 8 5 5Para os dados apresentados, temos que o delineamento foi o Blocos balanceados Incompletos, no qual queremos analisar o efeito do combustível.
\[y_{ij}=\mu+\tau_i+\beta_j+e_{ij} \begin{cases} i=1,...,5\\ j=1,...,5\\ \end{cases}\]
library(dplyr)## Warning: pacote 'dplyr' foi compilado no R versão 4.4.2##
## Anexando pacote: 'dplyr'## Os seguintes objetos são mascarados por 'package:stats':
##
## filter, lag## Os seguintes objetos são mascarados por 'package:base':
##
## intersect, setdiff, setequal, unionk = 4
a = 5
r = 4
lambda = r*(k-1)/(a-1)
# ANOVA
mod = aov(valores~bloc+trat)
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 31.20 7.800 8.566 0.00216 **
## trat 4 35.73 8.933 9.810 0.00125 **
## Residuals 11 10.02 0.911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativa
## TAUS
bloc_falt = c(1, 3, 2, 5, 4)
tot_tau = tapply(valores, trat, sum)
tot_bet = tapply(valores, bloc, sum)
Q1 = tot_tau[1] - sum(tot_bet[-bloc_falt[1]])/k
Q2 = tot_tau[2] - sum(tot_bet[-bloc_falt[2]])/k
Q3 = tot_tau[3] - sum(tot_bet[-bloc_falt[3]])/k
Q4 = tot_tau[4] - sum(tot_bet[-bloc_falt[4]])/k
Q5 = tot_tau[5] - sum(tot_bet[-bloc_falt[5]])/k
Q = c(Q1,Q2,Q3,Q4,Q5)
sum(k*Q^2/(lambda*a))## [1] 35.73333taus_head = k*Q/(lambda*a)
## Betas
mu = mean(valores)
B1 = (mu*k - taus_head[1]- tot_bet[1]) /k
B2 = (mu*k - taus_head[3]- tot_bet[2]) /k
B3 = (mu*k - taus_head[2]- tot_bet[3]) /k
B4 = (mu*k - taus_head[5]- tot_bet[4]) /k
B5 = (mu*k - taus_head[4]- tot_bet[5]) /k
B = c(B1,B2,B3,B4,B5)
## Verificando
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 31.20 7.800 8.566 0.00216 **
## trat 4 35.73 8.933 9.810 0.00125 **
## Residuals 11 10.02 0.911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1data.frame(valores,
mod$fitted.values,
mu + B[bloc] + taus_head[trat],
bloc,
trat
)## valores mod.fitted.values mu...B.bloc....taus_head.trat. bloc trat
## 1 17 15.650000 12.850000 2 1
## 2 14 14.383333 14.116667 3 1
## 3 13 14.250000 14.250000 4 1
## 4 12 11.716667 16.783333 5 1
## 5 14 13.783333 11.783333 1 2
## 6 14 14.183333 11.383333 2 2
## 7 13 12.783333 12.783333 4 2
## 8 10 10.250000 15.316667 5 2
## 9 12 12.850000 10.850000 1 3
## 10 13 11.983333 11.716667 3 3
## 11 12 11.850000 11.850000 4 3
## 12 9 9.316667 14.383333 5 3
## 13 13 12.116667 10.116667 1 4
## 14 11 12.516667 9.716667 2 4
## 15 11 11.250000 10.983333 3 4
## 16 12 11.116667 11.116667 4 4
## 17 11 11.250000 9.250000 1 5
## 18 12 11.650000 8.850000 2 5
## 19 10 10.383333 10.116667 3 5
## 20 8 7.716667 12.783333 5 5## Com contraste seria algo
#
# c1 <- c(1,1,-1,-1)
# c2 <- c(1,-1,0,0)
# c3 <- c(0,0,1,-1)
#
#
#
# sqc1 <- k*sum(c1*c(q1,q2,q3,q4))^2/(2*4*sum(c1^2))
# sqc2 <- k*sum(c2*c(q1,q2,q3,q4))^2/(2*4*sum(c2^2))
# sqc3 <- k*sum(c3*c(q1,q2,q3,q4))^2/(2*4*sum(c3^2))
#
# sqc1+sqc2+sqc3
#
# qmres <- summary(mod1)[[1]][3,3]
# # um dos casos
# 1-pf(sqc1/qmres,1,5)|medbloc <- tapply(valores,bloc,mean)
media <- mean(valores)
SS_B = sum( (medbloc-media)^2)*4
medtrat <- tapply(valores, trat, mean)
SS_T = sum( (medtrat- media)^2)*4
tot_tau## 1 2 3 4 5
## 56 51 46 47 41tot_bet## 1 2 3 4 5
## 50 54 48 50 39