Walmes Zeviani
##-----------------------------------------------------------------------------
## Definições da sessão.
## rm(list=ls())
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "plyr")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra doBy multcomp plyr
## TRUE TRUE TRUE TRUE TRUE
sessionInfo()
## R version 3.1.0 (2014-04-10)
## 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] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] plyr_1.8.1 multcomp_1.3-3 TH.data_1.0-3 mvtnorm_0.9-9997
## [5] doBy_4.5-10 MASS_7.3-33 survival_2.37-7 latticeExtra_0.6-26
## [9] RColorBrewer_1.0-5 lattice_0.20-29 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.0 lme4_1.1-6
## [5] Matrix_1.1-3 methods_3.1.0 minqa_1.2.3 nlme_3.1-117
## [9] Rcpp_0.11.2 RcppEigen_0.3.2.0.2 sandwich_2.3-0 stringr_0.6.2
## [13] tools_3.1.0 zoo_1.7-10
## trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/latencia.txt",
header=TRUE, sep="\t")
da <- na.omit(da)
str(da)
## 'data.frame': 45 obs. of 7 variables:
## $ cor0 : num 79.4 56.4 65 73.1 102.7 ...
## $ ida0 : num 1.51 1.36 1.2 1.41 1.57 1.43 1.13 1.63 1.55 1.31 ...
## $ ms0 : num 10.2 12.3 15 13.2 11.4 ...
## $ ss0 : num 7.01 5.7 5.81 7.26 7.57 7.09 6.06 7.5 6.17 6.34 ...
## $ f0 : num 5.89 5.55 6.07 5.93 5.57 5.91 5.53 5.98 5.63 5.48 ...
## $ b0 : num 12.3 11.8 13.6 13.4 12.2 ...
## $ lat48: int 0 1 1 0 0 1 1 0 0 1 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:11] 2 8 10 19 23 26 31 36 48 51 ...
## .. ..- attr(*, "names")= chr [1:11] "2" "8" "10" "19" ...
##-----------------------------------------------------------------------------
## Ver.
splom(da, type=c("p","smooth"))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- glm(lat48~b0+ida0+ms0+ss0, data=da, family="binomial")
## Diagnóstico nos resíduos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Estimativas dos parâmetros.
summary(m0)
##
## Call:
## glm(formula = lat48 ~ b0 + ida0 + ms0 + ss0, family = "binomial",
## data = da)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.129 -0.689 -0.265 0.715 1.811
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.730 6.511 2.42 0.0157 *
## b0 -1.937 0.661 -2.93 0.0034 **
## ida0 -3.226 3.243 -0.99 0.3198
## ms0 0.527 0.245 2.15 0.0313 *
## ss0 0.900 0.676 1.33 0.1834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 60.571 on 44 degrees of freedom
## Residual deviance: 40.542 on 40 degrees of freedom
## AIC: 50.54
##
## Number of Fisher Scoring iterations: 5
##-----------------------------------------------------------------------------
## Refinando para um número menor de variáveis.
m1 <- glm(lat48~(b0+ms0)^2, data=da, family="binomial")
summary(m1)
##
## Call:
## glm(formula = lat48 ~ (b0 + ms0)^2, family = "binomial", data = da)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.705 -0.705 -0.341 0.720 2.186
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.6998 11.4438 1.02 0.31
## b0 -1.4673 1.0168 -1.44 0.15
## ms0 0.6951 0.8501 0.82 0.41
## b0:ms0 -0.0151 0.0702 -0.22 0.83
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 60.571 on 44 degrees of freedom
## Residual deviance: 42.523 on 41 degrees of freedom
## AIC: 50.52
##
## Number of Fisher Scoring iterations: 5
m2 <- glm(lat48~b0+ms0, data=da, family="binomial")
summary(m2)
##
## Call:
## glm(formula = lat48 ~ b0 + ms0, family = "binomial", data = da)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.703 -0.743 -0.326 0.729 2.212
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 13.872 5.283 2.63 0.0086 **
## b0 -1.654 0.531 -3.12 0.0018 **
## ms0 0.517 0.173 2.98 0.0029 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 60.571 on 44 degrees of freedom
## Residual deviance: 42.571 on 42 degrees of freedom
## AIC: 48.57
##
## Number of Fisher Scoring iterations: 5
##-----------------------------------------------------------------------------
## Predição da curva de probabilidade de latência.
pred <- with(da,
expand.grid(ms0=seq(min(ms0),max(ms0),l=20),
b0=seq(min(b0),max(b0),l=20)))
pred$lat48 <- predict(m2, newdata=pred, type="response")
##-----------------------------------------------------------------------------
## Gráfico.
p1 <-
wireframe(lat48~ms0+b0, data=pred, drape=TRUE, scales=list(arrows=FALSE),
screen=list(z=60, x=-60), zlim=c(0,1))
p2 <- levelplot(lat48~ms0+b0, data=pred, aspect=1)+
as.layer(contourplot(lat48~ms0+b0, data=pred, cuts=20))
plot(p1, split=c(1,1,2,1), more=TRUE)
plot(p2, split=c(2,1,2,1), more=FALSE)
##-----------------------------------------------------------------------------
## Aprimorando o gráfico.
panel.3d.contour <- function(x, y, z, rot.mat, distance, type="on",
nlevels=20, zlim.scaled, col.contour=1, ...){
clines <- contourLines(x, y, matrix(z, nrow=length(x), byrow=TRUE), nlevels=nlevels)
if(any(type%in%c("bottom"))){
for(ll in clines){
n <- ltransform3dto3d(rbind(ll$x, ll$y, zlim.scaled[1]), rot.mat, distance)
panel.lines(n[1,], n[2,], col=col.contour, lty=1, lwd=1)
}}
panel.3dwire(x, y, z, rot.mat, distance, zlim.scaled=zlim.scaled, ...)
if(any(type%in%c("on"))){
for(ll in clines){
n <- ltransform3dto3d(rbind(ll$x, ll$y, ll$level), rot.mat, distance)
panel.lines(n[1,], n[2,], col=col.contour, lty=1, lwd=1)
}}
if(any(type%in%c("top"))){
for(ll in clines){
n <- ltransform3dto3d(rbind(ll$x, ll$y, zlim.scaled[2]), rot.mat, distance)
panel.lines(n[1,], n[2,], col=col.contour, lty=1, lwd=1)
}}
}
colr <- brewer.pal(11, "RdBu")
colr <- colorRampPalette(rev(colr), space="rgb")
wireframe(lat48~ms0+b0, data=pred, drape=TRUE, scales=list(arrows=FALSE),
screen=list(z=60, x=-60), zlim=c(0,1), col.regions=colr(100),
nlevels=60, col=rgb(0.7,0.7,0.7,1), col.contour="gray30",
panel.3d.wireframe="panel.3d.contour", type=c("bottom"),
par.settings=list(regions=list(alpha=0.75)),
xlab="Matéria seca\n na inoculação",
ylab="Grau brix\n na inoculação",
zlab=list("Probabilidade de latência antes de 48 horas", rot=90))
