#=============================================================================== #. Name : gexp_quant_fe #. Author : J.C.Faria # Date : 06/10/2022 01:33:35 # Aim : Gerar/simular e analisar experimentos quantitativos fatoriais nos # delineamentos experimentais básicos: # - Del. Inteiramente ao Acaso - DIC # - Del. Blocos Casualizados - DBC # - Del. Quadrado Latino - DQL # Mail : <<< joseclaudio.faria@gmail.com >>> #=============================================================================== #.. Instalação do pacote 'gexp' #... do GitHub: via devtools #install.packages('devtools', dep=TRUE) #library(devtools) #install_github('ivanalaman/gexp') #... do GitHub: via pacman #install.packages('pacman', dep=TRUE) #p_load_gh('ivanalaman/gexp') #... do CRAN # install.packages('gexp') #.. Anexando o pacote #... Necessita gexp >= 1.0.2 library(gexp) packageVersion('gexp') >= '1.0.2' #.. ___FE/DIC___ set.seed(1) fe_dic <- gexp(mu=10, r=3, fe=list(f1=c(2, 3), f2=c(5, # b1* -10, # b2* 0, # b3 0)), # b4 fl=list(TIP=paste0('t', 1:2), X=seq(0, 40, by=10)), inte=c(1, 15, rep(1, 6)), type='FE') str(fe_dic) summary(fe_dic) par(mfrow=c(1, 1)) plot(fe_dic, main='FE sob DIC') # Gráfico dos efeitos par(mfrow=c(2, 1)) with(fe_dic$dfm, interaction.plot(TIP, X, Y1, xlab='TIP', ylab='Y1', col='blue3', lwd=2, ylim=c(min(Y1), max(Y1)))) with(fe_dic$dfm, interaction.plot(X, TIP, Y1, xlab='X', ylab='Y1', col='blue3', lwd=2, ylim=c(min(Y1), max(Y1)))) # ANOVA av_fe_dic <- aov(Y1 ~ TIP + factor(X) + TIP*factor(X), fe_dic$dfm) summary(av_fe_dic) ## Calculando as médias de tratamentos (por TIP) (fe_dic_m <- aggregate(Y1 ~ TIP + X, data=fe_dic$dfm, mean)) ## Subsets por 'TIP' (fe_dic_m_t1 <- subset(fe_dic_m, sub= TIP == 't1')) (fe_dic_m_t2 <- subset(fe_dic_m, sub= TIP == 't2')) ## Diagrama de dispersão com as médias do níveis de 'X' par(mfrow=c(1, 1)) with(fe_dic_m, plot(Y1 ~ X, ylab='m(Y1)', ylim=c(min(Y1), 110/100*max(Y1)), pch=c(8, 19), col=c('red', 'blue'))) ## Ajustando modelo linear grau II e visualizando os resultados reg_t1 <- lm(Y1 ~ X + I(X^2), data=fe_dic_m_t1) summary(reg_t1) reg_t2 <- lm(Y1 ~ X + I(X^2), data=fe_dic_m_t2) summary(reg_t2) ## Visualizar médias e modelo linear ajustado lines(spline(fe_dic_m_t1$X, fitted(reg_t1), n=1e3), # n -> argumento de suavização lwd=2, col='red') lines(spline(fe_dic_m_t2$X, fitted(reg_t2), n=1e3), # n -> argumento de suavização lwd=2, col='blue') legend('topleft', title='TIP', pch=c(8, 19), lty = c(2, 1), legend=c('t1', 't2'), col=c('red', 'blue')) #.. ___FE/DBC___ set.seed(2) fe_dbc <- gexp(mu=10, r=1, fe=list(f1=c(2, 3), f2=c(5, # b1* -10, # b2* 0, # b3 0)), # b4 fl=list(TIP=paste0('t', 1:2), X=seq(0, 40, by=10)), inte=c(1, 15, rep(1, 6)), blke=c(1, 2, 3), blkl=list(BLO=paste0('b', 1:3)), type='FE', des='RCBD') str(fe_dbc) summary(fe_dbc) par(mfrow=c(1, 1)) plot(fe_dbc, main='FE sob DBC') # Gráfico dos efeitos par(mfrow=c(2, 1)) with(fe_dbc$dfm, interaction.plot(TIP, X, Y1, xlab='TIP', ylab='Y1', col='blue3', lwd=2, ylim=c(min(Y1), max(Y1)))) with(fe_dbc$dfm, interaction.plot(X, TIP, Y1, xlab='X', ylab='Y1', col='blue3', lwd=2, ylim=c(min(Y1), max(Y1)))) # ANOVA av_fe_dbc <- aov(Y1 ~ BLO + TIP + factor(X) + TIP*factor(X), fe_dbc$dfm) summary(av_fe_dbc) ## Calculando as médias de tratamentos (por TIP) (fe_dbc_m <- aggregate(Y1 ~ TIP + X, data=fe_dbc$dfm, mean)) ## Subsets por 'TIP' (fe_dbc_m_t1 <- subset(fe_dbc_m, sub= TIP == 't1')) (fe_dbc_m_t2 <- subset(fe_dbc_m, sub= TIP == 't2')) ## Diagrama de dispersão com as médias do níveis de 'X' par(mfrow=c(1, 1)) with(fe_dbc_m, plot(Y1 ~ X, ylab='m(Y1)', ylim=c(min(Y1), 110/100*max(Y1)), pch=c(8, 19), col=c('red', 'blue'))) ## Ajustando modelo linear grau II e visualizando os resultados reg_t1 <- lm(Y1 ~ X + I(X^2), data=fe_dbc_m_t1) summary(reg_t1) reg_t2 <- lm(Y1 ~ X + I(X^2), data=fe_dbc_m_t2) summary(reg_t2) ## Visualizar médias e modelo linear ajustado lines(spline(fe_dbc_m_t1$X, fitted(reg_t1), n=1e3), # n -> argumento de suavização lwd=2, col='red') lines(spline(fe_dbc_m_t2$X, fitted(reg_t2), n=1e3), # n -> argumento de suavização lwd=2, col='blue') legend('topleft', title='TIP', pch=c(8, 19), lty = c(2, 1), legend=c('t1', 't2'), col=c('red', 'blue')) #.. ___FE/DQL___ set.seed(3) fe_dql <- gexp(mu=10, r=1, fe=list(f1=c(2, 3), f2=c(5, # b1* -10, # b2* 0, # b3 0)), # b4 fl=list(TIP=paste0('t', 1:2), X=seq(0, 40, by=10)), inte=c(1, 15, rep(1, 6)), rowe=rep(1, 10), rowl=list(LIN=paste0('l', 1:10)), cole=rep(1, 10), coll=list(COL=paste0('c', 1:10)), type='FE', des='LSD') str(fe_dql) summary(fe_dql) par(mfrow=c(1, 1)) plot(fe_dql, main='FE sob DQL') # Gráfico dos efeitos par(mfrow=c(2, 1)) with(fe_dql$dfm, interaction.plot(TIP, X, Y1, xlab='TIP', ylab='Y1', col='blue3', lwd=2, ylim=c(min(Y1), max(Y1)))) with(fe_dql$dfm, interaction.plot(X, TIP, Y1, xlab='X', ylab='Y1', col='blue3', lwd=2, ylim=c(min(Y1), max(Y1)))) fe_dql$dfm$X <- as.numeric(fe_dql$dfm$X) str(fe_dql) # ANOVA av_fe_dql <- aov(Y1 ~ LIN + COL + TIP + factor(X) + TIP*factor(X), fe_dql$dfm) summary(av_fe_dql) ## Calculando as médias de tratamentos (por TIP) (fe_dql_m <- aggregate(Y1 ~ TIP + X, data=fe_dql$dfm, mean)) ## Subsets por 'TIP' (fe_dql_m_t1 <- subset(fe_dql_m, sub= TIP == 't1')) (fe_dql_m_t2 <- subset(fe_dql_m, sub= TIP == 't2')) ## Diagrama de dispersão com as médias do níveis de 'X' par(mfrow=c(1, 1)) with(fe_dql_m, plot(Y1 ~ X, ylab='m(Y1)', ylim=c(min(Y1), 110/100*max(Y1)), pch=c(8, 19), col=c('red', 'blue'))) ## Ajustando modelo linear grau II e visualizando os resultados reg_t1 <- lm(Y1 ~ X + I(X^2), data=fe_dql_m_t1) summary(reg_t1) reg_t2 <- lm(Y1 ~ X + I(X^2), data=fe_dql_m_t2) summary(reg_t2) ## Visualizar médias e modelo linear ajustado lines(spline(fe_dql_m_t1$X, fitted(reg_t1), n=1e3), # n -> argumento de suavização lwd=2, col='red') lines(spline(fe_dql_m_t2$X, fitted(reg_t2), n=1e3), # n -> argumento de suavização lwd=2, col='blue') legend('topleft', title='TIP', pch=c(8, 19), lty = c(2, 1), legend=c('t1', 't2'), col=c('red', 'blue'))