##====================================================================== ## Controle de Processos Industriais www.leg.ufpr.br/ce074 ## Curso de Estatística - 2015/2 ## Prof. Fernando Mayer LEG - DEST - UFPR ##====================================================================== ##====================================================================== ## Experimentos fatoriais 2^2 simples ## Montgomery, Cap. 14, p. 366- ##====================================================================== ##---------------------------------------------------------------------- ## Fatorial 1 fat1 <- data.frame(A = c("A1", "A1", "A2", "A2"), B = c("B1", "B2", "B1", "B2"), Resp = c(10, 20, 30, 40)) fat1 levels(fat1$A) ## Fatorial 2 fat2 <- data.frame(A = c("A1", "A1", "A2", "A2"), B = c("B1", "B2", "B1", "B2"), Resp = c(10, 20, 30, 0)) fat2 ## Visualizacao par(mfrow = c(1, 2)) with(fat1, interaction.plot(A, B, Resp)) with(fat2, interaction.plot(A, B, Resp)) par(mfrow = c(1, 1)) wireframe(Resp ~ A * B, data = fat1) wireframe(Resp ~ A * B, data = fat2) ##---------------------------------------------------------------------- ## Estimativas dos efeitos ## Contraste soma zero m3 <- lm(Resp ~ A * B, data = fat1, contrasts = list(A = "contr.sum", B = "contr.sum")) print(coef(m3), digits = 0) model.matrix(m3) ## Media geral mean(fat1$Resp) ## A mean(fat1$Resp) + coef(m3)[2] # media de A1 mean(fat1$Resp[fat1$A == "A1"]) mean(fat1$Resp) - coef(m3)[2] # media de A2 mean(fat1$Resp[fat1$A == "A2"]) ## B mean(fat1$Resp) + coef(m3)[3] # media de B1 mean(fat1$Resp[fat1$B == "B2"]) mean(fat1$Resp) - coef(m3)[3] # media de B2 mean(fat1$Resp[fat1$B == "B2"]) ## A:B mean(fat1$Resp) + coef(m3)[4] # media de A1B1:A2B2 mean(c(fat1$Resp[fat1$A == "A1" & fat1$B == "B1"], fat1$Resp[fat1$A == "A2" & fat1$B == "B2"])) mean(fat1$Resp) - coef(m3)[4] # media de A2B1:A1B2 mean(c(fat1$Resp[fat1$A == "A2" & fat1$B == "B1"], fat1$Resp[fat1$A == "A1" & fat1$B == "B2"])) ## Na mao, E invertendo a ordem do soma zero para # A1 = -A2 (ao inves de A2 = -A1) # B1 = -B2 (ao inves de B2 = -B1) y <- fat1$Resp X <- rbind(c(1, -1, -1, 1), c(1, -1, 1, -1), c(1, 1, -1, -1), c(1, 1, 1, 1)) (Xt <- t(X)) Xt %*% X solve(Xt %*% X) MASS::fractions(solve(Xt %*% X)) Xt %*% y ## Soma total sum(fat1$Resp) ## Soma catalisador sum(fat1$Resp[fat1$A == "A2"]) - sum(fat1$Resp[fat1$A == "A1"]) ## Soma temperatura sum(fat1$Resp[fat1$B == "B2"]) - sum(fat1$Resp[fat1$B == "B1"]) ## Soma interacao sum(fat1$Resp[fat1$A == "A1" & fat1$B == "B1"], fat1$Resp[fat1$A == "A2" & fat1$B == "B2"]) - sum(fat1$Resp[fat1$A == "A1" & fat1$B == "B2"], fat1$Resp[fat1$A == "A2" & fat1$B == "B1"]) solve(Xt %*% X) %*% Xt %*% y (solve(Xt %*% X) %*% Xt %*% y) * 2 ## Media geral mean(fat1$Resp) ## Media catalisador (mean(fat1$Resp[fat1$A == "A1"]) - mean(fat1$Resp[fat1$A == "A2"]))/2 ## Media temperatura (mean(fat1$Resp[fat1$B == "B1"]) - mean(fat1$Resp[fat1$B == "B2"]))/2 ## Media interacao ((mean(c(fat1$Resp[fat1$A == "A1" & fat1$B == "B1"], fat1$Resp[fat1$A == "A2" & fat1$B == "B2"]))) - (mean(c(fat1$Resp[fat1$A == "A1" & fat1$B == "B2"], fat1$Resp[fat1$A == "A2" & fat1$B == "B1"]))))/2 ##====================================================================== ## Dados dos experimentos com as bolhas (em sala) ##====================================================================== url <- "http://www.leg.ufpr.br/~fernandomayer/ensino/ce074-2015-02/" url.esq <- "dados/experimento_bolhas/dados_grupo_esquerda.csv" esq <- read.csv(paste0(url, url.esq)) url.dir <- "dados/experimento_bolhas/dados_grupo_direita.csv" dir <- read.csv(paste0(url, url.dir)) str(esq) str(dir) ##====================================================================== ## Dados agregados ## Esquerda esq.agg <- aggregate(TEMPO ~ CONCENTRACAO + TIPO, data = esq, FUN = function(x) round(mean(x))) esq.agg par(mfrow = c(1, 2)) with(esq.agg, interaction.plot(TIPO, CONCENTRACAO, TEMPO, ylim = c(40, 70), main = "Esquerda")) with(esq.agg, interaction.plot(CONCENTRACAO, TIPO, TEMPO, ylim = c(40, 70), main = "Esquerda")) par(mfrow = c(1, 1)) ## Direita dir.agg <- aggregate(TEMPO ~ CONCENTRACAO + TIPO, data = dir, FUN = function(x) round(mean(x))) dir.agg par(mfrow = c(1, 2)) with(dir.agg, interaction.plot(TIPO, CONCENTRACAO, TEMPO, ylim = c(40, 70), main = "Direita")) with(dir.agg, interaction.plot(CONCENTRACAO, TIPO, TEMPO, ylim = c(40, 70), main = "Direita")) par(mfrow = c(1, 1)) ##---------------------------------------------------------------------- ## Modelos ## Esquerda me <- lm(TEMPO ~ TIPO * CONCENTRACAO, data = esq.agg, contrasts = list(TIPO = "contr.sum", CONCENTRACAO = "contr.sum")) print(coef(me)) model.matrix(me) ## Efeitos calculados ## TIPO (60+60-62-46)/2 ## Passando do TIPO A para o TIPO B, existe um aumento de 6s ## CONCENTRACAO (62+60-46-60)/2 ## Passando de CONC baixa para CONC ALTA, existe um aumento de 8s ## Interacao (divergente) (60+46-60-62)/2 ## Porque os resultados nao sao os mesmos? ##====================================================================== ## Efeitos e coeficientes ## Efeito ## ef = contraste/2 ## O coeficiente de regressão é metade da estimativa do efeito, visto ## que os coeficientes de regeressão medem o efeito de uma variação ## unitária em X sobre a média de Y, e a estimativa do efeito está ## baseada na variação de duas unidades (entre -1 e +1). ## Portanto, ## Coeficiente do efeito ## beta = efeito/2 ##====================================================================== ## Direita md <- lm(TEMPO ~ TIPO * CONCENTRACAO, data = dir.agg, contrasts = list(TIPO = "contr.sum", CONCENTRACAO = "contr.sum")) print(coef(md)) model.matrix(md) ## Efeitos calculados ## TIPO (66+51-63-61)/2 ## CONCENTRACAO (63+51-66-61)/2 ## Passando de CONC baixa para CONC ALTA, existe um aumento de 8s ## Interacao (divergente) (51+61-63-66)/2