Importação e análise exploratória
#-----------------------------------------------------------------------
# Pacotes.
# Instalar pacotes antes.
library(tidyverse) # Manipulação e visualização de dados.
library(emmeans) # Médias marginais e desdobramentos.#-----------------------------------------------------------------------
# Leitura e preparo.
tb <- read_tsv("dados.txt", locale = locale(decimal_mark = ","))
attr(tb, "spec") <- NULL
head(tb)## # A tibble: 6 x 11
## iba sexo viva calo morta folha broto comprimento numero enraizada
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 Macho 0 0 1 0 0 0 0 0
## 2 0 Macho 1 0 0 1 1 0 0 0
## 3 0 Macho 0 0 1 0 0 0 0 0
## 4 0 Macho 0 0 1 0 0 0 0 0
## 5 0 Macho 0 0 1 0 0 0 0 0
## 6 0 Macho 0 0 1 0 0 0 0 0
## # … with 1 more variable: rept <chr># Combinação dos níveis dos fatores.
tb %>%
xtabs(~sexo + iba, .)## iba
## sexo 0 1500 3000 4500 6000
## Fêmea 64 64 64 64 64
## Macho 64 64 64 64 64# Desce preenchendo a informação da repetição.
tb <- tb %>%
mutate(sexo = factor(sexo),
iba = iba/1000) %>%
fill(rept)
# Estrutura.
str(tb)## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 640 obs. of 11 variables:
## $ iba : num 0 0 0 0 0 0 0 0 0 0 ...
## $ sexo : Factor w/ 2 levels "Fêmea","Macho": 2 2 2 2 2 2 2 2 2 2 ...
## $ viva : num 0 1 0 0 0 0 0 1 1 1 ...
## $ calo : num 0 0 0 0 0 0 0 0 0 0 ...
## $ morta : num 1 0 1 1 1 1 1 0 0 0 ...
## $ folha : num 0 1 0 0 0 0 0 1 0 1 ...
## $ broto : num 0 1 0 0 0 0 0 0 0 0 ...
## $ comprimento: num 0 0 0 0 0 0 0 0 0 0 ...
## $ numero : num 0 0 0 0 0 0 0 0 0 0 ...
## $ enraizada : num 0 0 0 0 0 0 0 0 0 0 ...
## $ rept : chr "R1" "R1" "R1" "R1" ...#-----------------------------------------------------------------------
# Obtém agregação.
# Faz agregação por unidade experimental (= 16 estacas).
tbu <- tb %>%
group_by(iba, sexo, rept) %>%
mutate(n = 1L) %>%
summarise_at(vars(viva:n), sum) %>%
ungroup()
str(tbu)## Classes 'tbl_df', 'tbl' and 'data.frame': 40 obs. of 12 variables:
## $ iba : num 0 0 0 0 0 0 0 0 1.5 1.5 ...
## $ sexo : Factor w/ 2 levels "Fêmea","Macho": 1 1 1 1 2 2 2 2 1 1 ...
## $ rept : chr "R1" "R2" "R3" "R4" ...
## $ viva : num 0 0 2 1 6 6 0 1 0 0 ...
## $ calo : num 0 0 2 1 1 0 0 0 0 0 ...
## $ morta : num 16 16 12 14 10 10 16 15 16 15 ...
## $ folha : num 0 0 0 1 4 4 0 0 0 0 ...
## $ broto : num 0 0 2 2 2 4 0 1 0 1 ...
## $ comprimento: num 0 0 0 0 0 0 0 0 0 0 ...
## $ numero : num 0 0 0 0 0 0 0 0 0 0 ...
## $ enraizada : num 0 0 0 0 0 0 0 0 0 0 ...
## $ n : int 16 16 16 16 16 16 16 16 16 16 ...# Empilha as respostas para fazer gráfico com todas elas.
tbs <- tbu %>%
mutate_at(vars(viva:enraizada), ~ ./n) %>%
gather(key = "variable", value = "value", viva:enraizada) %>%
mutate(variable = factor(variable,
levels = c("morta", "viva", "calo",
"broto", "folha", "enraizada",
"numero", "comprimento")))
#-----------------------------------------------------------------------
# Gráficos.
ggplot(tbs, aes(x = iba,
y = value,
group = sexo,
color = sexo)) +
facet_wrap(facets = ~variable, scale = "free_y") +
geom_jitter(width = 0.1, height = 0, pch = 1) +
stat_summary(geom = "line", fun.y = "mean")
#-----------------------------------------------------------------------Análise do experimento
Estacas mortas
Para as variáveis binárias considerou-se o modelo quasi-binomial. Este modelo é baseado no modelo binomial com a diferença de que é estimado um parâmetro de dispersão, tornando o modelo mais flexível para acomodar o comum cenário de superdispersão.
O modelo considerado é \[ \begin{align*} y_{ijk} &\sim \text{Quasi-Binomial}(p_{ij}, \phi, n_{ijk})\\ \text{logit}(p_{ij}) &= \mu + \tau_i + \lambda_j + \theta_{ij} &\quad (1)\\ \text{logit}(p_{ij}) &= \mu + \tau_i + \beta_1\text{iba}_j + \beta_{1i}\text{iba}_j &\quad (2)\\ \end{align*} \] em que:
- \(y_{ijk}\) é o número observado de estacas na condição considerada sucesso, e.g. morrer, viver, apresentar calo, etc.
- \(n_{ijk}\) é o total de estacas da unidade experimental, no caso 16 em todas elas.
- \(p_{ij}\) é a probabilidade de sucesso para a condição experimental \(ij\).
- \(\phi\) é o parâmetro de dispersão.
- \(\text{logit(.)}\) é a função de ligação logit.
- \(\mu\) é uma constante comum à todas as celas experimentais.
- \(\tau_i\) (\(i = 1, 2\)) acomoda o efeito principal de sexo.
- \(\lambda_j\) (\(j = 1, \ldots, 5\)) acomoda o efeito principal de iba, considerado como fator qualitativo no modelo com preditor (1).
- \(\theta_{ij}\) acomoda o efeito de interação entre sexo e iba.
- \(\beta_1\) acomoda o efeito linear de iba, considerado como fator quantitativo no modelo com preditor (2).
- \(\beta_{1i}\) acomoda o efeito de interação entre sexo e iba com termos linear por nível de sexo.
O modelo foi ajustado aos dados por mínimos quadrados ponderados interativo. Todas as análises foram feitas no software R de computação estatística (R Core Team, 2019).
R Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
# Ajuste do modelo.
m0 <- glm(cbind(morta, n - morta) ~ sexo * factor(iba),
data = tbu,
family = quasibinomial)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(morta, n - morta)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 91.962
## sexo 1 0.2150 38 91.747 0.0950 0.7601
## factor(iba) 4 4.7014 34 87.045 0.5192 0.7222
## sexo:factor(iba) 4 7.6566 30 79.389 0.8456 0.5075# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: cbind(morta, n - morta) ~ iba
## Model 2: cbind(morta, n - morta) ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 38 91.532
## 2 30 79.389 8 12.143 0.6705 0.7131# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = cbind(morta, n - morta) ~ iba, family = quasibinomial,
## data = tbu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.4121 -0.7486 0.1448 1.3062 2.2566
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.75736 0.28565 6.152 3.52e-07 ***
## iba 0.03586 0.07986 0.449 0.656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 2.129602)
##
## Null deviance: 91.962 on 39 degrees of freedom
## Residual deviance: 91.532 on 38 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Médias estimadas marginais para as combinações experimentais conforme
# o modelo maximal.
emmeans(m0,
specs = ~sexo + iba,
type = "response")## sexo iba prob SE df asymp.LCL asymp.UCL
## Fêmea 0.0 0.906 0.0548 Inf 0.732 0.972
## Macho 0.0 0.797 0.0757 Inf 0.611 0.907
## Fêmea 1.5 0.797 0.0757 Inf 0.611 0.907
## Macho 1.5 0.859 0.0654 Inf 0.679 0.946
## Fêmea 3.0 0.875 0.0622 Inf 0.697 0.955
## Macho 3.0 0.953 0.0398 Inf 0.780 0.991
## Fêmea 4.5 0.891 0.0587 Inf 0.714 0.964
## Macho 4.5 0.859 0.0654 Inf 0.679 0.946
## Fêmea 6.0 0.891 0.0587 Inf 0.714 0.964
## Macho 6.0 0.828 0.0710 Inf 0.645 0.928
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba prob SE df asymp.LCL asymp.UCL
## 0.0 0.853 0.0358 Inf 0.768 0.910
## 1.5 0.860 0.0244 Inf 0.805 0.901
## 3.0 0.866 0.0197 Inf 0.822 0.900
## 4.5 0.872 0.0238 Inf 0.818 0.912
## 6.0 0.878 0.0323 Inf 0.799 0.928
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# A equação ajustada que gerou os valores acima.
substitute(1/(1 + exp(-(b0 + b1 * iba))),
setNames(as.list(coef(m1)), c("b0", "b1")))## 1/(1 + exp(-(1.75735546277956 + 0.0358588903104317 * iba)))# Malha para a predição e construção do gráfico.
grid <- with(tb,
crossing(#sexo = levels(sexo),
iba = seq(min(iba),
max(iba),
length.out = 71)))
grid$morta <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.95, 0.05)
ggplot(data = tbu,
mapping = aes(x = iba,
y = morta/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = morta)) +
labs(x = "Concentração de IBA (milhar)",
y = "Proporção de estacas mortas",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Estacas vivas
# Ajuste do modelo.
m0 <- glm(cbind(viva, n - viva) ~ sexo * factor(iba),
data = tbu,
family = quasibinomial)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(viva, n - viva)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 81.896
## sexo 1 2.0332 38 79.863 1.1123 0.3000
## factor(iba) 4 5.0887 34 74.774 0.6959 0.6007
## sexo:factor(iba) 4 9.1207 30 65.653 1.2474 0.3123# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: cbind(viva, n - viva) ~ iba
## Model 2: cbind(viva, n - viva) ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 38 79.611
## 2 30 65.653 8 13.957 0.9544 0.4887# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = cbind(viva, n - viva) ~ iba, family = quasibinomial,
## data = tbu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9248 -1.6640 -0.3082 0.7913 2.7760
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.09762 0.31296 -6.703 6.23e-08 ***
## iba -0.10204 0.09323 -1.094 0.281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.882562)
##
## Null deviance: 81.896 on 39 degrees of freedom
## Residual deviance: 79.611 on 38 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Médias estimadas marginais para as combinações experimentais conforme
# o modelo maximal.
emmeans(m0,
specs = ~sexo + iba,
type = "response")## sexo iba prob SE df asymp.LCL asymp.UCL
## Fêmea 0.0 0.0469 0.0357 Inf 0.01016 0.191
## Macho 0.0 0.2031 0.0680 Inf 0.10063 0.367
## Fêmea 1.5 0.0781 0.0454 Inf 0.02407 0.226
## Macho 1.5 0.0938 0.0493 Inf 0.03214 0.244
## Fêmea 3.0 0.0781 0.0454 Inf 0.02407 0.226
## Macho 3.0 0.0312 0.0294 Inf 0.00478 0.178
## Fêmea 4.5 0.0625 0.0409 Inf 0.01668 0.208
## Macho 4.5 0.1250 0.0559 Inf 0.04986 0.280
## Fêmea 6.0 0.0781 0.0454 Inf 0.02407 0.226
## Macho 6.0 0.0469 0.0357 Inf 0.01016 0.191
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba prob SE df asymp.LCL asymp.UCL
## 0.0 0.1093 0.0305 Inf 0.0623 0.185
## 1.5 0.0953 0.0191 Inf 0.0639 0.140
## 3.0 0.0829 0.0151 Inf 0.0577 0.118
## 4.5 0.0720 0.0175 Inf 0.0443 0.115
## 6.0 0.0624 0.0217 Inf 0.0312 0.121
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# A equação ajustada que gerou os valores acima.
substitute(1/(1 + exp(-(b0 + b1 * iba))),
setNames(as.list(coef(m1)), c("b0", "b1")))## 1/(1 + exp(-(-2.09761655953521 + -0.102038492730299 * iba)))# Malha para a predição e construção do gráfico.
grid$viva <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.95, 0.95)
ggplot(data = tbu,
mapping = aes(x = iba,
y = viva/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = viva)) +
labs(x = "Concentração de IBA (milhar)",
y = "Proporção de estacas vivas",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Estacas com folhas
# Ajuste do modelo.
m0 <- glm(cbind(folha, n - folha) ~ sexo * factor(iba),
data = tbu,
family = quasibinomial)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(folha, n - folha)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 76.401
## sexo 1 1.6178 38 74.783 0.9397 0.3401
## factor(iba) 4 5.9584 34 68.824 0.8652 0.4961
## sexo:factor(iba) 4 8.6336 30 60.191 1.2537 0.3098# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: cbind(folha, n - folha) ~ iba
## Model 2: cbind(folha, n - folha) ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 38 76.311
## 2 30 60.191 8 16.12 1.1704 0.3487# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = cbind(folha, n - folha) ~ iba, family = quasibinomial,
## data = tbu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7441 -1.6533 -0.3815 0.5261 3.2188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.4247 0.3424 -7.081 1.91e-08 ***
## iba 0.0199 0.0917 0.217 0.829
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.900235)
##
## Null deviance: 76.401 on 39 degrees of freedom
## Residual deviance: 76.311 on 38 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Médias estimadas marginais para as combinações experimentais conforme
# o modelo maximal.
emmeans(m0,
specs = ~sexo + iba,
type = "response")## sexo iba prob SE df asymp.LCL asymp.UCL
## Fêmea 0.0 0.0156 0.0203 Inf 0.00119 0.175
## Macho 0.0 0.1250 0.0542 Inf 0.05127 0.274
## Fêmea 1.5 0.1406 0.0570 Inf 0.06095 0.292
## Macho 1.5 0.1094 0.0512 Inf 0.04200 0.256
## Fêmea 3.0 0.0625 0.0397 Inf 0.01736 0.201
## Macho 3.0 0.0312 0.0285 Inf 0.00506 0.170
## Fêmea 4.5 0.0469 0.0347 Inf 0.01063 0.184
## Macho 4.5 0.1250 0.0542 Inf 0.05127 0.274
## Fêmea 6.0 0.0938 0.0478 Inf 0.03320 0.238
## Macho 6.0 0.1094 0.0512 Inf 0.04200 0.256
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba prob SE df asymp.LCL asymp.UCL
## 0.0 0.0813 0.0256 Inf 0.0433 0.148
## 1.5 0.0836 0.0185 Inf 0.0537 0.128
## 3.0 0.0859 0.0153 Inf 0.0603 0.121
## 4.5 0.0883 0.0189 Inf 0.0576 0.133
## 6.0 0.0907 0.0273 Inf 0.0495 0.160
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# A equação ajustada que gerou os valores acima.
substitute(1/(1 + exp(-(b0 + b1 * iba))),
setNames(as.list(coef(m1)), c("b0", "b1")))## 1/(1 + exp(-(-2.42471047238679 + 0.0198980624818498 * iba)))# Malha para a predição e construção do gráfico.
grid$folha <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.95, 0.95)
ggplot(data = tbu,
mapping = aes(x = iba,
y = folha/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = folha)) +
labs(x = "Concentração de IBA (milhar)",
y = "Proporção de estacas com folhas",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Estacas com calos
# Ajuste do modelo.
m0 <- glm(cbind(calo, n - calo) ~ sexo * factor(iba),
data = tbu,
family = quasibinomial)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(calo, n - calo)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 62.375
## sexo 1 1.0363 38 61.339 0.7803 0.3841
## factor(iba) 4 20.4822 34 40.857 3.8554 0.0121 *
## sexo:factor(iba) 4 1.4563 30 39.401 0.2741 0.8923
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: cbind(calo, n - calo) ~ iba
## Model 2: cbind(calo, n - calo) ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 38 52.444
## 2 30 39.401 8 13.043 1.2275 0.3173# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = cbind(calo, n - calo) ~ iba, family = quasibinomial,
## data = tbu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4055 -0.8257 -0.5533 -0.4042 4.5285
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.7540 0.5096 -5.404 3.74e-06 ***
## iba -0.4202 0.2275 -1.847 0.0726 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 2.338452)
##
## Null deviance: 62.375 on 39 degrees of freedom
## Residual deviance: 52.444 on 38 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Médias estimadas marginais para as combinações experimentais conforme
# o modelo maximal.
emmeans(m0,
specs = ~sexo + iba,
type = "response")## sexo iba prob SE df asymp.LCL asymp.UCL
## Fêmea 0.0 0.0469 3.04e-02 Inf 0.01277 0.158
## Macho 0.0 0.0156 1.79e-02 Inf 0.00163 0.134
## Fêmea 1.5 0.0938 4.20e-02 Inf 0.03778 0.214
## Macho 1.5 0.0469 3.04e-02 Inf 0.01277 0.158
## Fêmea 3.0 0.0156 1.79e-02 Inf 0.00163 0.134
## Macho 3.0 0.0312 2.51e-02 Inf 0.00633 0.140
## Fêmea 4.5 0.0000 1.83e-06 Inf 0.00000 1.000
## Macho 4.5 0.0000 1.83e-06 Inf 0.00000 1.000
## Fêmea 6.0 0.0000 1.83e-06 Inf 0.00000 1.000
## Macho 6.0 0.0000 1.83e-06 Inf 0.00000 1.000
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba prob SE df asymp.LCL asymp.UCL
## 0.0 0.05986 0.02868 Inf 0.022915 0.1474
## 1.5 0.03279 0.01239 Inf 0.015524 0.0679
## 3.0 0.01773 0.00919 Inf 0.006378 0.0483
## 4.5 0.00952 0.00753 Inf 0.002007 0.0439
## 6.0 0.00509 0.00562 Inf 0.000582 0.0430
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# A equação ajustada que gerou os valores acima.
substitute(1/(1 + exp(-(b0 + b1 * iba))),
setNames(as.list(coef(m1)), c("b0", "b1")))## 1/(1 + exp(-(-2.75397448474482 + -0.42016514230657 * iba)))# Malha para a predição e construção do gráfico.
grid$calo <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.95, 0.95)
ggplot(data = tbu,
mapping = aes(x = iba,
y = calo/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = calo)) +
labs(x = "Concentração de IBA (milhar)",
y = "Proporção de estacas com calo",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Estacas com broto
# Ajuste do modelo.
m0 <- glm(cbind(broto, n - broto) ~ sexo * factor(iba),
data = tbu,
family = quasibinomial)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(broto, n - broto)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 66.248
## sexo 1 0.3910 38 65.857 0.3030 0.5861
## factor(iba) 4 9.1157 34 56.741 1.7661 0.1618
## sexo:factor(iba) 4 10.8336 30 45.908 2.0989 0.1057# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: cbind(broto, n - broto) ~ iba
## Model 2: cbind(broto, n - broto) ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 38 59.689
## 2 30 45.908 8 13.781 1.335 0.2648# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = cbind(broto, n - broto) ~ iba, family = quasibinomial,
## data = tbu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9247 -1.2670 -0.0298 0.4917 3.2083
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.09775 0.28629 -7.328 8.9e-09 ***
## iba -0.19323 0.09394 -2.057 0.0466 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.476975)
##
## Null deviance: 66.248 on 39 degrees of freedom
## Residual deviance: 59.689 on 38 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Médias estimadas marginais para as combinações experimentais conforme
# o modelo maximal.
emmeans(m0,
specs = ~sexo + iba,
type = "response")## sexo iba prob SE df asymp.LCL asymp.UCL
## Fêmea 0.0 0.0625 0.0344 Inf 0.02068 0.174
## Macho 0.0 0.1094 0.0443 Inf 0.04793 0.231
## Fêmea 1.5 0.1562 0.0516 Inf 0.07923 0.285
## Macho 1.5 0.0781 0.0381 Inf 0.02917 0.193
## Fêmea 3.0 0.1094 0.0443 Inf 0.04793 0.231
## Macho 3.0 0.0156 0.0176 Inf 0.00168 0.130
## Fêmea 4.5 0.0156 0.0176 Inf 0.00168 0.130
## Macho 4.5 0.0781 0.0381 Inf 0.02917 0.193
## Fêmea 6.0 0.0312 0.0247 Inf 0.00647 0.138
## Macho 6.0 0.0312 0.0247 Inf 0.00647 0.138
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba prob SE df asymp.LCL asymp.UCL
## 0.0 0.1093 0.0279 Inf 0.0654 0.1770
## 1.5 0.0841 0.0157 Inf 0.0580 0.1205
## 3.0 0.0643 0.0122 Inf 0.0442 0.0928
## 4.5 0.0489 0.0132 Inf 0.0287 0.0823
## 6.0 0.0371 0.0142 Inf 0.0173 0.0776
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# A equação ajustada que gerou os valores acima.
substitute(1/(1 + exp(-(b0 + b1 * iba))),
setNames(as.list(coef(m1)), c("b0", "b1")))## 1/(1 + exp(-(-2.09775447778747 + -0.193234200108692 * iba)))# Malha para a predição e construção do gráfico.
grid$broto <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.95, 0.95)
ggplot(data = tbu,
mapping = aes(x = iba,
y = broto/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = broto)) +
labs(x = "Concentração de IBA (milhar)",
y = "Proporção de estacas com broto",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Estacas enraizadas
# Ajuste do modelo.
m0 <- glm(cbind(enraizada, n - enraizada) ~ sexo * factor(iba),
data = tbu,
family = quasibinomial)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(enraizada, n - enraizada)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 45.841
## sexo 1 0.4902 38 45.351 0.6549 0.424723
## factor(iba) 4 16.2496 34 29.101 5.4281 0.002071 **
## sexo:factor(iba) 4 5.0845 30 24.017 1.6984 0.176408
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: cbind(enraizada, n - enraizada) ~ iba
## Model 2: cbind(enraizada, n - enraizada) ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 38 31.547
## 2 30 24.017 8 7.5301 1.2577 0.3018# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = cbind(enraizada, n - enraizada) ~ iba, family = quasibinomial,
## data = tbu)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.59105 -0.79081 -0.38865 -0.04462 1.45096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.3535 0.6644 -8.058 9.56e-10 ***
## iba 0.4761 0.1336 3.562 0.00101 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.8590547)
##
## Null deviance: 45.841 on 39 degrees of freedom
## Residual deviance: 31.547 on 38 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Médias estimadas marginais para as combinações experimentais conforme
# o modelo maximal.
emmeans(m0,
specs = ~sexo + iba,
type = "response")## sexo iba prob SE df asymp.LCL asymp.UCL
## Fêmea 0.0 0.0000 1.38e-06 Inf 0.00000 1.0000
## Macho 0.0 0.0000 1.38e-06 Inf 0.00000 1.0000
## Fêmea 1.5 0.0156 1.34e-02 Inf 0.00287 0.0806
## Macho 1.5 0.0156 1.34e-02 Inf 0.00287 0.0806
## Fêmea 3.0 0.0312 1.88e-02 Inf 0.00945 0.0983
## Macho 3.0 0.0156 1.34e-02 Inf 0.00287 0.0806
## Fêmea 4.5 0.0469 2.29e-02 Inf 0.01772 0.1182
## Macho 4.5 0.0156 1.34e-02 Inf 0.00287 0.0806
## Fêmea 6.0 0.0312 1.88e-02 Inf 0.00945 0.0983
## Macho 6.0 0.1250 3.58e-02 Inf 0.06999 0.2133
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba prob SE df asymp.LCL asymp.UCL
## 0.0 0.00471 0.00311 Inf 0.00129 0.0171
## 1.5 0.00957 0.00455 Inf 0.00376 0.0241
## 3.0 0.01935 0.00597 Inf 0.01054 0.0353
## 4.5 0.03875 0.00819 Inf 0.02552 0.0584
## 6.0 0.07606 0.01964 Inf 0.04543 0.1246
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# A equação ajustada que gerou os valores acima.
substitute(1/(1 + exp(-(b0 + b1 * iba))),
setNames(as.list(coef(m1)), c("b0", "b1")))## 1/(1 + exp(-(-5.35354649867744 + 0.476068749113278 * iba)))# Malha para a predição e construção do gráfico.
grid$enraizada <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.05, 0.95)
ggplot(data = tbu,
mapping = aes(x = iba,
y = enraizada/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = enraizada)) +
labs(x = "Concentração de IBA (milhar)",
y = "Proporção de estacas enraizadas",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Número de raízes
O número de raízes por estaca é uma variável de contagem, tipo para o qual habitualmente assume-se a distribuição Poisson. Será considerado a especificação com o modelo Quasi-Poisson que é mais flexível por ter uma parâmetro de dispersão. Como os valores de contagem foram agregados usando a média para a unidade experimental, serão utilizados os tamanhos amostrais como peso.
O modelo associado a essa reposta é \[ \begin{align*} y_{ijk} &\sim \text{Quasi-Poisson}(\mu_{ij}, \phi, n_{ijk})\\ \log(\mu_{ij}) &= \mu + \tau_i + \lambda_j + \theta_{ij} &\quad (1)\\ \log(\mu_{ij}) &= \mu + \tau_i + \beta_1\text{iba}_j + \beta_{1i}\text{iba}_j &\quad (2)\\ \end{align*} \] em que:
- \(y_{ijk}\) é o número médio de raízes observado na unidade experimental \(ijk\).
- \(n_{ijk}\) é o total de estacas na unidade experimental \(ijk\) a ser usado como peso.
- \(log(.)\) é a função de ligação log.
- \(\mu_{ij}\) é o número médio de raízes para a condição experimental \(ij\).
- \(\phi\) é o parâmetro de dispersão.
- \(\mu\) é uma constante comum à todas as celas experimentais.
- \(\tau_i\) (\(i = 1, 2\)) acomoda o efeito principal de sexo.
- \(\lambda_j\) (\(j = 1, \ldots, 5\)) acomoda o efeito principal de iba, considerado como fator qualitativo no modelo com preditor (1).
- \(\theta_{ij}\) acomoda o efeito de interação entre sexo e iba.
- \(\beta_1\) acomoda o efeito linear de iba, considerado como fator quantitativo no modelo com preditor (2).
- \(\beta_{1i}\) acomoda o efeito de interação entre sexo e iba com termos linear por nível de sexo.
O modelo foi ajustado aos dados por mínimos quadrados ponderados interativo.
ATENÇÃO
Deve-se decidir se a análise será conduzida com todas as estacas ou apenas com as estacas viáveis ou as que apresentaram enraizamento. São viáveis aquelas não consideradas mortas.
Na primeira opção não é feita distinção entre processos geradores de 0, ou seja, 0 raízes de uma estaca morta e 0 raízes de uma estaca viável são considerados produtos de um mesmo processo biológico. Ao considerar, por outro lado, a contagem de raízes apenas das estacas viáveis, o 0 indica que não houve enraizamento para uma estaca em condições de enraizar. Com isso, tem-se uma interpretação não ambígua do processo gerador de 0 e a interpretação se torna clara.
Para a análise que irá seguir com as variáveis resposta número de raízes e comprimento de raízes, será considerado apenas o subconjunto de estacas viáveis para avaliar o número de raízes e o comprimento de raízes.
# Filtra e agrega considerando as estacas não mortas.
tbv <- tb %>%
filter(morta == 0) %>%
mutate(n = 1) %>%
group_by(sexo, iba, rept) %>%
summarise_at(vars(numero, comprimento, n), sum)
# Número de estacas por cela experimental.
tbv %>%
xtabs(n ~ sexo + iba, .)## iba
## sexo 0 1.5 3 4.5 6
## Fêmea 6 13 8 7 7
## Macho 13 9 3 9 11# Ajuste do modelo.
m0 <- glm(numero/n ~ sexo * factor(iba),
weights = n,
data = tbv,
family = quasipoisson)
# m0 <- glm(numero ~ sexo * factor(iba),
# offset = log(n),
# data = tbv,
# family = quasipoisson)
# ATTENTION: para modelos baseados no Possion com log-link, o uso do
# `offset` ou `weights` resulta no mesmo modelo.
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de deviance.
anova(m0, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: numero/n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 29 224.799
## sexo 1 20.893 28 203.906 5.0511 0.036052 *
## factor(iba) 4 98.650 24 105.256 5.9625 0.002509 **
## sexo:factor(iba) 4 31.684 20 73.572 1.9150 0.147250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(m0, test = "F", scope = . ~ .)## Single term deletions
##
## Model:
## numero/n ~ sexo * factor(iba)
## Df Deviance F value Pr(>F)
## <none> 73.572
## sexo 1 73.572 0.0000 1.0000
## factor(iba) 4 100.385 1.8222 0.1642
## sexo:factor(iba) 4 105.256 2.1532 0.1116# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0, test = "F")## Analysis of Deviance Table
##
## Model 1: numero/n ~ iba
## Model 2: numero/n ~ sexo * factor(iba)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 28 124.613
## 2 20 73.572 8 51.041 1.5425 0.205# Estimativa dos parâmetros.
summary(m1)##
## Call:
## glm(formula = numero/n ~ iba, family = quasipoisson, data = tbv,
## weights = n)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.3231 -1.5696 -0.7552 -0.2661 5.1898
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.801 1.214 -2.307 0.02869 *
## iba 0.644 0.226 2.850 0.00811 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 7.330216)
##
## Null deviance: 224.80 on 29 degrees of freedom
## Residual deviance: 124.61 on 28 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba)),
type = "response")## iba rate SE df asymp.LCL asymp.UCL
## 0.0 0.0607 0.0738 Inf 0.00562 0.656
## 1.5 0.1596 0.1421 Inf 0.02787 0.914
## 3.0 0.4193 0.2447 Inf 0.13361 1.316
## 4.5 1.1016 0.3784 Inf 0.56189 2.160
## 6.0 2.8943 1.0243 Inf 1.44639 5.792
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale# As equações ajustadas que gerou os valores acima.
pars <- setNames(as.list(coef(m1)), c("b0", "b1"))
substitute(exp(b0 + b1 * iba), pars)## exp(-2.80108568149161 + 0.643969387343502 * iba)# Malha para a predição e construção do gráfico.
grid <- with(tb,
crossing(#sexo = levels(sexo),
iba = seq(min(iba),
max(iba),
length.out = 71)))
grid$numero <- predict(m1, newdata = grid, type = "response")
# Valores observados e ajustados.
leg <- c(0.05, 0.95)
ggplot(data = tbv,
mapping = aes(x = iba,
y = numero/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = numero)) +
labs(x = "Concentração de IBA (milhar)",
y = "Número médio de raízes por estaca",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------Comprimento médio de raízes
O modelo associado a essa reposta é \[ \begin{align*} t(y_{ijk}) &\sim \text{Normal}(\mu_{ij}, \sigma^2)\\ \mu_{ij} &= \mu + \tau_i + \lambda_j + \theta_{ij} &\quad (1)\\ \mu_{ij} &= \mu + \tau_i + \beta_1\text{iba}_j + \beta_{1i}\text{iba}_j &\quad (2)\\ \end{align*} \] em que:
- \(y_{ijk}\) é o comprimento médio de raízes observado na unidade experimental \(ijk\).
- \(t(.)\) é transformação aplicada na variável resposta visando mitigar o afastamento dos pressupostos. No caso, foi usada a transformação \(y_t = t(y) = \log(y + \alpha)\), sendo \(\alpha = 0.01\) o valor utilizado baseado no perfil da função de log-verossimilhança em \(\alpha\).
- \(\mu_{ij}\) é o comprimento médio de raízes para a condição experimental \(ij\).
- \(\sigma^2\) é a variância do erro experimental.
- \(\mu\) é uma constante comum à todas as celas experimentais.
- \(\tau_i\) (\(i = 1, 2\)) acomoda o efeito principal de sexo.
- \(\lambda_j\) (\(j = 1, \ldots, 5\)) acomoda o efeito principal de iba, considerado como fator qualitativo no modelo com preditor (1).
- \(\theta_{ij}\) acomoda o efeito de interação entre sexo e iba.
- \(\beta_1\) acomoda o efeito linear de iba, considerado como fator quantitativo no modelo com preditor (2).
- \(\beta_{1i}\) acomoda o efeito de interação entre sexo e iba com termos linear por nível de sexo.
O modelo foi ajustado aos dados por mínimos quadrados.
# Ajuste do modelo.
m0 <- lm(comprimento/n ~ sexo * factor(iba),
data = tbv)
MASS::logtrans(m0, alpha = seq(0.1, 3, length.out = 31))
abline(v = 0.01, col = "red")
# Usa resposta transformada.
m0 <- update(m0, log(comprimento/n + 0.01) ~ .)
# Inspeção dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variâcia para variável transformada.
anova(m0)## Analysis of Variance Table
##
## Response: log(comprimento/n + 0.01)
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 1 0.110 0.1104 0.0161 0.90027
## factor(iba) 4 79.659 19.9146 2.9045 0.04794 *
## sexo:factor(iba) 4 39.244 9.8109 1.4309 0.26033
## Residuals 20 137.132 6.8566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo final baseado no abandono de termos não relevantes e com termo
# linear para o efeito de iba.
m1 <- update(m0, formula = . ~ iba)
# Teste F entre modelos encaixados.
anova(m1, m0)## Analysis of Variance Table
##
## Model 1: log(comprimento/n + 0.01) ~ iba
## Model 2: log(comprimento/n + 0.01) ~ sexo + factor(iba) + sexo:factor(iba)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 182.30
## 2 20 137.13 8 45.172 0.8235 0.5917# Estimativa dos parâmetros.
summary(m1)##
## Call:
## lm(formula = log(comprimento/n + 0.01) ~ iba, data = tbv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5937 -1.2849 -0.1819 1.8021 5.7344
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.4232 0.8307 -5.325 1.14e-05 ***
## iba 0.7353 0.2183 3.368 0.00222 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.552 on 28 degrees of freedom
## Multiple R-squared: 0.2883, Adjusted R-squared: 0.2629
## F-statistic: 11.34 on 1 and 28 DF, p-value: 0.00222# Resultados com o modelo final com termo linear para efeito de iba.
emmeans(m1,
specs = ~iba,
at = list(iba = unique(tb$iba))) %>%
as.data.frame() %>%
mutate_at(c(2, 5, 6), ~exp(.) - 0.01)## iba emmean SE df lower.CL upper.CL
## 1 0.0 0.001995454 0.8306973 28 -0.0078121452 0.05576804
## 2 1.5 0.026142527 0.5889114 28 0.0008172253 0.11075946
## 3 3.0 0.098898113 0.4670131 28 0.0318369434 0.27345280
## 4 4.5 0.318112055 0.5512817 28 0.0960706909 1.00496012
## 5 6.0 0.978607775 0.7773337 28 0.1911401269 4.84902713# As equações ajustadas que gerou os valores acima.
pars <- setNames(as.list(coef(m1)), c("b0", "b1"))
substitute(exp(b0 + b1 * iba) - 0.01, pars)## exp(-4.42322754756973 + 0.735294989056067 * iba) - 0.01# Malha para a predição e construção do gráfico.
grid <- with(tb,
crossing(#sexo = levels(sexo),
iba = seq(min(iba),
max(iba),
length.out = 71)))
grid$log_comprimento <- predict(m1, newdata = grid)
grid$comprimento <- exp(grid$log_comprimento) - 0.01
# Valores observados e ajustados.
leg <- c(0.05, 0.95)
ggplot(data = tbv,
mapping = aes(x = iba,
# y = log(comprimento/n + 0.01),
y = comprimento/n,
color = sexo)) +
# geom_point() +
geom_jitter(width = 0.1, height = 0) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = iba, y = comprimento)) +
labs(x = "Concentração de IBA (milhar)",
y = "Comprimento médio de raízes",
color = "Tipo de estaca") +
theme(legend.position = leg,
legend.justification = round(leg))
#-----------------------------------------------------------------------