CE 063 - Planejamento e análise de experimentos irregulares

Universidade Federal do Paraná
Prof. Dr. Walmes M. Zeviani
Curso de Estatística - 2015/1
Lab. de Estatística e Geoinformação - LEG
Departamento de Estatística - UFPR

Revisão de modelos lineares

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.