Modelo Linear-Platô aos Tempos de Prova de Fórmula 1

Carregando os Dados

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))

Modelo Linear-Platô aos Tempos de Prova de Fórmula 1

Walmes Marques Zeviani

Carregando Pacotes

Carregando os Dados

Informações da Sessão