Priscila Santos
Walmes M. Zeviani
Análise preliminar.
Incluír o básico.
##----------------------------------------------------------------------------
## Session Definition
pkg <- c("lattice",
"latticeExtra",
"gridExtra",
"doBy",
"multcomp",
"plyr",
"MASS",
"wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
##-----------------------------------------------------------------------------
## Trellis plots settings.
lattice.options(default.theme=modifyList(standard.theme(color=FALSE),
list(strip.background=list(col="transparent"))))
## trellis.device(color=FALSE)
##----------------------------------------------------------------------------
## Read the data.
## list.files()
da <- read.table("molhamento_foliar.txt", header=TRUE,
sep="\t", stringsAsFactors=TRUE)
## str(da)
xt <- xtabs(~temp+dpm, data=da)
addmargins(xt)
## dpm
## temp 4 8 12 16 20 24 28 Sum
## 10 18 18 18 18 18 18 18 126
## 15 18 18 18 18 18 18 18 126
## 20 18 18 18 18 18 18 18 126
## 25 18 18 18 18 18 18 18 126
## 30 18 18 18 18 18 18 18 126
## Sum 90 90 90 90 90 90 90 630
dn <- dimnames(xt)
fl <- c(outer(dn$dpm, dn$temp, paste, sep="."))
da <- transform(da,
trat=interaction(dpm, temp))
da$trat <- factor(da$trat, levels=fl)
str(da)
## 'data.frame': 630 obs. of 5 variables:
## $ trat : Factor w/ 35 levels "4.10","8.10",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ rep : int 1 2 3 4 5 6 7 8 9 1 ...
## $ temp : int 10 10 10 10 10 10 10 10 10 10 ...
## $ dpm : int 4 4 4 4 4 4 4 4 4 8 ...
## $ lesao: num 0 0 0 0 0 0 0 0 0 0 ...
##-----------------------------------------------------------------------------
## Gráficos.
xyplot(lesao~dpm|temp, data=da)
## Criar a variável ocorrência = {0,1}.
da$oco <- sign(da$lesao)
## Gráfico de ocorrência ou não de lesão.
xyplot(oco~dpm|temp, data=da)
## Gráfico do tamanho da lesão quando houve ocorrência.
xyplot(lesao~dpm|temp, data=subset(da, oco==1))
##-----------------------------------------------------------------------------
## Modelo binomial para ocorrência.
m0 <- glm(oco~temp+I(temp^2)+dpm+I(dpm^2)+temp:dpm, data=da,
family=binomial)
par(mfrow=c(2,2)); plot(m0); layout(1)
summary(m0)
##
## Call:
## glm(formula = oco ~ temp + I(temp^2) + dpm + I(dpm^2) + temp:dpm,
## family = binomial, data = da)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.38881 -0.25189 -0.01194 -0.00074 2.36401
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -61.03995 12.35179 -4.942 7.74e-07 ***
## temp 3.31682 0.62644 5.295 1.19e-07 ***
## I(temp^2) -0.06560 0.01043 -6.291 3.15e-10 ***
## dpm 1.82894 0.64304 2.844 0.00445 **
## I(dpm^2) -0.02909 0.01134 -2.564 0.01034 *
## temp:dpm -0.01257 0.01378 -0.913 0.36150
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 405.19 on 629 degrees of freedom
## Residual deviance: 200.63 on 624 degrees of freedom
## AIC: 212.63
##
## Number of Fisher Scoring iterations: 10
##-----------------------------------------------------------------------------
## Modelo GAM.
require(mgcv)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-3. For overview type 'help("mgcv-package")'.
m0 <- gam(oco~s(temp, k=3)+s(dpm, k=3), data=da,
family=quasibinomial)
summary(m0)
##
## Family: quasibinomial
## Link function: logit
##
## Formula:
## oco ~ s(temp, k = 3) + s(dpm, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.7165 0.8974 -9.713 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(temp) 1.991 2.000 55.57 < 2e-16 ***
## s(dpm) 1.953 1.998 37.35 4.16e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.405 Deviance explained = 51.2%
## GCV = 0.31867 Scale est. = 0.32436 n = 630
##-----------------------------------------------------------------------------
## Splines.
## require(splines)
## m0 <- glm(oco~bs(temp, df=4)+bs(dpm, df=4), data=da,
## family=binomial)
## m0 <- glm(oco~bs(temp)+bs(dpm)+bs(I(temp*dpm)), data=da,
## family=quasibinomial)
## m1 <- glm(oco~bs(temp)+bs(dpm), data=da,
## family=quasibinomial)
## anova(m0, m1, test="F")
## summary(m0)
##-----------------------------------------------------------------------------
## Gráfico da probabilidade de ocorrência.
pred <- with(da,
expand.grid(
dpm=seq(min(dpm), max(dpm), l=50),
temp=seq(min(temp), max(temp), l=50)
))
pred$p <- predict(m0, newdata=pred, type="response")
## levelplot(p~temp+dpm, data=pred)
wireframe(p~dpm+temp, data=pred,
zlim=c(0,1),
scales=list(arrows=FALSE),
panel.3d.wireframe=panel.3d.contour,
levels=seq(0,1,by=0.05),
type=c("on", "top"),
drape=TRUE)
##----------------------------------------------------------------------------
## Versions.
sessionInfo()
## R version 3.1.2 (2014-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines grid stats graphics grDevices utils datasets
## [9] base
##
## other attached packages:
## [1] mgcv_1.8-3 nlme_3.1-118 wzRfun_0.5 MASS_7.3-35
## [5] plyr_1.8.1 multcomp_1.3-8 TH.data_1.0-4 mvtnorm_1.0-1
## [9] doBy_4.5-12 survival_2.37-7 gridExtra_0.9.1 latticeExtra_0.6-26
## [13] RColorBrewer_1.0-5 lattice_0.20-29 rmarkdown_0.3.10 knitr_1.8
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.4 Matrix_1.1-4
## [6] Rcpp_0.11.3 sandwich_2.3-2 stringr_0.6.2 tools_3.1.2 yaml_2.1.13
## [11] zoo_1.7-11
##-----------------------------------------------------------------------------
## Citation.
citation()
##
## To cite R in publications use:
##
## R Core Team (2014). R: A language and environment for statistical computing. R
## Foundation for Statistical Computing, Vienna, Austria. URL
## http://www.R-project.org/.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2014},
## url = {http://www.R-project.org/},
## }
##
## We have invested a lot of time and effort in creating R, please cite it when
## using it for data analysis. See also 'citation("pkgname")' for citing R
## packages.
##-----------------------------------------------------------------------------
## Last modification.
Sys.time()
## [1] "2014-12-25 01:07:46 BRST"