1 Introducao

Este projeto apresenta a modelagem da diferenca de gols entre mandantes e visitantes no Campeonato Brasileiro. A variavel resposta considerada e a diferenca entre o total de gols marcados pelos mandantes e o total de gols marcados pelos visitantes em cada rodada.

Como essa diferenca assume valores inteiros positivos, nulos ou negativos, foram avaliadas distribuicoes discretas capazes de representar esse tipo de variavel. Os modelos considerados foram:

Normal discreta;
skew-normal discreta;
Skellam;
diferenca de duas Weibull discretas.

Os parametros foram estimados por maxima verossimilhanca e os modelos foram comparados pelos criterios AIC e BIC, alem da inspecao grafica das probabilidades ajustadas.

2 Pacotes e configuracoes

knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.align = "center"
)

options(OutDec = ",", scipen = 99, digits = 4)

library(extraDistr)
library(sn)
library(moments)
library(dplyr)
library(knitr)

3 Importacao e organizacao dos dados

A base utilizada deve estar no mesmo diretorio do arquivo R Markdown. Neste script, considera-se o arquivo tudo.txt, com informacoes agregadas por ano e rodada.

Foram mantidos apenas os campeonatos a partir de 2006, periodo em que o Campeonato Brasileiro passou a apresentar formato mais estavel, com 20 equipes e 38 rodadas.

dados <- read.table("tudo.txt", header = TRUE)

dados <- subset(dados, ano >= 2006)

head(dados)

##      ano rodata mandante_Placar visitante_Placar soma_diferenca jogos
## 134 2006      1              11                7              4    10
## 135 2006      2              21                5             16    10
## 136 2006      3              18               15              3    10
## 137 2006      4              15               10              5    10
## 138 2006      5              14               11              3    10
## 139 2006      6              15               15              0    10

str(dados)

## 'data.frame':    760 obs. of  6 variables:
##  $ ano             : int  2006 2006 2006 2006 2006 2006 2006 2006 2006 2006 ...
##  $ rodata          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ mandante_Placar : int  11 21 18 15 14 15 25 11 12 12 ...
##  $ visitante_Placar: int  7 5 15 10 11 15 13 10 9 12 ...
##  $ soma_diferenca  : int  4 16 3 5 3 0 12 1 3 0 ...
##  $ jogos           : int  10 10 10 10 10 10 10 10 10 10 ...

A variavel de interesse e a diferenca de gols entre mandante e visitante. No codigo original, essa variavel se encontra na quinta coluna da base.

x <- as.numeric(as.character(dados[, 5]))

n <- length(x)
xs <- seq(min(x), max(x), by = 1)

summary(x)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -13,00    2,00    5,00    4,89    8,00   22,00

4 Analise descritiva

Antes do ajuste dos modelos, sao obtidas medidas descritivas da diferenca de gols por rodada.

descritiva_x <- data.frame(
  n = length(x),
  media = mean(x),
  mediana = median(x),
  desvio_padrao = sd(x),
  minimo = min(x),
  primeiro_quartil = quantile(x, 0.25),
  terceiro_quartil = quantile(x, 0.75),
  maximo = max(x),
  assimetria = skewness(x),
  curtose = kurtosis(x)
)

kable(descritiva_x, digits = 3)

	n	media	mediana	desvio_padrao	minimo	primeiro_quartil	terceiro_quartil	maximo	assimetria	curtose
25%	760	4,888	5	4,912	-13	2	8	22	0,017	3,355

5 Definicao das distribuicoes

Nesta etapa, sao definidas as funcoes de probabilidade utilizadas no ajuste dos modelos.

5.1 Normal discreta

A normal discreta e obtida pela diferenca entre valores acumulados da distribuicao normal continua:

\[ P(X=x)=\Phi\left(\frac{x-\mu+1}{\sigma}\right) - \Phi\left(\frac{x-\mu}{\sigma}\right). \]

ddnorm2 <- function(x, mu, sigma) {
  pnorm((x - mu + 1) / sigma) -
    pnorm((x - mu) / sigma)
}

5.2 Skew-normal discreta

