Dispersão de esporos com GLM Poisson
Louise Larissa May De Mio & Walmes Marques Zeviani
Preparação do ambiente R
rm(list = objects())
# Pacotes.
library(gdata)
library(latticeExtra)
library(wzRfun)
# names(trellis.par.get()$box.3d)
# names(trellis.par.get()$box.dot)
mycol <- c("#E41A1C", "#377EB8", "#4DAF4A",
"#984EA3", "#FF7F00", "#FFFF33")
mypch <- c(1, 5)
mycex <- 0.8 * c(1, 0.9)
ps <- list(box.rectangle = list(col = 1, fill = c("gray70")),
box.umbrella = list(col = 1, lty = 1),
box.dot = list(pch = 1),
dot.symbol = list(col = 1, pch = 19),
dot.line = list(col = "gray50", lty = 3),
plot.symbol = list(col = 1, cex = mycex[1]),
plot.line = list(col = 1),
plot.polygon = list(col = "gray95"),
# superpose.line = list(col = mycol, lty = 1),
superpose.line = list(col = 1, lty = 1:2),
# superpose.symbol = list(col = mycol, pch = 1, fill = mycol),
superpose.symbol = list(col = 1, pch = mypch, cex = mycex),
# superpose.region = list(col = mycol, pch = 1),
# superpose.polygon = list(col = mycol),
strip.background = list(col = c("gray80", "gray50")),
axis.text = list(cex = 0.8))
trellis.par.set(ps)
lattice.options(default.args = list(as.table = TRUE))
# Para ter nomes dos meses em inglês.
Sys.setlocale("LC_TIME", "en_US.UTF-8")
# Informações da sessão.
# sessionInfo()
devtools::session_info()Safra de 2011
Importação e preparo
#-----------------------------------------------------------------------
# Lê e organiza a tabela de dados.
f <- "Planilha total dados 2011 e 2013.xlsx"
da1 <- read.xls(f,
sheet = 1,
header = TRUE,
stringsAsFactors = FALSE)
# Mantém só colunas com informações úteis. Renomeia variáveis.
da1 <- da1[, 1:4]
names(da1) <- c("ts", "dia", "dis", "numesp")
# Tranforma para data.
da1$ts <- as.Date(da1$ts)
# Estrutura.
str(da1)## 'data.frame': 40 obs. of 4 variables:
## $ ts : Date, format: "2011-12-12" "2011-12-12" ...
## $ dia : int 0 0 0 0 0 4 4 4 4 4 ...
## $ dis : int 0 92 104 118 130 0 92 104 118 130 ...
## $ numesp: int 113 31 25 27 18 107 24 11 10 4 ...Análise exploratória
#-----------------------------------------------------------------------
xyplot(numesp ~ ts | factor(dis),
data = da1,
type = "o",
xlab = "Date",
ylab = "Number of spores in the trap",
main = list("Season 2011/12", font = 1),
as.table = TRUE)
xyplot(numesp ~ dis | factor(ts),
data = da1,
type = "o",
xlab = "Trap distance from the orchard border (m)",
ylab = "Number of spores in the trap",
main = list("Season 2011/12", font = 1),
as.table = TRUE)
Ajuste do modelo linear generalizado
Devido a resposta ser uma variável de contagem, assumiu-se um modelo baseado em quasi-verossimilhança que tem os dois primeiros momentos iguais ao da distribuição de Poisson e estima um parâmetro extra de dispersão, o que é necessário para acomodar situações de super/sub dispersão. A função de ligação foi a logarÃtmica (assim a função inversa é a exponencial) e a função de variância foi a identidade, ou seja, a relação média variância é dada por \(\textrm{Var}(Y|x) = \phi \textrm{E}(Y|x)\). O modelo para explicar a média do número de esporos considerou, inicialmente, os produtos de termos lineares e quadráticos das variáveis distância da borda do pomar (dis) e perÃodo após a primeira coleta (dia). Os termos não significativos, de acordo com a estatÃstica \(z\) do teste de Wald, foram abandonados para obtenção de um modelo mais parcimonioso.
#-----------------------------------------------------------------------
# Ajuste com GLM de famÃlia quasi-Poisson.
# Declarando um modelo para a contagem usando quasi-Poisson.
g0 <- glm(numesp ~ poly(dia, 2) * poly(dis, 2),
data = da1,
family = quasipoisson)
# Reajuste com polinômio cru (trata-se ainda do mesmo modelo).
g0 <- update(g0, ~(dia + I(dia^2)) * (dis + I(dis^2)))
# Análise de diagnóstico.
par(mfrow = c(2, 2))
plot(g0)
layout(1)
# Teste para os termos de efeito fixo do modelo via estatÃstica F da
# diferença de deviances.
anova(g0, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: numesp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 39 1185.03
## dia 1 581.19 38 603.85 375.6570 < 2.2e-16 ***
## I(dia^2) 1 20.16 37 583.69 13.0280 0.001067 **
## dis 1 524.73 36 58.96 339.1645 < 2.2e-16 ***
## I(dis^2) 1 1.47 35 57.50 0.9481 0.337735
## dia:dis 1 6.77 34 50.73 4.3759 0.044734 *
## dia:I(dis^2) 1 0.36 33 50.36 0.2359 0.630612
## I(dia^2):dis 1 1.22 32 49.14 0.7882 0.381496
## I(dia^2):I(dis^2) 1 0.70 31 48.45 0.4499 0.507370
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A análise dos resÃduos não indica fuga nos pressupostos do modelo, sendo portanto o modelo considerado adequadamente ajustado aos dados.
De acordo com o quadro de deviance acima, pela estatÃstica F, não se tem evidência contra a hipótese de não haver interação entre os fatores. O modelo reduzido, porém, apontou interação entre os termos lineares ao nÃvel de 5%.
# Resumo do ajuste do modelo maximal.
summary(g0)##
## Call:
## glm(formula = numesp ~ dia + I(dia^2) + dis + I(dis^2) + dia:dis +
## dia:I(dis^2) + I(dia^2):dis + I(dia^2):I(dis^2), family = quasipoisson,
## data = da1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2373 -0.8710 -0.1429 0.2724 2.5063
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.731e+00 9.870e-02 47.935 < 2e-16 ***
## dia -1.147e-02 1.690e-02 -0.679 0.50214
## I(dia^2) -1.672e-03 5.157e-04 -3.241 0.00284 **
## dis -1.343e-02 8.663e-03 -1.550 0.13118
## I(dis^2) -1.572e-05 7.616e-05 -0.206 0.83787
## dia:dis 8.954e-04 1.704e-03 0.525 0.60299
## dia:I(dis^2) -1.200e-05 1.517e-05 -0.791 0.43485
## I(dia^2):dis -2.768e-05 5.315e-05 -0.521 0.60626
## I(dia^2):I(dis^2) 3.184e-07 4.698e-07 0.678 0.50294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.547118)
##
## Null deviance: 1185.032 on 39 degrees of freedom
## Residual deviance: 48.446 on 31 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6# Modelo reduzido obtido com o abandono de termos irrelevantes.
g1 <- update(g0, ~dia + dis + I(dia^2))
# Análise de diagnóstico.
par(mfrow = c(2, 2))
plot(g1)
layout(1)
# Teste para os termos de efeito fixo abadonados do modelo via
# estatÃstica F da diferença de deviances.
anova(g1, g0, test = "F")## Analysis of Deviance Table
##
## Model 1: numesp ~ dia + dis + I(dia^2)
## Model 2: numesp ~ dia + I(dia^2) + dis + I(dis^2) + dia:dis + dia:I(dis^2) +
## I(dia^2):dis + I(dia^2):I(dis^2)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 36 58.963
## 2 31 48.446 5 10.517 1.3596 0.2663# Resumo do ajuste.
summary(g1)##
## Call:
## glm(formula = numesp ~ dia + dis + I(dia^2), family = quasipoisson,
## data = da1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5092 -1.0794 -0.4905 0.5374 3.1572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8285351 0.0841405 57.387 < 2e-16 ***
## dia -0.0272389 0.0134481 -2.025 0.05028 .
## dis -0.0174778 0.0009651 -18.110 < 2e-16 ***
## I(dia^2) -0.0013982 0.0004146 -3.372 0.00179 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.523216)
##
## Null deviance: 1185.032 on 39 degrees of freedom
## Residual deviance: 58.963 on 36 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Equação estimada para previsão.
beta <- coef(g1)
eqn <- data.frame(sign = ifelse(beta >= 0, "+", "-"),
abs = abs(beta))
eqn$sign[1] <- ifelse(eqn$sign[1] == "+", "", "-")
eqn## sign abs
## (Intercept) 4.828535098
## dia - 0.027238891
## dis - 0.017477768
## I(dia^2) - 0.001398151O modelo reduzido, considerando teste F para a diferença de deviance entre modelos encaixados, não difere em qualidade do modelo maximal. Dessa forma, esse modelo será considerado para a interpretação dos efeitos e previsão de quantidade de esporos.
Ess modelo final tem a seguinte equação ajustada para o número médio de esporos (\(\hat{y}\)) \[ \hat{y} = \exp\{ 4.829 - 0.0272 \times \texttt{dia} - 0.0014 \times \texttt{dia}^2 - 0.0175 \times \texttt{dis} \} \]
Predição ao longo das datas fixando as distâncias
#-----------------------------------------------------------------------
# Predição.
# Malha para predição em cada distância como função do tempo.
grid <- with(da1,
expand.grid(dis = sort(unique(dis)),
dia = seq(min(dia),
max(dia),
length.out = 30)))
# Predição na escala da resposta.
grid$z <- predict(g0, newdata = grid, type = "response")
# Visualizando o ajuste.
xy1_diaINdis <-
xyplot(numesp ~ dia | factor(dis),
data = da1,
xlab = "Time after first assessement (days)",
ylab = "Number of spores in the trap",
main = list("Season 2011/12", font = 1),
as.table = TRUE) +
as.layer(xyplot(z ~ dia | dis,
data = grid,
type = "l"))
xy1_diaINdis
Predição ao longo da distância fixando as datas
#-----------------------------------------------------------------------
# Predição.
# Malha para predição em cada distância como função do tempo.
grid <- with(da1,
expand.grid(dis = seq(min(dis),
max(dis),
length.out = 30),
dia = sort(unique(dia))[1:5]))
# Predição na escala da resposta.
grid$z <- predict(g0, newdata = grid, type = "response")
grid <- merge(grid, unique(da1[, c("dia", "ts")]),
by = "dia", all.x = TRUE)
grid$ts2 <- strftime(grid$ts, format = "%b %d, %Y")
da1$ts2 <- strftime(da1$ts, format = "%b %d, %Y")
# Visualizando o ajuste.
xy1_disINdia <-
xyplot(numesp ~ dis | factor(ts2),
data = subset(da1, dia <= max(grid$dia)),
xlab = "Distance (m)",
ylab = "Number of spores in the trap",
main = list("Season 2011/12", font = 1),
as.table = TRUE) +
as.layer(xyplot(z ~ dis | factor(ts2),
data = grid,
type = "l"))
xy1_disINdia
Predição usando distância e data
#-----------------------------------------------------------------------
# Predição em um grid para gerar tabela.
grid <- with(da1,
expand.grid(dis = c(0, 10, 25, 50, 100),
dia = seq(0, 60, by = 15)))
grid$z <- predict(g0, newdata = grid, type = "response")
tb <- reshape2::dcast(grid, formula = dis ~ dia, value.var = "z")
cap <-
"Número de esporos como função da distância (m) e do perÃodo após a primeira coleta (dias)."
knitr::kable(tb, caption = cap)| dis | 0 | 15 | 30 | 45 | 60 |
|---|---|---|---|---|---|
| 0 | 113.45746 | 65.57520 | 17.862473 | 2.2931839 | 0.1387495 |
| 10 | 99.04208 | 60.85600 | 15.783938 | 1.7280473 | 0.0798591 |
| 25 | 80.30460 | 51.93232 | 12.679416 | 1.1687593 | 0.0406738 |
| 50 | 55.73353 | 35.21799 | 8.050711 | 0.6657721 | 0.0199177 |
| 100 | 25.30901 | 10.16961 | 2.322929 | 0.3016267 | 0.0222641 |
#-----------------------------------------------------------------------
# Predição para gerar a superficÃe de resposta.
grid_lv <- with(da1,
expand.grid(dis = seq(0, 200, length.out = 151),
dia = seq(0, 50, length.out = 151)))
grid_lv$z <- predict(g0, newdata = grid_lv, type = "response")
at <- c(Inf, 90, 70, 50, 30, 10, 1, 0)
lv1 <- levelplot(z ~ dia + dis,
data = grid_lv,
xlab = "Time after the first assessement (days)",
ylab = "Distance from orchard border (m)",
main = list("Season 2011/12", font = 1),
# col.regions = gray.colors,
col.regions = colorRampPalette(
colors = c("white", "gray60")), #(length(at) + 1),
at = at,
contour = FALSE) +
as.layer(
contourplot(z ~ dia + dis,
data = grid_lv,
# label.style = "flat",
at = at))
lv1
grid_wf <- with(da1,
expand.grid(dis = seq(0, 200, length.out = 31),
dia = seq(0, 50, length.out = 31)))
grid_wf$z <- predict(g0, newdata = grid_wf, type = "response")
# Visualizando o ajuste.
w1 <- wireframe(z ~ dia + dis,
data = grid_wf,
drape = TRUE,
xlab = list("Time after the\nfirst assessement (days)",
rot = 32),
ylab = list("Distance from\norchard border (m)",
rot = -32),
zlab = list("Number of spores", rot = 90),
main = list("Season 2011/12", font = 1),
alpha.regions = 0.25,
col.regions =
# colorRampPalette(colors = c("white", "black"))(101),
colorRampPalette(colors = c("orange", "blue"))(101),
screen = list(x = -110, z = -20, y = -140),
scpos = list(x = 9, y = 5, z = 2),
type = "on",
panel.3d.wireframe = panel.3d.contour,
scales = list(arrows = FALSE)) +
layer(data = da1,
panel.cloud(z = numesp,
x = dia,
y = dis,
type = "p",
# pch = 19,
# col = "red",
...))
print(w1)
# rp.wire(w)Safra 2012/2013
Importação e preparo
#-----------------------------------------------------------------------
# Leitura dos dados.
# Linha de armadilhas I.
da2 <- read.xls(f, sheet = 2, header = TRUE, stringsAsFactors = FALSE)
da2 <- da2[, 1:4]
names(da2) <- c("ts", "dia", "dis", "numesp")
da2$line <- "I"
# Linha de armadilhas II.
da3 <- read.xls(f, sheet = 3, header = TRUE, stringsAsFactors = FALSE)
da3 <- da3[, 1:4]
names(da3) <- c("ts", "dia", "dis", "numesp")
da3$line <- "II"
# Justa as linhas.
da <- rbind(da2, da3)
da$line <- factor(da$line)
da$ts <- as.Date(da$ts, format = "%d/%m/%y")
str(da)## 'data.frame': 480 obs. of 5 variables:
## $ ts : Date, format: "2012-12-18" "2012-12-18" ...
## $ dia : int 0 0 0 0 0 0 0 0 7 7 ...
## $ dis : int 0 4 8 16 32 64 128 276 0 4 ...
## $ numesp: int NA 3 5 11 2 3 2 1 10 9 ...
## $ line : Factor w/ 2 levels "I","II": 1 1 1 1 1 1 1 1 1 1 ...# A origem é na borda do pomar, aos 16 metros.
da$dis <- da$dis - 16
da <- subset(da, dis >= 0)
str(da)## 'data.frame': 288 obs. of 5 variables:
## $ ts : Date, format: "2012-12-18" "2012-12-18" ...
## $ dia : int 0 0 0 0 0 7 7 7 7 7 ...
## $ dis : num 0 16 48 112 260 0 16 48 112 260 ...
## $ numesp: int 11 2 3 2 1 28 3 2 2 1 ...
## $ line : Factor w/ 2 levels "I","II": 1 1 1 1 1 1 1 1 1 1 ...Análise exploratória
#-----------------------------------------------------------------------
# Visualização.
# # A partir de qual data os resgistros são todos nulos?
# head(aggregate(cbind(y = numesp > 0) ~ dia,
# data = da,
# FUN = sum),
# n = 13)
# Descartado a porção nula trazida com o aumento da data.
# db <- subset(da, dia <= 80)
db <- subset(da, dia <= 50)
xyplot(numesp ~ ts | factor(dis),
groups = line,
data = db,
subset = is.finite(numesp),
type = "o",
xlab = "Date",
ylab = "Number of spores in the trap",
main = list("Season 2012/13", font = 1),
as.table = TRUE)
xyplot(numesp ~ dis | factor(ts),
groups = line,
data = db,
subset = is.finite(numesp),
type = "o",
xlab = "Trap distance from the orchard border (m)",
ylab = "Number of spores in the trap",
main = list("Season 2012/13", font = 1),
as.table = TRUE)
Ajuste do modelo linear generalizado
O modelo foi especificado de forma equivalente a safra anterior, porém houve a inclusão do efeito de linha de coleta de esporos. As linhas funcionam como blocos e são incluÃdas no modelo para acomodar variações locais que podem mudar a intensidade de esporos, como por exemplo, pela influência da direção do vento predominante.
#-----------------------------------------------------------------------
# Declarando um modelo para a contagem usando quasi-Poisson.
g0 <- glm(numesp ~ line * poly(dia, 2) * poly(dis, 2),
# family = quasi(link = "log", variance = "mu"),
# family = quasi(link = "log", variance = "mu^2"),
family = quasipoisson,
data = db)
# Reajuste com polinômio cru.
g0 <- update(g0, ~line * (dia + I(dia^2)) * (dis + I(dis^2)))
# ResÃduos.
par(mfrow = c(2, 2))
plot(g0)
layout(1)
# Resumo do ajuste.
anova(g0, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: numesp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 66 302.95
## line 1 10.995 65 291.96 3.7793 0.0576441 .
## dia 1 86.279 64 205.68 29.6572 1.65e-06 ***
## I(dia^2) 1 19.000 63 186.68 6.5310 0.0137545 *
## dis 1 42.827 62 143.85 14.7211 0.0003571 ***
## I(dis^2) 1 9.448 61 134.40 3.2476 0.0776799 .
## line:dia 1 0.000 60 134.40 0.0000 0.9955563
## line:I(dia^2) 1 2.618 59 131.78 0.8998 0.3474857
## line:dis 1 0.342 58 131.44 0.1177 0.7330101
## line:I(dis^2) 1 0.054 57 131.39 0.0185 0.8924889
## dia:dis 1 2.983 56 128.41 1.0252 0.3162560
## dia:I(dis^2) 1 4.344 55 124.06 1.4932 0.2275667
## I(dia^2):dis 1 1.495 54 122.57 0.5139 0.4768473
## I(dia^2):I(dis^2) 1 0.437 53 122.13 0.1502 0.7000111
## line:dia:dis 1 3.216 52 118.91 1.1055 0.2982205
## line:dia:I(dis^2) 1 1.376 51 117.54 0.4731 0.4947919
## line:I(dia^2):dis 1 0.436 50 117.10 0.1500 0.7002036
## line:I(dia^2):I(dis^2) 1 0.000 49 117.10 0.0001 0.9918406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo reduzido.
g1 <- update(g0, ~line + dia + dis)
# Teste para os termos abandonados.
anova(g1, g0, test = "F")## Analysis of Deviance Table
##
## Model 1: numesp ~ line + dia + dis
## Model 2: numesp ~ line + dia + I(dia^2) + dis + I(dis^2) + line:dia +
## line:I(dia^2) + line:dis + line:I(dis^2) + dia:dis + dia:I(dis^2) +
## I(dia^2):dis + I(dia^2):I(dis^2) + line:dia:dis + line:dia:I(dis^2) +
## line:I(dia^2):dis + line:I(dia^2):I(dis^2)
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 63 163.12
## 2 49 117.10 14 46.024 1.13 0.3574# ResÃduos.
par(mfrow = c(2, 2))
plot(g1)
layout(1)
# Resumo do ajuste.
summary(g1)##
## Call:
## glm(formula = numesp ~ line + dia + dis, family = quasipoisson,
## data = db)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1385 -1.2767 -0.6671 0.3074 5.4808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.464412 0.251475 9.800 2.74e-14 ***
## lineII -0.860631 0.335914 -2.562 0.01281 *
## dia -0.054583 0.011676 -4.675 1.60e-05 ***
## dis -0.009001 0.003014 -2.987 0.00401 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.184985)
##
## Null deviance: 302.95 on 66 degrees of freedom
## Residual deviance: 163.13 on 63 degrees of freedom
## (5 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 6beta <- coef(g1)
beta[1] <- beta[1] + 0.5 * beta[2]
beta <- beta[-2]
eqn <- data.frame(sign = ifelse(beta >= 0, "+", "-"),
abs = abs(beta))
eqn$sign[1] <- ifelse(eqn$sign[1] == "+", "", "-")
eqn## sign abs
## (Intercept) 2.034096195
## dia - 0.054583409
## dis - 0.009001397Para esta safra, não houve interação de distância com tempo. Portanto, o modelo final é um modelo aditivo nos termos linha, distância e tempo. Considerando o efeito médio de linha, a seguinte equação representa o número médio de esporos (\(\hat{y}\)) \[ \hat{y} = \exp\{ 2.034 - 0.0546 \times \texttt{dia} - 0.009 \times \texttt{dis} \} \]
Predição ao longo das datas fixando as distâncias
#-----------------------------------------------------------------------
# Predição.
# Predição média para as linhas.
grid <- with(db,
expand.grid(line = factor(levels(line)[1], levels(line)),
dis = sort(unique(dis)),
dia = seq(min(dia),
max(dia),
length.out = 30)))
X <- model.matrix(formula(g1)[-2], data = grid)
X[, "lineII"] <- 0.5
nu <- X %*% coef(g1)
grid$z <- g1$family$linkinv(nu)
grid$line <- NULL
# Visualizando o ajuste.
xy2_diaINdis <-
xyplot(numesp ~ dia | factor(dis),
groups = line,
auto.key = list(corner = c(0.975, 0.80),
title = "Line", cex.title = 1),
data = db,
xlab = "Time after first assessement (days)",
ylab = "Number of spores in the trap",
main = list("Season 2012/13", font = 1),
as.table = TRUE) +
as.layer(xyplot(z ~ dia | dis,
data = grid,
type = "l"))
xy2_diaINdis
Predição ao longo da distância fixando as datas
#-----------------------------------------------------------------------
# Predição.
# Predição com efeito médio das linhas.
grid <- with(db,
expand.grid(line = factor(unique(line)[1], levels(line)),
dis = seq(0, 250, length.out = 51),
dia = sort(unique(dia))[1:5],
KEEP.OUT.ATTRS = FALSE))
X <- model.matrix(formula(g1)[-2], data = grid)
X[, "lineII"] <- 0.5
nu <- X %*% coef(g1)
grid$z <- g1$family$linkinv(nu)
grid$line <- NULL
grid <- merge(grid, unique(db[, c("dia", "ts")]),
by = "dia", all.x = TRUE)
grid$ts2 <- strftime(grid$ts, format = "%b %d, %Y")
db$ts2 <- strftime(db$ts, format = "%b %d, %Y")
# Visualizando o ajuste.
xy2_disINdia <-
xyplot(numesp ~ dis | factor(ts2),
group = line,
auto.key = list(corner = c(0.975, 0.8),
title = "Line", cex.title = 1),
data = subset(db, dia <= max(grid$dia)),
xlab = "Distance (m)",
ylab = "Number of spores in the trap",
main = list("Season 2012/13", font = 1),
as.table = TRUE) +
as.layer(xyplot(z ~ dis | factor(ts2),
data = grid,
type = "l"))
xy2_disINdia
Predição usando distância e data
# Predição com efeito médio das linhas.
grid <- with(db,
expand.grid(line = factor(levels(line)[1],
levels = levels(line)),
dis = c(0, 10, 25, 50, 100, 300),
dia = seq(0, 60, by = 15),
KEEP.OUT.ATTRS = FALSE))
X <- model.matrix(formula(g1)[-2], data = grid)
colnames(X)## [1] "(Intercept)" "lineII" "dia" "dis"X[, "lineII"] <- 0.5
nu <- X %*% coef(g1)
grid$z <- g1$family$linkinv(nu)
grid$line <- NULL
tb <- reshape2::dcast(grid, formula = dis ~ dia, value.var = "z")
cap <-
"Número de esporos como função da distância (m) e do perÃodo após a primeira coleta (dias)."
knitr::kable(tb, caption = cap)| dis | 0 | 15 | 30 | 45 | 60 |
|---|---|---|---|---|---|
| 0 | 7.6453391 | 3.3714572 | 1.4867521 | 0.6556310 | 0.2891215 |
| 10 | 6.9872162 | 3.0812369 | 1.3587701 | 0.5991932 | 0.2642335 |
| 25 | 6.1047141 | 2.6920693 | 1.1871542 | 0.5235137 | 0.2308601 |
| 50 | 4.8745430 | 2.1495859 | 0.9479288 | 0.4180196 | 0.1843391 |
| 100 | 3.1079288 | 1.3705408 | 0.6043839 | 0.2665224 | 0.1175316 |
| 300 | 0.5135937 | 0.2264856 | 0.0998761 | 0.0440436 | 0.0194224 |
# Predição com efeito médio das linhas.
grid_lv <- with(db,
expand.grid(line = levels(line),
dis = seq(0, 200, length.out = 151),
dia = seq(0, 50, length.out = 151),
KEEP.OUT.ATTRS = FALSE))
# grid$z <- predict(g1, newdata = grid, type = "response")
X <- model.matrix(formula(g1)[-2], data = grid_lv)
colnames(X)## [1] "(Intercept)" "lineII" "dia" "dis"X[, "lineII"] <- 0.5
nu <- X %*% coef(g1)
grid_lv$z <- g1$family$linkinv(nu)
grid_lv$line <- NULL
# range(grid$z)
# at <- seq(0, 10, by = 1)
at <- c(Inf, 7, 5, 2, 3, 1, 0)
lv2 <- levelplot(z ~ dia + dis,
data = grid_lv,
xlab = "Time after the first assessement (days)",
ylab = "Distance from orchard border (m)",
main = list("Season 2012/13", font = 1),
col.regions = colorRampPalette(
colors = c("white", "gray60")), #(length(at) + 1),
at = at,
contour = FALSE) +
as.layer(
contourplot(z ~ dia + dis,
data = grid_lv,
label.style = "flat",
at = at))
print(lv2)