#-----------------------------------------------------------------------
# Load packages.
library(lattice)
library(latticeExtra)
library(car)
library(nlstools)
library(nls2)
library(rootSolve) # gradient().
library(XML) # readHTMLTable().
library(lubridate) # hms(), hour(), minute(), second()
# url <- "http://nls2.googlecode.com/svn-history/r4/trunk/R/as.lm.R"
# download.file(url, dest = basename(url))
# path <- ifelse(Sys.info()["user"] == "walmes", basename(url), url)
# source(path)
source("~/Dropbox/CursoR/GeneticaEsalq/as.lm.R")
Os dados utilizados nessas análises foram extraídos do site http://www.4mula1.ro/. Esse site divulga resultados de corridas de fórmula 1. Os dados referem-se ao tempo do primeiro colocado da corrida em Monza desde 1950.
#-----------------------------------------------------------------------
# Tempo do primeiro lugar nas provas de Monza desde 1950.
# Lendo tabela direto do site.
da <- readHTMLTable("http://www.4mula1.ro/s/?cdx=9b0870e&pos=1&idt=44")
da <- da[[1]]
# str(da)
head(da)
## # Date GrandPrix Driver Team Pos Pol
## 1 1 1950-09-03 Monza Giuseppe Farina Alfa Romeo 1 3
## 2 2 1951-09-16 Monza Alberto Ascari Ferrari 1 3
## 3 3 1952-09-07 Monza Alberto Ascari Ferrari 1 1
## 4 4 1953-09-13 Monza Juan Manuel Fangio Maserati 1 2
## 5 5 1954-09-05 Monza Juan Manuel Fangio Mercedes 1 1
## 6 6 1955-09-11 Monza Juan Manuel Fangio Mercedes 1 1
## Points FL Time
## 1 8.00 2:51:17.400
## 2 8.00 2:42:39.300
## 3 8.50 2:50:45.600
## 4 9.00 2:49:45.900
## 5 8.00 2:47:47.900
## 6 8.00 2:25:04.400
# Provas em que Ayrton Senna ganhou.
subset(da, grepl("Senna", Driver))
## # Date GrandPrix Driver Team Pos Pol Points FL
## 40 40 1990-09-09 Monza Ayrton Senna McLaren 1 1 9.00
## 42 42 1992-09-13 Monza Ayrton Senna McLaren 1 2 10.00
## Time
## 40 1:17:57.878
## 42 1:18:15.349
# Anos das corridas.
da$ano <- as.integer(sub("^(\\d{4})-.*$", "\\1", da$Date))
# Tempo do primeiro colocado na prova.
x <- hms(da$Time)
da$tpro <- 60 * hour(x) + minute(x) + second(x)/60
# da <- read.table("monza.txt", header = TRUE, sep = ";",
# stringsAsFactors = FALSE, dec = ",")
# str(da)
#-----------------------------------------------------------------------
xyplot(tpro ~ ano, data = da, type = c("p", "smooth"))
# Cria uma variável que começa em 0.
da$anoc <- da$ano - min(da$ano)
plot(tpro ~ anoc, data = da)
start <- list(b0 = 160, b1 = -3, xb = 25)
with(start, {
curve(b0 + b1 * x * (x <= xb) + b1 * xb * (x > xb),
add = TRUE, col = 2)
})
# Ajuste do modelo não linear.
n0 <- nls(tpro ~ b0 +
b1 * anoc * (anoc <= xb) +
xb * b1 * (anoc > xb),
data = da, start = start)
# Resumo do ajuste.
summary(n0)
##
## Formula: tpro ~ b0 + b1 * anoc * (anoc <= xb) + xb * b1 * (anoc > xb)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 173.3307 2.7459 63.12 <2e-16 ***
## b1 -3.7756 0.1884 -20.04 <2e-16 ***
## xb 25.0115 0.7878 31.75 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.204 on 63 degrees of freedom
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 4.761e-09
# Intervalos de confiança.
# confint(n0)
confint.default(n0)
## 2.5 % 97.5 %
## b0 167.948845 178.712505
## b1 -4.144778 -3.406397
## xb 23.467393 26.555556
# Curva ajustada.
plot(tpro ~ anoc, data = da)
with(as.list(coef(n0)), {
curve(b0 + b1 * x * (x <= xb) + b1 * xb * (x > xb),
add = TRUE,
col = 4)
})
#-----------------------------------------------------------------------
# Análise dos resíduos.
m0 <- as.lm(n0)
par(mfrow = c(2, 2))
plot(m0)
layout(1)
#-----------------------------------------------------------------------
# Bandas de confiança.
# Modelo escrito como função dos parâmetros (theta).
f <- function(theta, a) {
with(as.list(theta),
b0 +
b1 * a * (a <= xb) +
xb * b1 * (a > xb))
}
range(da$anoc)
## [1] 0 66
pred <- data.frame(anoc = sort(c(coef(n0)["xb"] + c(0, 0.1),
seq(min(da$anoc), max(da$anoc),
length.out = 30))))
# pred$fit <- f(theta = coef(n0), a = pred$anoc)
pred$fit <- predict(n0, newdata = pred)
der <- gradient(f, x = coef(n0), a = pred$anoc)
str(der)
## num [1:32, 1:3] 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "b0" "b1" "xb"
F <- der
U <- chol(vcov(n0))
pred$se <- sqrt(apply(F %*% t(U), 1, function(x) sum(x^2)))
tval <- qt(p = c(lwr = 0.025, upr = 0.975), df = df.residual(n0))
me <- outer(pred$se, tval, "*")
pred <- cbind(pred, sweep(me, 1, pred$fit, "+"))
str(pred)
## 'data.frame': 32 obs. of 5 variables:
## $ anoc: num 0 2.28 4.55 6.83 9.1 ...
## $ fit : num 173 165 156 148 139 ...
## $ se : num 2.75 2.39 2.06 1.77 1.55 ...
## $ lwr : num 168 160 152 144 136 ...
## $ upr : num 179 170 160 151 142 ...
# source(paste0("https://raw.githubusercontent.com/",
# "walmes/wzRfun/master/R/panel.cbH.R"))
library(wzRfun)
## ----------------------------------------------------------------------
## wzRfun: Walmes Zeviani functions
##
## For collaboration, support or bug report, visit:
## https://github.com/walmes/wzRfun
##
## wzRfun version 0.70 (build in 2016-07-29) was loaded.
## ----------------------------------------------------------------------
xyplot(tpro ~ anoc, data = da) +
as.layer(xyplot(fit ~ anoc, data = pred,
type = "l",
prepanel = prepanel.cbH,
cty = "bands",
ly = pred$lwr,
uy = pred$upr,
panel = panel.cbH))
## Updated in 2016-10-10.
##
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods stats graphics grDevices utils datasets
## [7] base
##
## other attached packages:
## [1] wzRfun_0.70 lubridate_1.5.6 XML_3.98-1.4
## [4] rootSolve_1.6.6 nls2_0.2 proto_0.3-10
## [7] nlstools_1.0-2 car_2.1-2 latticeExtra_0.6-28
## [10] RColorBrewer_1.1-2 lattice_0.20-33 rmarkdown_1.0
## [13] knitr_1.14
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.7 formatR_1.4 nloptr_1.0.4
## [4] tools_3.3.1 digest_0.6.10 lme4_1.1-12
## [7] evaluate_0.9 nlme_3.1-128 mgcv_1.8-13
## [10] Matrix_1.2-6 rpanel_1.1-3 yaml_2.1.13
## [13] parallel_3.3.1 mvtnorm_1.0-5 SparseM_1.7
## [16] stringr_1.0.0 MatrixModels_0.4-1 grid_3.3.1
## [19] nnet_7.3-12 survival_2.39-4 multcomp_1.4-6
## [22] TH.data_1.0-7 minqa_1.2.4 magrittr_1.5
## [25] codetools_0.2-14 htmltools_0.3.5 MASS_7.3-45
## [28] splines_3.3.1 pbkrtest_0.4-6 quantreg_5.26
## [31] sandwich_2.3-4 stringi_1.1.1 doBy_4.5-15
## [34] zoo_1.7-13