A skew-normal discreta e definida de forma analoga, pela diferenca entre valores acumulados da distribuicao skew-normal continua.

ddskewnorm <- function(x, xi, omega, alpha) {
  psn(x + 1, xi = xi, omega = omega, alpha = alpha) -
    psn(x, xi = xi, omega = omega, alpha = alpha)
}

5.3 Skellam

A distribuicao Skellam representa a diferenca entre duas variaveis aleatorias independentes com distribuicao de Poisson. Se \(X \sim Poisson(\lambda_1)\) e \(Y \sim Poisson(\lambda_2)\), entao \(D = X - Y\) segue distribuicao Skellam.

dskellam <- function(x, lambda1, lambda2) {
  exp(-(lambda1 + lambda2)) *
    (lambda1 / lambda2)^(x / 2) *
    besselI(2 * sqrt(lambda1 * lambda2), abs(x))
}

5.4 Diferenca de Weibull discretas

Tambem foi considerada a distribuicao da diferenca entre duas variaveis independentes com distribuicao Weibull discreta.

dweibull_discrete <- function(x, q, beta) {
  if (any(x < 0)) {
    stop("x deve ser >= 0")
  }

  q^(x^beta) - q^((x + 1)^beta)
}

dweibull_diff <- function(d, q1, beta1, q2, beta2, max_k = 200) {
  sapply(d, function(di) {
    if (di >= 0) {
      sum(
        dweibull_discrete(0:max_k + di, q1, beta1) *
          dweibull_discrete(0:max_k, q2, beta2)
      )
    } else {
      sum(
        dweibull_discrete(0:max_k, q1, beta1) *
          dweibull_discrete(0:max_k - di, q2, beta2)
      )
    }
  })
}

6 Funcoes de log-verossimilhanca

As funcoes abaixo retornam a log-verossimilhanca negativa. Dessa forma, todas podem ser minimizadas diretamente pela funcao optim.

loglik_norm <- function(par) {
  mu <- par[1]
  sigma <- exp(par[2])

  p <- ddnorm2(x, mu, sigma)

  -sum(log(p + 1e-12))
}

loglik_skew <- function(par) {
  xi <- par[1]
  omega <- abs(par[2])
  alpha <- par[3]

  p <- ddskewnorm(x, xi, omega, alpha)

  -sum(log(p + 1e-12))
}

loglik_skellam <- function(par) {
  lambda1 <- exp(par[1])
  lambda2 <- exp(par[2])

  p <- dskellam(x, lambda1, lambda2)

  -sum(log(p + 1e-12))
}

loglik_weibull_diff <- function(par) {
  q1 <- 1 / (1 + exp(-par[1]))
  beta1 <- exp(par[2])
  q2 <- 1 / (1 + exp(-par[3]))
  beta2 <- exp(par[4])

  p <- dweibull_diff(x, q1, beta1, q2, beta2)

  -sum(log(p + 1e-12))
}

7 Ajuste dos modelos

Os parametros foram estimados por maxima verossimilhanca. Para garantir as restricoes dos parametros, algumas transformacoes foram utilizadas:

sigma = exp(par[2]), garantindo \(\sigma > 0\);
lambda = exp(par), garantindo \(\lambda > 0\);
q = 1/(1 + exp(-par)), garantindo \(0 < q < 1\);
beta = exp(par), garantindo \(\beta > 0\).

fit_norm <- optim(
  par = c(mean(x), log(sd(x))),
  fn = loglik_norm,
  method = "BFGS"
)

fit_skew <- optim(
  par = c(mean(x), sd(x), 0),
  fn = loglik_skew,
  method = "BFGS"
)

lambda1_ini <- (var(x) + mean(x)) / 2
lambda2_ini <- (var(x) - mean(x)) / 2

lambda1_ini <- max(lambda1_ini, 0.1)
lambda2_ini <- max(lambda2_ini, 0.1)

fit_skellam <- optim(
  par = c(log(lambda1_ini), log(lambda2_ini)),
  fn = loglik_skellam,
  method = "BFGS"
)

fit_weibull <- optim(
  par = c(0, 0, 0, 0),
  fn = loglik_weibull_diff,
  method = "BFGS"
)

