Especificação, estimação e parametrização
##-----------------------------------------------------------------------------
## Pacotes.
pkg <- c("lattice", "latticeExtra", "plyr", "reshape")
sapply(pkg, require, character.only=TRUE)
##-----------------------------------------------------------------------------
## Matrizes de delineamentos.
incid <- function(vec){
stopifnot(is.factor(vec))
sapply(levels(vec), function(i) as.integer(vec==i))
}
## DIC.
da <- data.frame(trat=gl(3, 4, labels=c("A","B","C")))
incid(da$trat)
## A B C
## [1,] 1 0 0
## [2,] 1 0 0
## [3,] 1 0 0
## [4,] 1 0 0
## [5,] 0 1 0
## [6,] 0 1 0
## [7,] 0 1 0
## [8,] 0 1 0
## [9,] 0 0 1
## [10,] 0 0 1
## [11,] 0 0 1
## [12,] 0 0 1
## DBC.
da <- expand.grid(bloc=gl(3,1, labels=c("I","II","III")),
trat=gl(4,1, labels=LETTERS[1:4]))
incid(da$bloc)
## I II III
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
## [4,] 1 0 0
## [5,] 0 1 0
## [6,] 0 0 1
## [7,] 1 0 0
## [8,] 0 1 0
## [9,] 0 0 1
## [10,] 1 0 0
## [11,] 0 1 0
## [12,] 0 0 1
incid(da$trat)
## A B C D
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1
##-----------------------------------------------------------------------------
## Estimação.
url <- "http://www.leg.ufpr.br/~walmes/data/pimentel_batatinha.txt"
da <- read.table(url, header=TRUE, sep="\t")
str(da)
## 'data.frame': 32 obs. of 3 variables:
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ variedade: Factor w/ 8 levels "B 116-51","B 1-52",..: 7 6 8 5 3 2 1 4 7 6 ...
## $ producao : num 9.2 21.1 22.6 15.4 12.7 20 23.1 18 13.4 27 ...
names(da) <- substr(names(da), 1, 4)
str(da)
## 'data.frame': 32 obs. of 3 variables:
## $ bloc: Factor w/ 4 levels "I","II","III",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ vari: Factor w/ 8 levels "B 116-51","B 1-52",..: 7 6 8 5 3 2 1 4 7 6 ...
## $ prod: num 9.2 21.1 22.6 15.4 12.7 20 23.1 18 13.4 27 ...
da <- droplevels(subset(da, da$vari%in%levels(da$vari)[1:3]))
da <- arrange(da, vari, bloc)
str(da)
## 'data.frame': 12 obs. of 3 variables:
## $ bloc: Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ vari: Factor w/ 3 levels "B 116-51","B 1-52",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ prod: num 23.1 24.2 26.4 16.3 20 21.1 20 28 12.7 18 ...
X <- cbind(1, incid(da$bloc), incid(da$vari))
X
## I II III IV B 116-51 B 1-52 B 25-50 E
## [1,] 1 1 0 0 0 1 0 0
## [2,] 1 0 1 0 0 1 0 0
## [3,] 1 0 0 1 0 1 0 0
## [4,] 1 0 0 0 1 1 0 0
## [5,] 1 1 0 0 0 0 1 0
## [6,] 1 0 1 0 0 0 1 0
## [7,] 1 0 0 1 0 0 1 0
## [8,] 1 0 0 0 1 0 1 0
## [9,] 1 1 0 0 0 0 0 1
## [10,] 1 0 1 0 0 0 0 1
## [11,] 1 0 0 1 0 0 0 1
## [12,] 1 0 0 0 1 0 0 1
t(X)%*%X
## I II III IV B 116-51 B 1-52 B 25-50 E
## 12 3 3 3 3 4 4 4
## I 3 3 0 0 0 1 1 1
## II 3 0 3 0 0 1 1 1
## III 3 0 0 3 0 1 1 1
## IV 3 0 0 0 3 1 1 1
## B 116-51 4 1 1 1 1 4 0 0
## B 1-52 4 1 1 1 1 0 4 0
## B 25-50 E 4 1 1 1 1 0 0 4
crossprod(X, X)
## I II III IV B 116-51 B 1-52 B 25-50 E
## 12 3 3 3 3 4 4 4
## I 3 3 0 0 0 1 1 1
## II 3 0 3 0 0 1 1 1
## III 3 0 0 3 0 1 1 1
## IV 3 0 0 0 3 1 1 1
## B 116-51 4 1 1 1 1 4 0 0
## B 1-52 4 1 1 1 1 0 4 0
## B 25-50 E 4 1 1 1 1 0 0 4
crossprod(X) ## É possível interpretar tais números?
## I II III IV B 116-51 B 1-52 B 25-50 E
## 12 3 3 3 3 4 4 4
## I 3 3 0 0 0 1 1 1
## II 3 0 3 0 0 1 1 1
## III 3 0 0 3 0 1 1 1
## IV 3 0 0 0 3 1 1 1
## B 116-51 4 1 1 1 1 4 0 0
## B 1-52 4 1 1 1 1 0 4 0
## B 25-50 E 4 1 1 1 1 0 0 4
XlX <- crossprod(X)
## solve(XlX)
require(MASS)
## Loading required package: MASS
XlXi <- ginv(XlX) ## Inversa generalizada.
XlXi
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0.03324 0.00831 0.00831 0.00831 0.00831 0.01108 0.01108 0.01108
## [2,] 0.00831 0.25208 -0.08126 -0.08126 -0.08126 0.00277 0.00277 0.00277
## [3,] 0.00831 -0.08126 0.25208 -0.08126 -0.08126 0.00277 0.00277 0.00277
## [4,] 0.00831 -0.08126 -0.08126 0.25208 -0.08126 0.00277 0.00277 0.00277
## [5,] 0.00831 -0.08126 -0.08126 -0.08126 0.25208 0.00277 0.00277 0.00277
## [6,] 0.01108 0.00277 0.00277 0.00277 0.00277 0.17036 -0.07964 -0.07964
## [7,] 0.01108 0.00277 0.00277 0.00277 0.00277 -0.07964 0.17036 -0.07964
## [8,] 0.01108 0.00277 0.00277 0.00277 0.00277 -0.07964 -0.07964 0.17036
beta <- XlXi%*%crossprod(X, da$prod)
##-----------------------------------------------------------------------------
## Restrições paramétricas ou restrições na solução.
## Zera nível 1 nível.
model.matrix(~bloc+vari, data=da)
## (Intercept) blocII blocIII blocIV variB 1-52 variB 25-50 E
## 1 1 0 0 0 0 0
## 2 1 1 0 0 0 0
## 3 1 0 1 0 0 0
## 4 1 0 0 1 0 0
## 5 1 0 0 0 1 0
## 6 1 1 0 0 1 0
## 7 1 0 1 0 1 0
## 8 1 0 0 1 1 0
## 9 1 0 0 0 0 1
## 10 1 1 0 0 0 1
## 11 1 0 1 0 0 1
## 12 1 0 0 1 0 1
## attr(,"assign")
## [1] 0 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$bloc
## [1] "contr.treatment"
##
## attr(,"contrasts")$vari
## [1] "contr.treatment"
## Zera efeito do nível k (último).
model.matrix(~bloc+vari, data=da,
contrasts.arg=list(bloc="contr.SAS"))
## (Intercept) blocI blocII blocIII variB 1-52 variB 25-50 E
## 1 1 1 0 0 0 0
## 2 1 0 1 0 0 0
## 3 1 0 0 1 0 0
## 4 1 0 0 0 0 0
## 5 1 1 0 0 1 0
## 6 1 0 1 0 1 0
## 7 1 0 0 1 1 0
## 8 1 0 0 0 1 0
## 9 1 1 0 0 0 1
## 10 1 0 1 0 0 1
## 11 1 0 0 1 0 1
## 12 1 0 0 0 0 1
## attr(,"assign")
## [1] 0 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$bloc
## [1] "contr.SAS"
##
## attr(,"contrasts")$vari
## [1] "contr.treatment"
## A soma dos efeito dá zero, com isso o efeito k é o nagativo da soma
## dos demais.
model.matrix(~bloc+vari, data=da,
contrasts.arg=list(bloc="contr.sum"))
## (Intercept) bloc1 bloc2 bloc3 variB 1-52 variB 25-50 E
## 1 1 1 0 0 0 0
## 2 1 0 1 0 0 0
## 3 1 0 0 1 0 0
## 4 1 -1 -1 -1 0 0
## 5 1 1 0 0 1 0
## 6 1 0 1 0 1 0
## 7 1 0 0 1 1 0
## 8 1 -1 -1 -1 1 0
## 9 1 1 0 0 0 1
## 10 1 0 1 0 0 1
## 11 1 0 0 1 0 1
## 12 1 -1 -1 -1 0 1
## attr(,"assign")
## [1] 0 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$bloc
## [1] "contr.sum"
##
## attr(,"contrasts")$vari
## [1] "contr.treatment"
model.matrix(~bloc+vari, data=da,
contrasts.arg=list(bloc="contr.helmert"))
## (Intercept) bloc1 bloc2 bloc3 variB 1-52 variB 25-50 E
## 1 1 -1 -1 -1 0 0
## 2 1 1 -1 -1 0 0
## 3 1 0 2 -1 0 0
## 4 1 0 0 3 0 0
## 5 1 -1 -1 -1 1 0
## 6 1 1 -1 -1 1 0
## 7 1 0 2 -1 1 0
## 8 1 0 0 3 1 0
## 9 1 -1 -1 -1 0 1
## 10 1 1 -1 -1 0 1
## 11 1 0 2 -1 0 1
## 12 1 0 0 3 0 1
## attr(,"assign")
## [1] 0 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$bloc
## [1] "contr.helmert"
##
## attr(,"contrasts")$vari
## [1] "contr.treatment"
##-----------------------------------------------------------------------------
## Como gerar a matriz do modelo.
## Matriz de incidência "pura". Modelo de parâmetros por cela ou sem
## restrição.
X <- incid(da$vari)
X <- model.matrix(~0+vari, data=da)
X
## variB 116-51 variB 1-52 variB 25-50 E
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 0 1 0
## 6 0 1 0
## 7 0 1 0
## 8 0 1 0
## 9 0 0 1
## 10 0 0 1
## 11 0 0 1
## 12 0 0 1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$vari
## [1] "contr.treatment"
## Monte a matriz de identifique os contrastes.
ctr <- C(da$vari, contr=contr.treatment)
ctr <- cbind(1, attr(ctr, "contrasts"))
ctr
## 2 3
## B 116-51 1 0 0
## B 1-52 1 1 0
## B 25-50 E 1 0 1
## Multiplique.
X%*%ctr
## 2 3
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 1 0
## 6 1 1 0
## 7 1 1 0
## 8 1 1 0
## 9 1 0 1
## 10 1 0 1
## 11 1 0 1
## 12 1 0 1
model.matrix(~vari, data=da,
contrasts.arg=list(vari=contr.treatment))
## (Intercept) vari2 vari3
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 1 0
## 6 1 1 0
## 7 1 1 0
## 8 1 1 0
## 9 1 0 1
## 10 1 0 1
## 11 1 0 1
## 12 1 0 1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$vari
## 2 3
## B 116-51 0 0
## B 1-52 1 0
## B 25-50 E 0 1
## Outros tipos de contraste.
C(da$vari, contr=contr.SAS)
## [1] B 116-51 B 116-51 B 116-51 B 116-51 B 1-52 B 1-52 B 1-52 B 1-52
## [9] B 25-50 E B 25-50 E B 25-50 E B 25-50 E
## attr(,"contrasts")
## 1 2
## B 116-51 1 0
## B 1-52 0 1
## B 25-50 E 0 0
## Levels: B 116-51 B 1-52 B 25-50 E
C(da$vari, contr=contr.sum)
## [1] B 116-51 B 116-51 B 116-51 B 116-51 B 1-52 B 1-52 B 1-52 B 1-52
## [9] B 25-50 E B 25-50 E B 25-50 E B 25-50 E
## attr(,"contrasts")
## [,1] [,2]
## B 116-51 1 0
## B 1-52 0 1
## B 25-50 E -1 -1
## Levels: B 116-51 B 1-52 B 25-50 E
C(da$vari, contr=contr.helmert)
## [1] B 116-51 B 116-51 B 116-51 B 116-51 B 1-52 B 1-52 B 1-52 B 1-52
## [9] B 25-50 E B 25-50 E B 25-50 E B 25-50 E
## attr(,"contrasts")
## [,1] [,2]
## B 116-51 -1 -1
## B 1-52 1 -1
## B 25-50 E 0 2
## Levels: B 116-51 B 1-52 B 25-50 E
##-----------------------------------------------------------------------------
## Significado de cada contraste.
L <- vector(mode="list", length=5)
L[[1]] <- lm(prod~0+vari, data=da)
L[[2]] <- lm(prod~vari, data=da, contr=list(vari=contr.treatment))
L[[3]] <- lm(prod~vari, data=da, contr=list(vari=contr.SAS))
L[[4]] <- lm(prod~vari, data=da, contr=list(vari=contr.sum))
L[[5]] <- lm(prod~vari, data=da, contr=list(vari=contr.helmert))
names(L) <- c("cela", "trt", "SAS", "sum", "helm")
sapply(L, coef)
## cela trt SAS sum helm
## variB 116-51 22.50 22.500 16.500 20.425 20.4250
## variB 1-52 22.27 -0.225 6.000 2.075 -0.1125
## variB 25-50 E 16.50 -6.000 5.775 1.850 -1.9625
L[[6]] <- lm(prod~1, data=da)
coef(L[[6]])
## (Intercept)
## 20.42
coef(L[[6]])-coef(L[[1]])
## variB 116-51 variB 1-52 variB 25-50 E
## -2.075 -1.850 3.925
coef(L[[5]])
## (Intercept) vari1 vari2
## 20.4250 -0.1125 -1.9625
chelm <- coef(L[[1]])
(-chelm[1]+chelm[2])/2 ## 2 = (-1)^2 + 1^2
## variB 116-51
## -0.1125
(-chelm[1]-chelm[2]+2*chelm[3])/6 ## 6 = (-1)^2 + (-1)^2 + 2^2
## variB 116-51
## -1.962
## Cela: estimativas por cela.
## trt: 1 nível é referência, os demais são diferenças.
## sas: k nível é referência, os demais são diferenças.
## sum: intercept é a média geral, os demais são desvios da média.