\[\\[1mm]\]
Kůň by měl v modelu být jako faktorová proměnná, protože nedává smysl, aby tam byl jako číselná proměnná.
Opravená analýza (stejné hodnoty parametrů a generování odezvy jako minule):
| horse | technique | position | time | response |
|---|---|---|---|---|
| 1 | 0 | 0 | 0 | 62.57737 |
| 1 | 0 | 1 | 1 | 67.07907 |
| 2 | 0 | 0 | 1 | 61.68652 |
| 2 | 0 | 1 | 0 | 68.89723 |
| 3 | 1 | 0 | 0 | 48.13785 |
| 3 | 1 | 1 | 1 | 60.05408 |
| 4 | 1 | 0 | 1 | 56.93477 |
| 4 | 1 | 1 | 0 | 57.27467 |
| horse | technique | position | time | response |
|---|---|---|---|---|
| 1 | 0 | 0 | 0 | 60.75641 |
| 1 | 1 | 1 | 1 | 57.82523 |
| 2 | 1 | 0 | 0 | 48.84914 |
| 2 | 0 | 1 | 1 | 70.22613 |
| 3 | 0 | 0 | 1 | 68.05427 |
| 3 | 1 | 1 | 0 | 54.17742 |
| 4 | 1 | 0 | 1 | 57.12341 |
| 4 | 0 | 1 | 0 | 61.94710 |
Modely:
lm5.30 <- lm(response ~ technique + position + time + factor(horse),data=data5.30)
lm5.31 <- lm(response ~ technique + position + time + factor(horse),data=data5.31)Výstupy alias funkce:
alias(lm5.30,partial=TRUE)## Model :
## response ~ technique + position + time + factor(horse)
##
## Complete :
## (Intercept) technique position time factor(horse)2
## factor(horse)4 0 1 0 0 0
## factor(horse)3
## factor(horse)4 -1
##
## Partial :
## (Intercept) technique position time factor(horse)2
## (Intercept) 0 -0.5773503 -0.4082483 -0.4082483 -0.5773503
## technique 0 0.0000000 0.0000000 0.0000000 0.5000000
## position 0 0.0000000 0.0000000 0.0000000 0.0000000
## time 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)2 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)3 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)3
## (Intercept) 0.0
## technique -0.5
## position 0.0
## time 0.0
## factor(horse)2 0.0
## factor(horse)3 0.0alias(lm5.31,partial=TRUE)## Model :
## response ~ technique + position + time + factor(horse)
##
## Partial :
## (Intercept) technique position time factor(horse)2
## (Intercept) 0 -0.3779645 -0.3779645 -0.3779645 -0.5345225
## technique 0 0.0000000 0.0000000 0.0000000 0.0000000
## position 0 0.0000000 0.0000000 0.0000000 0.0000000
## time 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)2 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)3 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)4 0 0.0000000 0.0000000 0.0000000 0.0000000
## factor(horse)3 factor(horse)4
## (Intercept) -0.5345225 -0.5345225
## technique 0.0000000 0.0000000
## position 0.0000000 0.0000000
## time 0.0000000 0.0000000
## factor(horse)2 0.5000000 0.5000000
## factor(horse)3 0.0000000 0.5000000
## factor(horse)4 0.0000000 0.0000000Vidíme, že ve struktuře 5.30, kde každý kůň dostal na obě kopyta pouze jeden způsob léčení, se efekt posledního koně nedá spočítat, protože je lineárně závislý na způsobu léčení a efektu jiného koně. Ve struktuře 5.31 tomu tak není a jdou odhadnout efekty všech koní, lze mezi nimi ale pozorovat multikolinearitu.
\[\\[2mm]\]
Načítám data:
sausage <- factor(c(rep(c(1:12),4)))
judge <- factor(c(rep("A",12),rep("B",12),rep("C",12),rep("D",12)))
salt <- c(2.4,2.3,2.6,6.0,4.6,3.8,4.6,3.6,0.5,2.2,4.3,9.0,
2.5,2.0,2.5,5.5,4.1,2.3,4.5,3.1,0.3,1.8,4.5,7.0,
3.0,3.0,3.0,6.0,3.9,4.2,5.0,3.9,1.1,2.5,4.7,7.9,
2.7,2.8,2.6,5.8,4.9,4.4,4.9,4.0,1.3,2.3,4.6,8.5)
klobasky <- data.frame(sausage,judge,salt)Výpočet tabulky 2:
model <- lm(salt ~ sausage*judge, data = klobasky, contrasts = list(sausage="contr.sum",judge="contr.sum"))
anova(model)## Warning in anova.lm(model): ANOVA F-tests on an essentially perfect fit are
## unreliable## Analysis of Variance Table
##
## Response: salt
## Df Sum Sq Mean Sq F value Pr(>F)
## sausage 11 160.637 14.6034
## judge 3 3.938 1.3125
## sausage:judge 33 4.457 0.1351
## Residuals 0 0.000Model 5.2.1:
lme5.2.1 <- lme(salt~judge*sausage,random=~1|sausage,data=klobasky)## Error in MEEM(object, conLin, control$niterEM): Singularity in backsolve at level 0, block 1lmer5.2.1 <- lmer(salt~judge*sausage+(1|sausage),data=klobasky)## Error in is.nloptr(ret): objective in x0 returns NAR očividně nedokáže odhadnout parametry modelu 5.2.1 (nebo něco dělám špatně).
Testování 5.3:
z <- matrix(0,nrow=4,ncol=12)
judges <- c("A","B","C","D")
sausages <- c(1:12)
for (i in 1:4) {
for (j in 1:12) {
z[i,j] <- (klobasky[klobasky$sausage==sausages[j] & klobasky$judge==judges[i],]$salt-
mean(klobasky[klobasky$sausage==sausages[j],]$salt) -
mean(klobasky[klobasky$judge==judges[i],]$salt) +
mean(klobasky$salt))
}
}
ezzt <- eigen(z %*% t(z))
eztz <- eigen(t(z) %*% z)
p <- min((length(judges)-1),(length(sausages)-1))
n <- max((length(judges)-1),(length(sausages)-1))
I <- eztz$values[1:p]
U1 <- I[1]/sum(I)
U1## [1] 0.6041585
testová statistika \(U_1\) = 0.6041585 není signifikantní na 10% levelu.
\[\\[2mm]\]
Netuším, jak se dělají half-normal plots v Rku (google nepomohl).
\[\\[5mm]\]