8 Parametros estimados

mu_norm_hat <- fit_norm$par[1]
sigma_norm_hat <- exp(fit_norm$par[2])

xi_hat <- fit_skew$par[1]
omega_hat <- abs(fit_skew$par[2])
alpha_hat <- fit_skew$par[3]

lambda1_hat <- exp(fit_skellam$par[1])
lambda2_hat <- exp(fit_skellam$par[2])

par_hat <- fit_weibull$par

q1_hat <- 1 / (1 + exp(-par_hat[1]))
beta1_hat <- exp(par_hat[2])
q2_hat <- 1 / (1 + exp(-par_hat[3]))
beta2_hat <- exp(par_hat[4])

tabela_parametros <- data.frame(
  Modelo = c(
    "Normal discreta",
    "Skew-normal discreta",
    "Skellam",
    "Diferenca Weibull discreta"
  ),
  Parametro_1 = c(mu_norm_hat, xi_hat, lambda1_hat, q1_hat),
  Parametro_2 = c(sigma_norm_hat, omega_hat, lambda2_hat, beta1_hat),
  Parametro_3 = c(NA, alpha_hat, NA, q2_hat),
  Parametro_4 = c(NA, NA, NA, beta2_hat)
)

kable(tabela_parametros, digits = 4)

Modelo	Parametro_1	Parametro_2	Parametro_3	Parametro_4
Normal discreta	5,3882	4,900	NA	NA
Skew-normal discreta	3,2401	5,350	0,5822	NA
Skellam	14,4704	9,582	NA	NA
Diferenca Weibull discreta	0,9933	2,224	0,8279	1,252

9 Comparacao por AIC e BIC

Os criterios AIC e BIC foram utilizados para comparar os modelos ajustados. Menores valores indicam melhor equilibrio entre qualidade de ajuste e complexidade.

ll_norm <- -fit_norm$value
ll_skew <- -fit_skew$value
ll_skel <- -fit_skellam$value
ll_weibull <- -fit_weibull$value

k_norm <- 2
k_skew <- 3
k_skel <- 2
k_weibull <- 4

criterios <- data.frame(
  Modelo = c(
    "Normal discreta",
    "Skew-normal discreta",
    "Skellam",
    "Diferenca Weibull discreta"
  ),
  LogLik = c(ll_norm, ll_skew, ll_skel, ll_weibull),
  k = c(k_norm, k_skew, k_skel, k_weibull),
  AIC = c(
    -2 * ll_norm + 2 * k_norm,
    -2 * ll_skew + 2 * k_skew,
    -2 * ll_skel + 2 * k_skel,
    -2 * ll_weibull + 2 * k_weibull
  ),
  BIC = c(
    -2 * ll_norm + log(n) * k_norm,
    -2 * ll_skew + log(n) * k_skew,
    -2 * ll_skel + log(n) * k_skel,
    -2 * ll_weibull + log(n) * k_weibull
  )
)

criterios <- criterios %>%
  arrange(AIC)

kable(criterios, digits = 3)

Modelo	LogLik	k	AIC	BIC
Skellam	-2287	2	4578	4588
Normal discreta	-2288	2	4579	4588
Diferenca Weibull discreta	-2286	4	4580	4599
Skew-normal discreta	-2288	3	4581	4595

10 Probabilidades observadas e ajustadas

Nesta etapa, sao calculadas as frequencias relativas observadas e as probabilidades estimadas por cada modelo para os valores inteiros observados da diferenca de gols.

freq_obs <- table(factor(x, levels = xs))
prob_obs <- as.numeric(freq_obs) / length(x)

prob_norm <- ddnorm2(
  xs,
  mu = mu_norm_hat,
  sigma = sigma_norm_hat
)

prob_skew <- ddskewnorm(
  xs,
  xi = xi_hat,
  omega = omega_hat,
  alpha = alpha_hat
)

prob_skel <- dskellam(
  xs,
  lambda1 = lambda1_hat,
  lambda2 = lambda2_hat
)

