Ajuda com análise conjunta de dois experimentos com parcelas subdividas
Caros colegas, Peço ajuda a vocês para conseguir entender e analisar o seguinte experimento: - são dois experimentos feitos lado a lado: um com plantio direto, outro com convencional. - em cada local, existe um experimento de parcelas subdividades, sendo o fator manejo na parcela e dose de nitrogenio na subparcela. Diante disso, gostaria de saber como expressar o modelo da ANOVA no R, a fim de fazer comparação de médias, e mais pra frente, análise multivariada via discriminante canonica. Sei que são necessárias duas coisas: - Dividir o erro dentro da abordagem da análise conjunta. - Dividir o erro na parcela subdividida. Tenho utilizado a seguinte expressão, mas não estou confiante quanto a expressão. Infelizmente nunca trabalhei com experimento mais complexos como este. aov_N <- aov(N_kg_ha ~ Manejo * Rotacao * Nitrogenio + Bloco + Error(Nitrogenio/Bloco/Manejo), data = dados) summary(aov_N) Tanto o croqui quanto os dados estão neste link: https://drive.google.com/drive/folders/1rjVv9KwtE5A7y-VPXTV9Y2GlilKLXK2F?usp... Agradeço se alguém puder me ajudar com alguma orientação. José Lucas.
Lucas,
Que bom que forneceu os dados e croqui.
O seu croqui indica que você está usando blocos incompletos balanceados,
pelo menos se considerar que o bloco é o retângulo 4 x 4 parcelas. É isso
mesmo? Interessante.
No entanto, a tabela de dados que você mandou parece estar considerando
como bloco a linha que contém 24 parcelas.
No meu modo de ver, o bloco com essas dimensões (tão comprido), não deve
apresentar a homogeneidade que se espera.
É bem provável que a parcela 1C tenha mais semelhança em condições com as
parcelas 4B e 8D, por exemplo, do que com as parcelas 33C e 41D.
Talvez você possa analisar das duas formas.
De qualquer forma, o modelo para o seu experimento fica conforme código a
seguir.
> library(lattice)
>
> url <- "Dados_Conj_2ExpSub.csv"
> tb <- read.csv2(url)
> str(tb)
'data.frame': 96 obs. of 6 variables:
$ FID : Factor w/ 96 levels "12A","12C","13A",..: 9 43 53 87 95 1 3
5 7 11 ...
$ Manejo : Factor w/ 2 levels "PC","PD": 2 2 2 2 2 2 2 2 2 2 ...
$ Rotacao : Factor w/ 6 levels "A","CJ","CO",..: 2 3 5 4 3 6 2 5 1 5 ...
$ Nitrogenio: Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
$ Bloco : int 1 2 3 4 1 2 3 4 1 2 ...
$ N_kg_ha : int 209 234 209 158 101 96 46 118 130 162 ...
>
> tb <- transform(tb,
+ Bloco = factor(Bloco))
>
> xtabs(~Manejo + Rotacao, data = tb)
Rotacao
Manejo A CJ CO CS P S
PC 8 8 8 8 8 8
PD 8 8 8 8 8 8
> ftable(xtabs(~Manejo + Bloco + Rotacao, data = tb))
Rotacao A CJ CO CS P S
Manejo Bloco
PC 1 2 2 2 2 2 2
2 2 2 2 2 2 2
3 2 2 2 2 2 2
4 2 2 2 2 2 2
PD 1 2 2 2 2 2 2
2 2 2 2 2 2 2
3 2 2 2 2 2 2
4 2 2 2 2 2 2
>
> xyplot(N_kg_ha ~ Rotacao | Manejo,
+ groups = Nitrogenio,
+ type = c("p", "a"),
+ auto.key = TRUE,
+ data = tb)
>
> # Modelo para DBC em diferentes locais:
> # ~ Local + Error(Local/Bloco) + Trat + Local:Trat.
> # Modelo de parcela subdividida em DBC:
> # ~ Bloco + Trat + Error(Bloco:Trat) + Subtrat + Trat:Subtrat
> #
> # Juntanto os dois termos de erro:
> # ~ Error(Local/Bloco) + Error(Bloco:Trat) = Error(Local/Bloco:Trat)
> #
> # Juntanto os termos sistemáticos.
> # ~ Local + Trat + Local:Trat + Subtrat + Trat:Subtrat +
> # Error(Local/Bloco:Trat)
>
> m0 <- aov(N_kg_ha ~
+ Manejo +
+ Manejo/Bloco +
+ Rotacao +
+ Manejo:Rotacao +
+ Nitrogenio +
+ Rotacao:Nitrogenio +
+ Error(Manejo/Bloco:Rotacao),
+ data = tb)
Warning message:
In aov(N_kg_ha ~ Manejo + Manejo/Bloco + Rotacao + Manejo:Rotacao + :
Error() model is singular
>
> summary(m0)
Error: Manejo
Df Sum Sq Mean Sq
Manejo 1 504.2 504.2
Error: Manejo:Bloco:Rotacao
Df Sum Sq Mean Sq F value Pr(>F)
Rotacao 5 16367 3273 1.521 0.213
Manejo:Bloco 6 16966 2828 1.314 0.281
Manejo:Rotacao 5 6213 1243 0.577 0.717
Residuals 30 64567 2152
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Nitrogenio 1 13728 13728 6.561 0.0141 *
Rotacao:Nitrogenio 5 5636 1127 0.539 0.7457
Residuals 42 87874 2092
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> # NOTE: não haverá uma estatística F apropriada para testar o efeito de
> # Manejo já que não existe repetição genuína para os níveis desse fator,
> # assim não existe um erro de extrato para ser o denominador da
> # estatística F.
>
> # Trocando Manejo/Bloco por Manejo:Bloco no termo de Error().
> m0 <- aov(N_kg_ha ~
+ Manejo +
+ Manejo/Bloco +
+ Rotacao +
+ Manejo:Rotacao +
+ Nitrogenio +
+ Rotacao:Nitrogenio +
+ Error(Manejo:Bloco:Rotacao),
+ data = tb)
Warning message:
In aov(N_kg_ha ~ Manejo + Manejo/Bloco + Rotacao + Manejo:Rotacao + :
Error() model is singular
>
> summary(m0)
Error: Manejo:Bloco:Rotacao
Df Sum Sq Mean Sq F value Pr(>F)
Manejo 1 504 504 0.234 0.632
Rotacao 5 16367 3273 1.521 0.213
Manejo:Bloco 6 16966 2828 1.314 0.281
Manejo:Rotacao 5 6213 1243 0.577 0.717
Residuals 30 64567 2152
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Nitrogenio 1 13728 13728 6.561 0.0141 *
Rotacao:Nitrogenio 5 5636 1127 0.539 0.7457
Residuals 42 87874 2092
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> # NOTE: com uma mudança sutil na fórmula, muda-se para uma forma não
> # condizente com a realidade, o delineamento que passa a indicar que as
> # combinações Manejo x Rotação foram casualizados às parcelas dos blocos
> # (fatorial completo na parcela) e dessa forma tem-se uma estatística F
> # para Manejo, mas que não é condizente com a forma como o experimento
> # foi de fato feito. Então está errada essa especificação.
>
> # Declarando o primeiro modelo com a `lme4`.
> library(lme4)
>
> mm0 <- lmer(N_kg_ha ~
+ Manejo +
+ Manejo/Bloco +
+ Rotacao +
+ Manejo:Rotacao +
+ Nitrogenio +
+ Rotacao:Nitrogenio +
+ (1 | Manejo/Bloco:Rotacao),
+ data = tb)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
>
> anova(mm0)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Manejo 1 391.0 391.0 0.1869
Rotacao 5 15911.7 3182.3 1.5210
Nitrogenio 1 13728.2 13728.2 6.5613
Manejo:Bloco 6 16494.0 2749.0 1.3139
Manejo:Rotacao 5 6040.2 1208.0 0.5774
Rotacao:Nitrogenio 5 5636.1 1127.2 0.5387
> VarCorr(mm0)
Groups Name Std.Dev.
Bloco:Rotacao:Manejo (Intercept) 5.4722
Manejo (Intercept) 3.3708
Residual 45.7415
À disposição.
Walmes.
Boa noite prezados,
Mais ou menos no contexto da dúvida do José, gostaria de ter uma luz com relação a análise dos dados experimentais que estou trabalhando.
Experimento em dois locais, ambos em delineamento em blocos casualizados em esquema de parcela subdivididas, com quatro repetições. O manejo foi basicamente o mesmo para os dois locais.- parcela: epoca de inicio da cultura (6)- subparcela: cultivar (6)- foram realizadas algumas colheitas antes de colheita final do experimento, com cerca de 15 meses (trabalho com cana também).
Estou utilizando os códigos em anexo para realizar uma anova, com a ajuda dos pacotes 'agricolae' e 'ExpDes'. Tambem utilizei o modelo que o Walmes sugeriu. Pelo que puder checar, os três estão chegando na mesma coisa, o que é bom e indica que talvez eu esteja no caminho certo.Estou rodando para cada local e época de colheita separados (veja os subsets).
As dúvidas, além de saber se estou ou não num caminho lúcido em termos da estatística:1) Em alguns casos houve violação de alguma pressuposição (estou testando norm. dos residuos e homogeneidade apenas), como para a variável resposta 3, no local II colhido na idade 8.- Procedo com alguma transformação ? tipo utilizar boxcox e ver qual melhor forma de transformar os dados ?- Ou Parto para outro tipo de abordagem ? alguma sugestão ?2) Há alguma forma de análise considerando as todas as idades ?3) Ao menos duas idades e 3 dos cultivares são comuns em ambos os locais, seria muito abusado tentar realizar uma análise conjunta apenas para essas cultivares? (com abordagem similar ao proposto para o caso do José)
Agradeço quem puder me dar um help!Até maisHenrique
Em sexta-feira, 20 de setembro de 2019 16:40:56 GMT-4, Walmes Zeviani por (R-br) <r-br@listas.c3sl.ufpr.br> escreveu:
Lucas,
Que bom que forneceu os dados e croqui.O seu croqui indica que você está usando blocos incompletos balanceados, pelo menos se considerar que o bloco é o retângulo 4 x 4 parcelas. É isso mesmo? Interessante.No entanto, a tabela de dados que você mandou parece estar considerando como bloco a linha que contém 24 parcelas.No meu modo de ver, o bloco com essas dimensões (tão comprido), não deve apresentar a homogeneidade que se espera.É bem provável que a parcela 1C tenha mais semelhança em condições com as parcelas 4B e 8D, por exemplo, do que com as parcelas 33C e 41D.Talvez você possa analisar das duas formas.
De qualquer forma, o modelo para o seu experimento fica conforme código a seguir.
> library(lattice)
>
> url <- "Dados_Conj_2ExpSub.csv"
> tb <- read.csv2(url)
> str(tb)
'data.frame': 96 obs. of 6 variables:
$ FID : Factor w/ 96 levels "12A","12C","13A",..: 9 43 53 87 95 1 3 5 7 11 ...
$ Manejo : Factor w/ 2 levels "PC","PD": 2 2 2 2 2 2 2 2 2 2 ...
$ Rotacao : Factor w/ 6 levels "A","CJ","CO",..: 2 3 5 4 3 6 2 5 1 5 ...
$ Nitrogenio: Factor w/ 2 levels "Com","Sem": 2 2 2 2 2 2 2 2 2 2 ...
$ Bloco : int 1 2 3 4 1 2 3 4 1 2 ...
$ N_kg_ha : int 209 234 209 158 101 96 46 118 130 162 ...
>
> tb <- transform(tb,
+ Bloco = factor(Bloco))
>
> xtabs(~Manejo + Rotacao, data = tb)
Rotacao
Manejo A CJ CO CS P S
PC 8 8 8 8 8 8
PD 8 8 8 8 8 8
> ftable(xtabs(~Manejo + Bloco + Rotacao, data = tb))
Rotacao A CJ CO CS P S
Manejo Bloco
PC 1 2 2 2 2 2 2
2 2 2 2 2 2 2
3 2 2 2 2 2 2
4 2 2 2 2 2 2
PD 1 2 2 2 2 2 2
2 2 2 2 2 2 2
3 2 2 2 2 2 2
4 2 2 2 2 2 2
>
> xyplot(N_kg_ha ~ Rotacao | Manejo,
+ groups = Nitrogenio,
+ type = c("p", "a"),
+ auto.key = TRUE,
+ data = tb)
>
> # Modelo para DBC em diferentes locais:
> # ~ Local + Error(Local/Bloco) + Trat + Local:Trat.
> # Modelo de parcela subdividida em DBC:
> # ~ Bloco + Trat + Error(Bloco:Trat) + Subtrat + Trat:Subtrat
> #
> # Juntanto os dois termos de erro:
> # ~ Error(Local/Bloco) + Error(Bloco:Trat) = Error(Local/Bloco:Trat)
> #
> # Juntanto os termos sistemáticos.
> # ~ Local + Trat + Local:Trat + Subtrat + Trat:Subtrat +
> # Error(Local/Bloco:Trat)
>
> m0 <- aov(N_kg_ha ~
+ Manejo +
+ Manejo/Bloco +
+ Rotacao +
+ Manejo:Rotacao +
+ Nitrogenio +
+ Rotacao:Nitrogenio +
+ Error(Manejo/Bloco:Rotacao),
+ data = tb)
Warning message:
In aov(N_kg_ha ~ Manejo + Manejo/Bloco + Rotacao + Manejo:Rotacao + :
Error() model is singular
>
> summary(m0)
Error: Manejo
Df Sum Sq Mean Sq
Manejo 1 504.2 504.2
Error: Manejo:Bloco:Rotacao
Df Sum Sq Mean Sq F value Pr(>F)
Rotacao 5 16367 3273 1.521 0.213
Manejo:Bloco 6 16966 2828 1.314 0.281
Manejo:Rotacao 5 6213 1243 0.577 0.717
Residuals 30 64567 2152
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Nitrogenio 1 13728 13728 6.561 0.0141 *
Rotacao:Nitrogenio 5 5636 1127 0.539 0.7457
Residuals 42 87874 2092
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> # NOTE: não haverá uma estatística F apropriada para testar o efeito de
> # Manejo já que não existe repetição genuína para os níveis desse fator,
> # assim não existe um erro de extrato para ser o denominador da
> # estatística F.
>
> # Trocando Manejo/Bloco por Manejo:Bloco no termo de Error().
> m0 <- aov(N_kg_ha ~
+ Manejo +
+ Manejo/Bloco +
+ Rotacao +
+ Manejo:Rotacao +
+ Nitrogenio +
+ Rotacao:Nitrogenio +
+ Error(Manejo:Bloco:Rotacao),
+ data = tb)
Warning message:
In aov(N_kg_ha ~ Manejo + Manejo/Bloco + Rotacao + Manejo:Rotacao + :
Error() model is singular
>
> summary(m0)
Error: Manejo:Bloco:Rotacao
Df Sum Sq Mean Sq F value Pr(>F)
Manejo 1 504 504 0.234 0.632
Rotacao 5 16367 3273 1.521 0.213
Manejo:Bloco 6 16966 2828 1.314 0.281
Manejo:Rotacao 5 6213 1243 0.577 0.717
Residuals 30 64567 2152
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Nitrogenio 1 13728 13728 6.561 0.0141 *
Rotacao:Nitrogenio 5 5636 1127 0.539 0.7457
Residuals 42 87874 2092
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> # NOTE: com uma mudança sutil na fórmula, muda-se para uma forma não
> # condizente com a realidade, o delineamento que passa a indicar que as
> # combinações Manejo x Rotação foram casualizados às parcelas dos blocos
> # (fatorial completo na parcela) e dessa forma tem-se uma estatística F
> # para Manejo, mas que não é condizente com a forma como o experimento
> # foi de fato feito. Então está errada essa especificação.
>
> # Declarando o primeiro modelo com a `lme4`.
> library(lme4)
>
> mm0 <- lmer(N_kg_ha ~
+ Manejo +
+ Manejo/Bloco +
+ Rotacao +
+ Manejo:Rotacao +
+ Nitrogenio +
+ Rotacao:Nitrogenio +
+ (1 | Manejo/Bloco:Rotacao),
+ data = tb)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
>
> anova(mm0)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Manejo 1 391.0 391.0 0.1869
Rotacao 5 15911.7 3182.3 1.5210
Nitrogenio 1 13728.2 13728.2 6.5613
Manejo:Bloco 6 16494.0 2749.0 1.3139
Manejo:Rotacao 5 6040.2 1208.0 0.5774
Rotacao:Nitrogenio 5 5636.1 1127.2 0.5387
> VarCorr(mm0)
Groups Name Std.Dev.
Bloco:Rotacao:Manejo (Intercept) 5.4722
Manejo (Intercept) 3.3708
Residual 45.7415
À disposição.Walmes.
_______________________________________________
R-br mailing list
R-br@listas.c3sl.ufpr.br
https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Henrique, Bom que você conferiu reprodutibilidade MAS o código de conduta da lista pede para que não sejam enviados anexos para a lista. Existem outras formas de conferir acesso aos dados sem que seja por anexo, ok?
As dúvidas, além de saber se estou ou não num caminho lúcido em termos da estatística: 1) Em alguns casos houve violação de alguma pressuposição (estou testando norm. dos residuos e homogeneidade apenas), como para a variável resposta 3, no local II colhido na idade 8. - Procedo com alguma transformação ? tipo utilizar boxcox e ver qual melhor forma de transformar os dados ? - Ou Parto para outro tipo de abordagem ? alguma sugestão ?
Nem uma coisa nem outra. Segue com a análise usual. Por que? Bem, para começar, você sabe que, se o nível de significância é de 5%, 1 em cada 20 análises (em média) irá rejeitar a hipótese nula mesmo que ela seja verdadeira. É o erro tipo I. Então, o que está acontecendo com você é basicamente isso. De tanto fazer testes, uma hora haverá uma rejeição que pode ser apenas casual. Não rodei seu script para saber se de fato a rejeição deve-se a alguma anomalia grave. Mas de qualquer forma, fica aí o raciocínio. Outro ponto, eu não gosto dessa abordagem "pilha de testes de hipóteses". Prefiro, ainda que seja considerado subjetiva, a análise gráfica dos resíduos porque a leitura conjunta dos gráficos subsidia alguma ação em caso de necessidade, coisa que o teste não faz porque ele diz apenas "sim" e "não". Um último ponto, sendo a mesma variável resposta vista em divisões dos dados (por local, por época, etc), não seria desejável aplicar uma transformação diferente em cada fração. O melhor seria encontrar uma transformação que seja adequada para todas as frações. Você pode examinar isso graficamente, fazer uma boa análise exploratória e depois de ajustar o modelo, fazer uma boa análise dos resíduos. A `VariavelResposta1`, por exemplo, mostra um aumento de variabilidade em função da idade, coisa que é comum em dados de crescimento (não se é o caso). No início as unidades experimentais de um mesmo tratamento são mais parecidas, mas com o passar do tempo na exposição aos várias fatores não controlados, as diferenças nas variáveis resposta tendem a aumentar, aumentando portanto a variância do erro. Talvez uma transformação log, raíz, etc, resolva isso. Mas nada impede de modelar a variância conjuntamente com a média (sem querer complicar).
2) Há alguma forma de análise considerando as todas as idades?
Se você considerar representar efeito do fator idade (fator quantitativo) por alguma função, por exemplo, polinômio, então mesmo não tendo as mesmas idades (mesmo sendo um fatorial incompleto) nas cultivares, você vai analisar tudo conjuntamente com facilidade. Se assumir idade como variável quantitativa (o que não é recomendável pela própria natureza da variável), aí você terá um fatorial incompleto que é coisa complicada de analisar por causa das celas ausentes.
3) Ao menos duas idades e 3 dos cultivares são comuns em ambos os locais, seria muito abusado tentar realizar uma análise conjunta apenas para essas cultivares? (com abordagem similar ao proposto para o caso do José)
Poderia ser feito. A idade não é preocupação se for usado regressão nela. Mas a análise exploratória que fiz indica que só existe uma época comum aos dois locais. Dessa forma não tem como declara interação local:epoca, pois pela ausência de mais combinações comuns, o efeito de interação não é estimável se época for considerado fator qualitativo. Se for quantitativo, assumindo uma função, então daria para estimar/representar a interação. Tudo isso vai do seu interesse e disposição em fazer a análise. Ao meu ver, o experimento parece não ter sido planejado para uma análise conjunta visando estimar interação etc. O jeito mais tranquilo é analisar os locais separados, até mesmo porque, não se tem como inferir sobre o efeito de local (não creio que essa seja a hipótese), e uma vez que estimar a interação local:epoca é complicado neste caso, não sei se compensa o esforço. A interação local:cultivar dá para estimar, mas apenas baseando-se em 3 cultivares comuns aos locais, mas tendo em mente que as épocas de cultivo foram diferentes (apenas uma comum), o que certamente incluirá um complicador para gerar confundimento. Já as idades, como em cada experimento elas são as mesmas, e o padrão é bem linear, você pode começar com o modelo linear simples pro efeito de idade e ir aumentando o grau conforme necessidade. Segue código usado para explorar os dados. # Lê direto da planilha. da <- gdata::read.xls("db_r-list.xlsx", stringsAsFactors = FALSE) str(da) factor_roman <- function(x) { u <- sort(as.roman(unique(x))) factor(x, levels = as.character(u)) } # Converte para fator com níveis romanos ordenados. da <- transform(da, Epoca = factor_roman(Epoca), Cultivar = factor_roman(Cultivar), Local = factor_roman(Local), Bloco = factor(Repeticao)) # Tabelas de frequência dos pontos experimentais. xtabs(~Epoca + Local, data = da) xtabs(~Cultivar + Local, data = da) xtabs(~Cultivar + Epoca, data = da) ftable(xtabs(~Local + Idade + Cultivar, data = da)) library(lattice) # Diagrama de dispersão. xyplot(VariavelResposta1 ~ Idade | Epoca + Local, groups = Cultivar, data = da, type = c("p", "a")) # Separando por local. gridExtra::grid.arrange( ncol = 1, xyplot(VariavelResposta1 ~ Idade | Epoca, groups = Cultivar, data = da, subset = Local == "I", layout = c(NA, 1), type = c("p", "a")), xyplot(VariavelResposta1 ~ Idade | Epoca, groups = Cultivar, data = da, subset = Local == "II", layout = c(NA, 1), type = c("p", "a"))) # Cultivares comuns aos dois locais. u <- as.array(xtabs(~Cultivar + Local, data = da)) i <- apply(u, 1, prod) > 0 common <- as.data.frame(u[i, ]) common$Freq <- NULL common # Pega apenas as combinações presentes. db <- droplevels(merge(common, da, all.x = TRUE, all.y = FALSE)) summary(db) # Diagrama de dispersão para a seleção. xyplot(VariavelResposta1 ~ Idade | Epoca + Local, groups = Cultivar, data = db, type = c("p", "a")) xyplot(VariavelResposta1 ~ Idade | Cultivar + Local, groups = Epoca, data = db, type = c("p", "a")) À disposição. Walmes.
Prezado Walmes, Primeiramente perdao pela falha. Faco parte da lista desde uns 3 anos atras, quando comecei a usar o R de fato.No inicio foi mais por manter antenado nas possibilidades e aprender com duvidas de colegas que usam ou usaram esse canal.O email do Jose me encorajou a enviar um pedido de ajuda.Comumente vejo emails sobre a questao de reprodutibilidade e dei muita atencao a isso e nao em uma ponto basico sobre o codigo de conduta que li la tras qdo entrei na lista, o que deveria ter percebido pq nao lembro de ter visto alguem anexando nada. Enfim, aprendido. Segundo, obrigado pelo compartilhamento da sua visao e do codigo. Com certeza elucidou muita coisa pra mim. Ate maisHenrique Em quarta-feira, 2 de outubro de 2019 10:05:37 GMT-4, Walmes Zeviani <walmeszeviani@gmail.com> escreveu: Henrique, Bom que você conferiu reprodutibilidade MAS o código de conduta da lista pede para que não sejam enviados anexos para a lista.Existem outras formas de conferir acesso aos dados sem que seja por anexo, ok?
As dúvidas, além de saber se estou ou não num caminho lúcido em termos da estatística: 1) Em alguns casos houve violação de alguma pressuposição (estou testando norm. dos residuos e homogeneidade apenas),> como para a variável resposta 3, no local II colhido na idade 8. - Procedo com alguma transformação ? tipo utilizar boxcox e ver qual melhor forma de transformar os dados ? - Ou Parto para outro tipo de abordagem ? alguma sugestão ? Nem uma coisa nem outra. Segue com a análise usual. Por que? Bem, para começar, você sabe que, se o nível de significância é de 5%, 1 em cada 20 análises (em média) irá rejeitar a hipótese nula mesmo que ela seja verdadeira. É o erro tipo I. Então, o que está acontecendo com você é basicamente isso. De tanto fazer testes, uma hora haverá uma rejeição que pode ser apenas casual. Não rodei seu script para saber se de fato a rejeição deve-se a alguma anomalia grave. Mas de qualquer forma, fica aí o raciocínio. Outro ponto, eu não gosto dessa abordagem "pilha de testes de hipóteses". Prefiro, ainda que seja considerado subjetiva, a análise gráfica dos resíduos porque a leitura conjunta dos gráficos subsidia alguma ação em caso de necessidade, coisa que o teste não faz porque ele diz apenas "sim" e "não". Um último ponto, sendo a mesma variável resposta vista em divisões dos dados (por local, por época, etc), não seria desejável aplicar uma transformação diferente em cada fração. O melhor seria encontrar uma transformação que seja adequada para todas as frações. Você pode examinar isso graficamente, fazer uma boa análise exploratória e depois de ajustar o modelo, fazer uma boa análise dos resíduos. A `VariavelResposta1`, por exemplo, mostra um aumento de variabilidade em função da idade, coisa que é comum em dados de crescimento (não se é o caso). No início as unidades experimentais de um mesmo tratamento são mais parecidas, mas com o passar do tempo na exposição aos várias fatores não controlados, as diferenças nas variáveis resposta tendem a aumentar, aumentando portanto a variância do erro. Talvez uma transformação log, raíz, etc, resolva isso. Mas nada impede de modelar a variância conjuntamente com a média (sem querer complicar).
2) Há alguma forma de análise considerando as todas as idades? Se você considerar representar efeito do fator idade (fator quantitativo) por alguma função, por exemplo, polinômio, então mesmo não tendo as mesmas idades (mesmo sendo um fatorial incompleto) nas cultivares, você vai analisar tudo conjuntamente com facilidade. Se assumir idade como variável quantitativa (o que não é recomendável pela própria natureza da variável), aí você terá um fatorial incompleto que é coisa complicada de analisar por causa das celas ausentes.
3) Ao menos duas idades e 3 dos cultivares são comuns em ambos os locais, seria muito abusado tentar realizar uma análise conjunta> apenas para essas cultivares? (com abordagem similar ao proposto para o caso do José) Poderia ser feito. A idade não é preocupação se for usado regressão nela. Mas a análise exploratória que fiz indica que só existe uma época comum aos dois locais. Dessa forma não tem como declara interação local:epoca, pois pela ausência de mais combinações comuns, o efeito de interação não é estimável se época for considerado fator qualitativo. Se for quantitativo, assumindo uma função, então daria para estimar/representar a interação. Tudo isso vai do seu interesse e disposição em fazer a análise. Ao meu ver, o experimento parece não ter sido planejado para uma análise conjunta visando estimar interação etc. O jeito mais tranquilo é analisar os locais separados, até mesmo porque, não se tem como inferir sobre o efeito de local (não creio que essa seja a hipótese), e uma vez que estimar a interação local:epoca é complicado neste caso, não sei se compensa o esforço. A interação local:cultivar dá para estimar, mas apenas baseando-se em 3 cultivares comuns aos locais, mas tendo em mente que as épocas de cultivo foram diferentes (apenas uma comum), o que certamente incluirá um complicador para gerar confundimento. Já as idades, como em cada experimento elas são as mesmas, e o padrão é bem linear, você pode começar com o modelo linear simples pro efeito de idade e ir aumentando o grau conforme necessidade.
Segue código usado para explorar os dados. # Lê direto da planilha. da <- gdata::read.xls("db_r-list.xlsx", stringsAsFactors = FALSE) str(da) factor_roman <- function(x) { u <- sort(as.roman(unique(x))) factor(x, levels = as.character(u)) } # Converte para fator com níveis romanos ordenados. da <- transform(da, Epoca = factor_roman(Epoca), Cultivar = factor_roman(Cultivar), Local = factor_roman(Local), Bloco = factor(Repeticao)) # Tabelas de frequência dos pontos experimentais. xtabs(~Epoca + Local, data = da) xtabs(~Cultivar + Local, data = da) xtabs(~Cultivar + Epoca, data = da) ftable(xtabs(~Local + Idade + Cultivar, data = da)) library(lattice) # Diagrama de dispersão. xyplot(VariavelResposta1 ~ Idade | Epoca + Local, groups = Cultivar, data = da, type = c("p", "a")) # Separando por local. gridExtra::grid.arrange( ncol = 1, xyplot(VariavelResposta1 ~ Idade | Epoca, groups = Cultivar, data = da, subset = Local == "I", layout = c(NA, 1), type = c("p", "a")), xyplot(VariavelResposta1 ~ Idade | Epoca, groups = Cultivar, data = da, subset = Local == "II", layout = c(NA, 1), type = c("p", "a"))) # Cultivares comuns aos dois locais. u <- as.array(xtabs(~Cultivar + Local, data = da)) i <- apply(u, 1, prod) > 0 common <- as.data.frame(u[i, ]) common$Freq <- NULL common # Pega apenas as combinações presentes. db <- droplevels(merge(common, da, all.x = TRUE, all.y = FALSE)) summary(db) # Diagrama de dispersão para a seleção. xyplot(VariavelResposta1 ~ Idade | Epoca + Local, groups = Cultivar, data = db, type = c("p", "a")) xyplot(VariavelResposta1 ~ Idade | Cultivar + Local, groups = Epoca, data = db, type = c("p", "a")) À disposição.Walmes.
participantes (3)
-
Henrique Boriolo Dias -
José Lucas Safanelli -
Walmes Zeviani