suppressMessages ( library ( dplyr ) )
suppressMessages ( library ( ggplot2 ) )
options(dplyr.summarise.inform = FALSE)
- Os dados a seguir correspondem a médias de uma variável quantitativa de acordo com o genótipo conjunto de dois locos de marcadores moleculares (do tipo Single Nucleotide Polimorphsms, SNP_A e SNP_B)), para 7 populações sob estudo. As ações gênicas nessas populações são denominadas como: P1: Ação gênica intermediária; P2: Dominância completa; P3: Ação gênica complementar; P4: Epistasia complexa; P5: Dominância parcial; P6: Superdominância; P7: Ação gênica duplicada.
a) Para cada população (P1 a P7) construa o gráfico de perfis de médias correspondente. Com base nesse gráfico, há indicação de efeito de interação entre os dois fatores genéticos (SNP_A e SNP_B)? Comente.
# dados para matrizes
medias1 = t ( matrix ( c ( 20 , 18 , 16 , 17 , 15 , 13 , 14 , 12 , 10) , ncol = 3 ) )
medias2 = t ( matrix ( c ( 20 , 20 , 16 , 20 , 20 , 16 , 14 , 14 , 10) , ncol = 3 ) )
medias3 = t ( matrix ( c ( 20 , 20 , 10 , 20 , 20 , 10 , 10 , 10 , 10) , ncol = 3 ) )
medias4 = t ( matrix ( c ( 38.7 , 6.6 , 3.2 , 4 , 2 , 1.3 , 2.7 , 3.4 , 1.6) , ncol = 3 ) )
medias5 = t ( matrix ( c ( 20 , 19 , 16 , 19 , 18 , 15 , 14 , 13 , 10) , ncol = 3 ) )
medias6 = t ( matrix ( c ( 20 , 21 , 16 , 22 , 23 , 18 , 14 , 15 , 10) , ncol = 3 ) )
medias7 = t ( matrix ( c ( 20 , 20 , 20 , 20 , 20 , 20 , 20 , 20 , 10) , ncol = 3 ) )
# Checagem
medias1
mediasTodas = c ( )
for ( i in 1:7 ){
a = matrix ( get ( paste0 ( "medias" , i ) ) , nrow = 9 , ncol = 1 )
mediasTodas = c ( mediasTodas , a )
}
df = as.data.frame ( cbind ( mediasTodas , rep ( c ( "aa" , "Aa" , "AA" ) , 9*7 ) , c (rep ("bb",3),rep("Bb",3),rep("BB",3) ) , paste0 ("P" ,sort ( rep ( seq ( 1 , 7 ) , 9 ) ) ) ) )
glimpse ( df )
colnames ( df ) = c ( "media" , "genA" , "genB" , "pop" )
df$media = as.numeric ( as.character ( df$media ) )
head ( df , 10 )
df %>% ggplot ( aes ( x = genA , y = media , group = genB , col = genB ) ) +
geom_line ( aes(linetype=genB) ) +
geom_point ( ) +
facet_wrap ( ~pop ) +
labs ( x = "Gene A" , y = "Média" , linetype = "Gene B" , col = "Gene B" , title = "Gráfico de perfil de médias" )
Pode-se observar descritivamente que:
- Não há evidencias de interação entre os genótipos nas populações 1, 2, 5 e 6;
- Há evidencias de interação entre os genótipos nas populações 3, 4 e 7.
b) Vamos agora gerar dados de Delineamentos Completamente Aleatorizados (DCA) com esquema Fatorial 3x3 (Fator SNP_A e fator SNP_B, ambos, em 3 níveis genotípicos). Considere $\sigma^2$ =1 e, com base nas 9 médias apresentadas em P1 a P7, gere $r=25$ valores para cada um dos 9 grupos, sob o modelo de distribuição Normal. Em cada caso, obtenha a tabela de ANOVA e interprete os resultados.
# Gerar os dados para os cruzamentos genotípicos para todas as populações
genA = c ( "aa" , "Aa" , "AA" )
genB = c ( "bb" , "Bb" , "BB" )
n = 25
sigma = 1
set.seed ( 202010 )
# loop nas populações
for ( i in 1:7 ){
# vetor no qual serão armazenados os valores dos cruzamentos genotipicos para cada população
pop = c ( )
# loop nas linhas de cada população
for ( indiceA in 1:length ( genA ) ){
# loop nas colunas de cada população
for ( indiceB in 1:length ( genB ) ){
# medias é a matriz da população i (de 1 a 7)
medias = get ( paste0 ( "medias" , i ) )
# geracao dos valores para a populacao i da combinacao genotipica (indiceA: 1 = aa , 2 = Aa , 3 = AA; indiceB: 1 = bb , 2 = Bb , 3 = BB)
obs = rnorm ( n , mean = medias [ indiceA , indiceB ] , sd = 1 )
# guardando valores no loop
observacoesConjuntas = cbind ( obs , genA [ indiceA ] , genB [ indiceB ] )
pop = rbind ( pop , observacoesConjuntas )
}
}
assign ( paste0 ( "dfCruzamento" , i) , as.data.frame ( pop ) )
assign ( paste0 ( "dfCruzamento" , i) , get ( paste0 ( "dfCruzamento" , i) ) %>% rename ( "Obs" = "obs" , "genA" = "V2" , "genB" = "V3" ) )
assign ( paste0 ( "dfCruzamento" , i) , get ( paste0 ( "dfCruzamento" , i) ) %>% mutate ( Obs = as.numeric ( as.character ( Obs ) ) ) )
assign ( paste0 ( "dfCruzamento" , i) , get ( paste0 ( "dfCruzamento" , i) ) %>% mutate ( genA = as.factor ( genA ) ) )
assign ( paste0 ( "dfCruzamento" , i) , get ( paste0 ( "dfCruzamento" , i) ) %>% mutate ( genB = as.factor ( genB ) ) )
}
# função para imprimir tabelas de anova
printTabela = function ( df ){
print ( anova ( aov ( df$Obs ~ df$genA * df$genB ) ) )
print ( summary ( lm ( df$Obs ~ df$genA * df$genB ) ) )
}
for ( i in 1:7 ){
print (paste ( "População" , i ) )
printTabela ( get ( paste0 ( "dfCruzamento" , i ) ) )
}
Assim como havíamos deduzido a partir dos gráficos, temos que nas populações 3, 4 e 7 há evidencias para se rejeitar a hipótese de que não há interação entre os genótipos, enquanto que nas demais não.
- Considere os dados do arquivo “map.txt” que traz informações de 13 variáveis fenotípicas ("SEX", "Peso", "Altura", "CIRCABD", "GLICOSE", "CTOTAL", "TRIG", "Idade", "mediaPAS", "mediaPAD", "MedPres", "MedCol", "MedDb"), além do genótipo de 102 SNPs ("rs12628452_A" "rs7289830_T", etc.) para uma amostra de 1675 indivíduos.
a) Seguindo o princípio da Aleatorização Mendeliana (Burgess and Thompson, 2015), suponha que o nível genotípico (aa=0, Aa=1, AA=2) de cada SNP foi aleatoriamente atribuído a cada indivíduo. Assim, escolha uma das variáveis quantitativas do banco de dados, e avalie o efeito de cada SNP sobre essa variável por meio de um modelo de ANOVA. Adote um esquema DCA com 1 fator em 3 níveis. Caso seja necessário, proponha uma transformação da variável para atender às premissas clássicas.
base2 = read.table ( "map.txt" , header = T , stringsAsFactors = F )
head ( base2 )
base2$Peso = base2$Peso/1000
head(base2)
fenotipos = c("SEX", "Peso", "Altura", "CIRCABD", "GLICOSE", "CTOTAL", "TRIG", "Idade", "mediaPAS", "mediaPAD", "MedPres", "MedCol", "MedDb")
snps = setdiff ( colnames ( base2 ) , fenotipos )
infosRelevantes = c ( )
for ( coluna in snps ){
baseAux = base2 [ (!is.na ( base2 [ , coluna ] ))&(!is.na(base2$Peso)) , ]
if ( nrow ( baseAux [ baseAux[ , coluna ] == 2 , ] ) == 1 ){
baseAux = baseAux [ baseAux [ , coluna ] < 2 , ]
}
# print ( nrow ( baseAux ) )
baseAux[ , coluna ] = as.factor ( baseAux[ , coluna ] )
# print ( nrow ( baseAux ) )
modelo = lm ( baseAux$Peso ~ baseAux[,coluna] )
# print ( length ( modelo$residuals ) )
obj = anova ( modelo )
transf = "N"
shapiroTest = shapiro.test ( modelo$residuals )
bart = bartlett.test ( modelo$residuals ~ baseAux[ , coluna ] )
ind = snpar::runs.test( modelo$residuals )
if ( shapiroTest$p.value < 0.01 ){
box = EnvStats::boxcox( lm ( baseAux$Peso ~ baseAux[,coluna] ) , objective.name = "Shapiro-Wilk" , optimize=T , lambda = c ( -2,2 ))
modelo = lm ( ((baseAux$Peso^box$lambda)-1)/box$lambda ~ baseAux[,coluna] )
transf = "S"
shapiroTest = shapiro.test ( modelo$residuals )
bart = bartlett.test ( modelo$residuals ~ baseAux[ , coluna ] )
ind = snpar::runs.test( modelo$residuals )
}
infosRelevantes = rbind ( infosRelevantes , c ( coluna , obj$Df[1] , obj$`Pr(>F)`[1] , nrow ( baseAux ) , shapiroTest$p.value , bart$p.value , ind$p.value , transf ) )
}
resultados = infosRelevantes %>% as.data.frame ()
glimpse ( resultados )
colnames ( resultados ) = c ( "snp" , "gl" , "valP" , "n" , "normalidade" , "homocedasticidade" , "indep" , "transf" )
resultados$gl = as.numeric( as.character ( resultados$gl ) )
resultados$valP = as.numeric( as.character ( resultados$valP ) )
resultados$n = as.numeric( as.character ( resultados$n ) )
resultados$normalidade = as.numeric( as.character ( resultados$normalidade ) )
resultados$homocedasticidade = as.numeric( as.character ( resultados$homocedasticidade ) )
resultados$indep = as.numeric( as.character ( resultados$indep ) )
glimpse ( resultados )
sig = 0.01
nrow ( resultados[resultados$indep > sig , ] )
Como pode-se observar, nenhum dos modelos apresenta a premissa da independencia (considerando o nível de significancia de 1%).
print ( paste ( "Existem" , nrow ( resultados [resultados$normalidade > sig & resultados$homocedasticidade > sig , ] ) , "modelos que atendem simultaneamente as premissas de normalidade e homocedasticidade." ) )
print ( "Dos quais apenas os modelos com os snps abaixo possuem efeito sobre a resposta (considerando nível de significancia a 5%) " )
resultados [ resultados$normalidade > sig & resultados$homocedasticidade > sig & resultados$valP < 0.05 , ]
b) Apresente os valores-p da estatística F do teste global da ANOVA em um “gráfico Manhattan” (gráfico de dispersão: -log10(valor-p) x índice do SNP). Interprete os resultados. Discuta sobre possíveis correções para os múltiplos (102) testes realizados simultaneamente.
Primeira das versões¶
# Atribuído cromossomo "aleatório" para os snps
resultados$Chr = 1
# ordenação por nome
resultados = resultados %>% mutate ( tam = nchar ( as.character ( resultados$snp ) ) ) %>% arrange ( tam , snp ) %>% mutate ( Ord = 1:nrow ( resultados ) )
resultados %>% head ( 4 )
# Ordenação do valor-P para calcular a significância por FDR e pelo método de Holm
correcoes = resultados %>% select ( snp , valP ) %>% arrange ( valP )
correcoes %>% head ( 5 )
# Cálculo da significância por FDR e pelo método de Holm
encontraSignif = function ( df , coluna ){
holm = 0
for ( indice in 1 : nrow ( df ) ){
if ( holm == 0 ){
if ( df [ indice , coluna ] < .05/( nrow ( df ) - indice + 1 ) ) {
df [ indice , "Holm" ] = "S"
} else {
df [ indice:nrow ( df ) , "Holm" ] = "N"
holm = 1
}
}
if ( df [ nrow (df) - indice + 1 , coluna ] < .05*(nrow (df) - indice + 1)/nrow ( df ) ){
df [ nrow (df) - indice + 1 , "FDR" ] = "S"
} else {
df [ nrow (df) - indice + 1 , "FDR" ] = "N"
}
}
return ( df )
}
correcoes = encontraSignif ( correcoes , "valP" )
correcoes %>% head ( 5 )
table ( correcoes$FDR , correcoes$Holm )
# Apenas um SNP teve o efeito significante e em ambas as correções
mergeDfs = merge ( resultados , correcoes , by = c ( "snp" , "valP" ) )
head ( mergeDfs )
# Correção de Bonferroni
bF = .05/nrow ( mergeDfs )
# Snps identificados por Holm e FDR
qqman::manhattan ( mergeDfs , chr = "Chr" , bp = "Ord" , p = "valP" , snp = "snp" , genomewideline = -log10(bF) , highlight = mergeDfs[ ( mergeDfs$Holm == "S" ) & ( mergeDfs$FDR == "S" ) , "snp"])
Conclusão¶
Por todos os critérios, o snp rs1041770_A foi o único que, considerando nível de significância a 5%, teve evidências de ter efeito sobre a variável resposta peso. Algumas considerações importantes:
- O critério de Bonferroni é o mais simples e aplicável dos métodos. Pode ser utilizado para testes (ou intervalos) que não necessariamente precisam ser independentes, ou do mesmo tipo, ou relacionados de alguma forma;
- O método de Holm é uma modificação do método de Bonferroni, mais poderoso no sentido de comparar apenas o menor valor-P com $\frac{\alpha}{k}$ (neste caso, deveríamos encontrar a mesma quantidade ou mais snp's significantes do que no método anterior). Porém não pode ser utilizado para se realizar intervalos de confiança;
- O método FDR controla a taxa de descobertas falsas (i.e., a quantidade esperada de rejeições de $H_0$ falsas em relação a todas as rejeições de $H_0$), além de ser menos conservador. Também existe o ponto de que este método é adequado apenas a testes estatisticamente independentes.
c) Para cada SNP, teste seu efeito aditivo, $H_0: \mu_{AA} - \mu_{aa} = 0$. Apresente os resultados desses testes por meio de um “gráfico vulcão” (gráfico de dispersão: -log10(valor-p) x estimativa do contraste ( $\mu_{AA} - \mu_{aa}$ ) ). Interprete os resultados
infosRelevantesContraste = c ( )
for ( coluna in snps ){
baseAux = base2 [ (!is.na ( base2 [ , coluna ] ))&(!is.na(base2$Peso)) , ]
if ( nrow ( baseAux [ baseAux[ , coluna ] == 2 , ] ) <= 1 ){
print ( paste ( "O snp" , coluna , "não será utilizado pois não possui homozigotos AA" ) )
} else{
baseAux[baseAux[,coluna] == 0 , "x1" ] = -1
baseAux[baseAux[,coluna] == 1 , "x1" ] = 0
baseAux[baseAux[,coluna] == 2 , "x1" ] = 1
baseAux[baseAux[,coluna] == 0 , "x2" ] = 0
baseAux[baseAux[,coluna] == 1 , "x2" ] = 1
baseAux[baseAux[,coluna] == 2 , "x2" ] = 0
modelo = lm ( baseAux$Peso ~ baseAux$x1 + baseAux$x2 )
# Aditividade
coefC = modelo$coefficients[2][[1]]
# Dominancia
coefD = modelo$coefficients[3][[1]]
sumMod = summary ( modelo )
# Aditividade
valPC = sumMod$coefficients[2,4]
# Dominancia
valPD = sumMod$coefficients[3,4]
infosRelevantesContraste = rbind ( infosRelevantesContraste , c ( coluna , coefC , valPC , coefD , valPD , nrow ( baseAux ) ) )
}
}
resultadosContraste = infosRelevantesContraste %>% as.data.frame ()
colnames ( resultadosContraste ) = c ( "snp" , "difMediasC" , "valPC" , "difMediasD" , "valPD" , "n" )
resultadosContraste$difMediasC = as.numeric( as.character ( resultadosContraste$difMediasC ) )
resultadosContraste$valPC = as.numeric( as.character ( resultadosContraste$valPC ) )
resultadosContraste$difMediasD = as.numeric( as.character ( resultadosContraste$difMediasD ) )
resultadosContraste$valPD = as.numeric( as.character ( resultadosContraste$valPD ) )
resultadosContraste$n = as.numeric( as.character ( resultadosContraste$n ) )
head ( resultadosContraste )
resultadosContrasteC = resultadosContraste %>% mutate ( logEscalaC = -log10(valPC) ) %>% arrange ( -logEscalaC )
resultadosContrasteC %>% head(10)
CorrecoesContrastesC = encontraSignif ( resultadosContrasteC , "valPC" )
CorrecoesContrastesC %>% head ( 10 )
bfContrastesC = -log10 ( 0.05/ nrow ( CorrecoesContrastesC ) )
CorrecoesContrastesC %>% ggplot ( aes ( x = difMediasC , y = logEscalaC , colour = Holm ) ) +
geom_point () +
geom_hline ( yintercept = bfContrastesC , col = "green" ) +
labs ( x = "Estimativa do efeito aditivo" , y = "-log10(p)" , title = "Gráfico de Vulcão - Aditividade - método de Holm" )
CorrecoesContrastesC %>% ggplot ( aes ( x = difMediasC , y = logEscalaC , colour = FDR ) ) +
geom_point () +
geom_hline ( yintercept = bfContrastesC , col = "green" ) +
labs ( x = "Estimativa do efeito aditivo" , y = "-log10(p)" , title = "Gráfico de Vulcão - Aditividade - método FDR" )
Conclusão¶
Pelos métodos de Bonferroni,Holm e FDR não encontramos evidências de efeito aditivo para nenhum dos snps considerando a variável peso como resposta.
d) Agora, para cada SNP, teste seu efeito de dominância, $H_0: \mu_{Aa} - \frac{(\mu_{AA} - \mu_{aa})}{2} = 0$. Apresente os resultados desses testes por meio de um “gráfico vulcão” (gráfico de dispersão: -log10(valor-p) x estimativa do contraste $H_0: \mu_{Aa} - \frac{(\mu_{AA} - \mu_{aa})}{2} = 0$. Interprete os resultados.
resultadosContrasteD = resultadosContraste %>% mutate ( logEscalaD = -log10(valPD) ) %>% arrange ( -logEscalaD )
resultadosContrasteD %>% head(10)
CorrecoesContrastesD = encontraSignif ( resultadosContrasteD , "valPD" )
CorrecoesContrastesD %>% head ( 10 )
bfContrastesD = -log10 ( 0.05/ nrow ( CorrecoesContrastesD ) )
CorrecoesContrastesD %>% ggplot ( aes ( x = difMediasD , y = logEscalaD , colour = Holm ) ) +
geom_point () +
geom_hline ( yintercept = bfContrastesD , col = "green" ) +
labs ( x = "Estimativa do efeito aditivo" , y = "-log10(p)" , title = "Gráfico de Vulcão - Dominância - método de Holm" )
CorrecoesContrastesD %>% ggplot ( aes ( x = difMediasD , y = logEscalaD , colour = FDR ) ) +
geom_point () +
geom_hline ( yintercept = bfContrastesD , col = "green" ) +
labs ( x = "Estimativa do efeito aditivo" , y = "-log10(p)" , title = "Gráfico de Vulcão - Dominância - método FDR" )
Conclusão¶
Pelos métodos de Bonferroni e Holm, não encontramos evidências de efeito aditivo para nenhum dos snps considerando a variável peso como resposta.
e) Hipoteticamente, assuma que há interesse, e é apropriado, realizar uma meta-análise dos resultados de significância do efeito aditivo desses 102 marcadores. Use as fórmulas apresentadas a seguir para combinar os resultados e obter um valor-p global.
Obs: A fórmula correta é $P = 2\phi(-|Z|)$ - é mais coerente com o que é dito no artigo.
infosNecessarias = CorrecoesContrastesC %>% select ( snp , difMediasC , valPC , n ) %>% mutate ( w = sqrt ( n ) , z = qnorm ( valPC/2 ) , sentido = sign ( difMediasC ) ) %>% mutate ( bigZNum = z*sentido*w )
head( infosNecessarias )
pnorm ( sum(infosNecessarias$bigZNum)/sum(infosNecessarias$w^2) )
2 * pnorm ( - abs ( sum(infosNecessarias$bigZNum)/sum(infosNecessarias$w^2) ) )
print ( paste ( "Z-score:" , sum(infosNecessarias$bigZNum)/sum(infosNecessarias$w^2)))
Baseado no valor calculado, conclui-se com o valor-p global obtido que não há significância de efeito aditivo para a variável peso.
f) Ainda, considerando o banco de dados “map.txt”, escolha um dentre os 102 SNPs e avalie o possível efeito de interação do SNP e do sexo do indivíduo sobre a variável que está sendo analisada. Para tanto, ajuste um modelo ANOVA para um esquema DCA Fatorial 3x2.
Tomando o SNP rs5993848_G para avaliar o possível efeito de interação com a variável sexo, temos que:
snpEscolhido = "rs5993848_G"
dfMod = base2 [ (!is.na ( base2 [ , snpEscolhido ] ))&(!is.na(base2$Peso)) , c ( "SEX" , "Peso" , snpEscolhido ) ]
head ( dfMod )
dfMod$rs5993848_G = as.factor ( dfMod$rs5993848_G )
dfMod$SEX = as.factor ( dfMod$SEX )
anovaMod = aov ( Peso ~ SEX*rs5993848_G , data = dfMod)
anova(anovaMod)
dfMod %>% ggplot ( ) + geom_boxplot ( aes ( y = Peso , x = rs5993848_G , fill = SEX ) ) + facet_wrap(~SEX)
A partir da tabela de ANOVA conclui-se que, apesar de não podermos rejeitar a hipótese de que a interação é nula, temos que tanto o efeito de SNP quanto o efeito de sexo existem sobre a variável resposta peso.
- Os dados a seguir se referem à contagem diferencial de um tipo de pulgão que ocorre em plantas, antes e depois de um tratamento por 10 dias. Três tratamentos, com água (T1), esporos (T2) e com óleo natural (T3) foram usados em 5 plantas para avaliar o controle da infestação pelo pulgão. No experimento, três folhas contaminadas foram selecionadas de cada uma das plantas, e os tratamentos foram então aleatoriamente atribuídos às folhas. A contagem de pulgões em quadrantes definidos nas folhas foi feita antes e depois do tratamento e a diferença registrada. Os dados estão indicados a seguir.
a) Complete a casela faltante (indicada como ?) e realize uma análise de variância destes dados. Que suposições foram adotadas na análise? Essas estão satisfeitas?
# Uma possibilidade
col1 = c ( -7 , 2 , 7 , 11 , 22 , 32 , 10 , 2 , 15 , NA , 2 , 16 , 4 , 5 , 10 )
col2 = rep ( c ( "Agua" , "Esporos" , "Oleo" ) , 5 )
col3 = sort ( rep ( c ( "P1" , "P2" , "P3" , "P4" , "P5" ) , 3 ) )
df3 = as.data.frame ( cbind ( col1 , col2 , col3 ) )
df3 %>% head ( 10 )
df3 = df3 %>% rename ( "Delta" = "col1" , "Trat" = "col2" , "Planta" = "col3" ) %>% mutate ( Delta = as.numeric ( as.character ( Delta ) ) )
head ( df3 )
mod3 = lm ( Delta ~ Trat + Planta , data = df3 )
valorMissing = predict ( object = mod3 , newdata = df3[10,] )
print ( paste ( "Uma possibilidade para o valor faltante é:" , valorMissing[[1]] ) )
Agora que completamos o valor faltante, vamos à ANOVA:
df3 [ 10 , "Delta" ] = valorMissing[[1]]
anova3 = aov ( Delta ~ Trat + Planta , data = df3 )
summary ( anova3 )
# Criando funcao para gerar grafico QQ
graficoQQ = function ( modelo , valoresInf , valoresSup ) {
residuosAjust = modelo$residuals / sigma ( modelo )
numObs = length ( residuosAjust )
quantisRes = quantile ( residuosAjust , c ( ( 1 : ( numObs - 1) )/numObs , ( numObs - 0.5 ) / numObs ) )
quantisNorm = qnorm ( c ( ( 1 : ( numObs - 1) )/numObs , ( numObs - 0.5 ) / numObs ) )
dfQQ = as.data.frame ( cbind ( quantisRes , quantisNorm , valoresInf , valoresSup ) )
print (
dfQQ %>% ggplot ( ) +
geom_point ( aes ( x = quantisNorm , y = quantisRes ) ) +
# geom_line ( aes ( x = quantisNorm , y = valoresInf ) , col = 'red' ) +
# geom_line ( aes ( x = quantisNorm , y = valoresSup ) , col = 'red' ) +
geom_abline ( intercept = 0 , slope = 1 ) +
ggtitle ( "Gráfico dos quantis de resíduos ajustados em relação a quantis normais" ) +
labs ( x = "Quantis normais" , y = "Quantis dos resíduos" )
)
}
# Homocedasticidade
graficoHomo = function ( modelo ) {
dfHomocedast = as.data.frame ( cbind ( modelo$residuals / sigma ( modelo ) , modelo$fitted.values ) )
names ( dfHomocedast ) = c ( "Residuo" , "Ajuste" )
print (
dfHomocedast %>% ggplot ( ) +
geom_point ( aes ( x = Ajuste , y = Residuo ) ) +
geom_hline ( yintercept = 0 , col = "red" ) +
ggtitle ( "Gráfico de resíduos ajustados em relação a valores ajustados" ) +
labs ( x = "Valores ajustados" , y = "Resíduos ajustados" )
)
}
# Residuos
graficoResiduos = function ( modelo ) {
residuosDf = as.data.frame ( modelo$residuals )
colnames ( residuosDf ) = c ( "residuo" )
residuosDf$indice = row ( residuosDf )
print (
residuosDf %>% ggplot ( ) +
geom_point ( aes ( x = indice , y = residuo ) ) +
geom_hline ( yintercept = 0 , col = "red" ) +
ggtitle ( "Gráfico dos resíduos" ) +
labs ( x = "Índice" , y = "Resíduo" )
)
}
# plota tudo
plotaGeral = function ( modelo , valoresInf , valoresSup ){
graficoQQ ( modelo , valoresInf , valoresSup )
graficoHomo ( modelo )
graficoResiduos ( modelo )
}
plotaGeral ( lm ( Delta ~ Trat + Planta , data = df3 ) , NULL , NULL )
print ( bartlett.test ( anova3$residuals ~ df3[ , "Trat" ] ) )
print ( snpar::runs.test( anova3$residuals ) )
print ( shapiro.test ( anova3$residuals ) )
Levando-se em consideração que existem poucas observações no experimento, não se percebe grandes desvios das premissas.
b) Compare as fontes de variação presentes no experimento.
bartlett.test ( anova3$residuals ~ df3[ , "Planta" ] )
Não há evidências para se rejeitar a hipótese de que exista diferença entre as fontes de variação.
c) Suspeitando que a hipótese de normalidade não é apropriada a estes dados, construa um Teste de Aleatorização para o efeito global dos tratamentos.
Uma possibilidade de teste seria criar um algoritmo para, fixados os valores observados, permutar-se os grupos aos quais os mesmos foram observados e a cada iteração, realizar o cálculo da estatística F (desta forma, teremos uma distribuição de comparações entre grupos). Assim, se o valor F que calcularmos for extremo nesta distribuição, teremos indícios de que o "grupo" aos quais as observações pertencem têm algum tipo de associação com a resposta.
nSim = 10^3
valF = c ( )
set.seed(123)
for ( i in 1:nSim ){
trat = sample ( df3$Trat , nrow ( df3 ) , replace = F )
resp = df3$Delta
tabela = anova ( aov ( resp ~ trat ) )
valF = c ( valF , tabela$`F value`[1] )
}
dfValf = valF %>% as.data.frame ()
colnames ( dfValf ) = c ( "fAleat" )
dfValf %>% head ( 5 )
tabelaOrig = anova ( aov ( resp ~ trat ) )
valFOrig = tabelaOrig$`F value`[1]
print ( paste ( "A estatística F dos dados originais é:" , round ( valFOrig , 3 ) ) )
dfValf %>% ggplot ( aes ( x = fAleat ) ) +
geom_histogram ( bins = 50 ) +
geom_vline ( xintercept = valFOrig ) +
geom_vline ( xintercept = 1 , col = "red" ) +
labs ( x = "Valores observados" , y = "Frequência absoluta" , title = "Histograma da distribuição de valores F do teste de aleatorização" )
valP = 100 * ( sum ( dfValf$fAleat < valFOrig ) + sum ( dfValf$fAleat > 1/valFOrig ) ) / nrow ( dfValf )
print ( paste0 ( "Temos que o valor-p observado para o teste bicaudal é dado por: " , valP , "%" ) )
Considerando o teste de aleatorização realizado, temos que não há evidências para se rejeitar a hipótese de que os tratamentos são iguais (valor-p dado por 58,9%).