prob_weib <- dweibull_diff(
  xs,
  q1 = q1_hat,
  beta1 = beta1_hat,
  q2 = q2_hat,
  beta2 = beta2_hat
)

11 Grafico dos modelos ajustados

O grafico compara as probabilidades observadas com as probabilidades ajustadas pelos modelos considerados.

simbolos <- c(16, 0, 1, 2, 3)

cores <- c(
  "black",
  "#0072B2",
  "#D55E00",
  "#009E73",
  "#CC79A7"
)

rotulos <- c(
  "Observado",
  "Normal discreta",
  "Skew-normal discreta",
  "Skellam",
  "Diferenca Weibull discreta"
)

probabilidades <- list(
  prob_obs,
  prob_norm,
  prob_skew,
  prob_skel,
  prob_weib
)

desenhar_grafico <- function() {
  par(
    family = "sans",
    mar = c(4.5, 5, 1, 1),
    oma = c(0, 0, 3, 0),
    mgp = c(3, 0.8, 0),
    las = 1,
    bty = "l",
    cex.axis = 1,
    cex.lab = 1.1
  )

  ymax <- max(unlist(probabilidades), na.rm = TRUE)

  plot(
    xs, prob_obs,
    type = "n",
    xlab = "Diferenca de gols por rodada",
    ylab = "Probabilidade",
    ylim = c(0, ymax * 1.10),
    yaxs = "i",
    xaxt = "n",
    yaxt = "n"
  )

  marcas_x <- pretty(range(xs), n = 7)
  marcas_x <- marcas_x[
    marcas_x >= min(xs) & marcas_x <= max(xs)
  ]

  marcas_y <- pretty(c(0, ymax * 1.10), n = 5)
  marcas_y <- marcas_y[
    marcas_y >= 0 & marcas_y <= ymax * 1.10
  ]

  abline(
    h = marcas_y,
    col = "grey90",
    lwd = 0.8
  )

  axis(1, at = marcas_x)

  axis(
    2,
    at = marcas_y,
    labels = formatC(
      marcas_y,
      format = "f",
      digits = 2,
      decimal.mark = ","
    )
  )

  for (i in 2:5) {
    points(
      xs, probabilidades[[i]],
      pch = simbolos[i],
      col = cores[i],
      cex = 1.15,
      lwd = 1.8
    )
  }

  points(
    xs, prob_obs,
    pch = simbolos[1],
    col = cores[1],
    cex = 0.85
  )

  legend(
    "top",
    inset = c(0, -0.12),
    legend = rotulos,
    pch = simbolos,
    col = cores,
    ncol = 3,
    pt.cex = 1,
    cex = 0.85,
    bty = "n",
    xpd = NA
  )
}

desenhar_grafico()

12 Exportacao do grafico

Os comandos abaixo exportam o grafico em PNG e PDF, em dimensoes adequadas para insercao em poster ou apresentacao.

png(
  filename = "grafico_ajuste_brasileirao.png",
  width = 39,
  height = 20,
  units = "cm",
  res = 300,
  pointsize = 26,
  bg = "white"
)

desenhar_grafico()
dev.off()

pdf(
  file = "grafico_ajuste_brasileirao.pdf",
  width = 39 / 2.54,
  height = 20 / 2.54,
  pointsize = 26,
  family = "Helvetica",
  bg = "white"
)

desenhar_grafico()
dev.off()

13 Consideracoes finais

A comparacao entre as distribuicoes permite avaliar qual modelo representa melhor a diferenca de gols por rodada no Campeonato Brasileiro. A Skellam e uma alternativa natural para esse problema, pois descreve a diferenca entre duas contagens independentes. Entretanto, a comparacao com outras distribuicoes discretas, como a normal discreta, a skew-normal discreta e a diferenca de Weibull discretas, permite avaliar se estruturas mais flexiveis oferecem melhor ajuste aos dados.

Os criterios AIC e BIC devem ser analisados em conjunto com a inspecao grafica, pois pequenas diferencas nos criterios podem indicar ajustes semelhantes em termos praticos.

Modelagem da Diferenca de Gols no Campeonato Brasileiro

Lucas Pereira Belo

12/10/2026