##============================================================================= ## AF-722: Modelagem e análise de dados experimentais com o programa ## computacional R ## www.leg.ufpr.br/af722 ## ## Programa de Pós Graduação em Agronomia - 2014/1 ## Prof Dr. Walmes Zeviani ## LEG - DEST - UFPR ## ## Script com algumas matrizes e operações matriciais relacioandas às ## restrições usualmente adotadas para estimação em modelos lineares. ## ##============================================================================= ##----------------------------------------------------------------------------- ## Usando o R para casualizar os níveis às unidades experimentais. ## Fator de 4 níveis com 4 repetições ## x <- factor(rep(c("a","b","c","d"), each=3)) x <- gl(4, 3, labels=c("a","b","c","d")) x levels(x) nlevels(x) table(x) ## Casualiação dos níveis para um DIC. sample(x) sample(x) sample(x) ##----------------------------------------------------------------------------- ## Tipo de parametrizações para níveis categóricos. ## Seja o vetor com as médias populacionais dos níveis. mu <- c(a=5, b=5, c=10, d=4) mean(mu) ## Número de níveis do fator x. k <- nlevels(x); k ## Modelo: \mu_i = \mu+\alpha_i, i=1,...,k, sujeito à uma restrição ## parametrica. ## Contraste tipo soma zero. ## A soma dos \alpha_i é zero. Assim, \mu representa a média geral. L <- cbind(1, contr.sum(k)) L ## [,1] [,2] [,3] [,4] ## 1 1 0 0 ## \mu+\alpha_1 ## 1 0 1 0 ## \mu+\alpha_2 ## 1 0 0 1 ## \mu+\alpha_3 ## 1 -1 -1 -1 ## \mu-(\alpha_1+\alpha_2+\alpha_3) solve(L) solve(L, mu) ## \mu, \alpha_1, \alpha_2, \alpha_3. ## Contraste que zera um nível 1 (primeiro), \alpha_1 = 0 (padrão do R). L <- cbind(1, contr.treatment(k)) ## 2 3 4 ## 1 0 0 0 ## \mu = \mu_1, (\alpha_1 = 0) ## 1 1 0 0 ## \mu+\alpha_2 (diferença de \mu_2 para \mu_1) ## 1 0 1 0 ## \mu+\alpha_3 (diferença de \mu_3 para \mu_1) ## 1 0 0 1 ## \mu+\alpha_4 (diferença de \mu_4 para \mu_1) solve(L) ## 1 2 3 4 ## 1 0 0 0 ## \mu = \mu_1 ## -1 1 0 0 ## contraste ## -1 0 1 0 ## contraste ## -1 0 0 1 ## contraste solve(L, mu) ## Contraste que zera um nível k (último), \alpha_k = 0 (padrão do SAS). L <- cbind(1, contr.SAS(k)) L solve(L) MASS::fractions(solve(L)) ## 1 2 3 4 ## 1/4 1/4 1/4 1/4 ## média geral ## -1/2 1/2 0 0 ## 1 vs 2 ## -1/6 -1/6 1/3 0 ## 1+2 vs 3 ## -1/12 -1/12 -1/12 1/4 ## 1+2+3 vs 4 solve(L, mu) ## Contraste de Helmert, alpha_i vs todos alpha_j onde j < i, ou seja, ## 1 vs 2, 1+2 vs 3, 1+2+3 vs 4, etc. L <- cbind(1, contr.helmert(k)) L solve(L) solve(L, mu) ##----------------------------------------------------------------------------- ## A model matriz faz a matriz do delineamento segundo a restrição escolhida X <- model.matrix(~x); X X <- model.matrix(~x, contrasts=list(x=contr.SAS)); X X <- model.matrix(~x, contrasts=list(x=contr.sum)); X X <- model.matrix(~x, contrasts=list(x=contr.helmert)); X ##----------------------------------------------------------------------------- ## Como montar a matriz de delineamento? ## Precisa da matriz de restrições L e da matriz de incidência (0 se a ## observação não pertence ao nível i, 1 se pertence). ## Matriz de pertence ou não. l <- levels(x); l w <- sapply(l, "==", x) w w*1 I <- w*1 L <- cbind(1, contr.helmert(k)) I%*%L model.matrix(~x, contrasts=list(x=contr.helmert)) L <- cbind(1, contr.sum(k)) I%*%L model.matrix(~x, contrasts=list(x=contr.sum)) ##----------------------------------------------------------------------------- ## Como a partir do vetor de efeitos obter o vetor de média? ## Médias. mu ## Vetor de efeitos sob a restrição. ef <- solve(L, mu); ef ## Médias a partir de L e o vetor de efeitos. L%*%ef ##-----------------------------------------------------------------------------