Estudo teórico da relação média variância

Em um primeiro momento, simplesmente iremos estudar a relação teórica entre \(V[Y]\) e \(E[Y]\). Para isso, basta estabelecer um grid de valores para \(E[Y]\) e exponenciar tais valores de acordo com o parâmetro de potência (\(p\)). Para tal exemplo, iremos utilizar:

• \(E[Y]\) = 10, 11, 12, …, 100;
• \(p\) = 0, 1, 1.5 e 2;
• \(\phi = 1\)

Não será utilizado \(p\) = 3 por conta da escala, que dificulta avaliar a relação média-variância das demais distribuições através de gráfico.

# Teórico
means <- 10:100
df <- data.frame(p0 = means^0,
                 p1 = means,
                 p1.5 = means^1.5,
                 p2 = means^2,
                 means = means)
df.l <- reshape2::melt(df, id.vars = "means")
levels(df.l$variable) <- c("p = 0 (normal)", "p = 1 (Poisson)", "p = 1.5 (Poisson-gamma)", "p = 2 (gamma)")
df.l$value.s <- sqrt(df.l$value)
df.l$means.s <- sqrt(df.l$means)

# Raiz quadrada no braço
g0 <- ggplot(data = df.l, aes(x = means.s, y = value.s, colour = variable)) +
        geom_line() +
        labs(x = expression(sqrt("E[Y]")),
             y = expression(sqrt("V[Y]")),
             title = "Relação média e variância teórico") +
        scale_y_continuous(breaks = c(0, 4, 6, 8, 10, 25, 31.6, 50, 75, 100, 105))

Forma da densidade

Uma vez que já entendemos a relação média variância, estudaremos a forma da curva da distribuição Tweedie, isto é, a sua densidade. Para tal, calcularemos o valor da densidade da distribuição Tweedie via tweedie::dtweedie (e quando \(p = 0\) utilizaremos dnorm) fixando \(\mu = 10\), dispersão (\(\phi = 3\)) e Y limitado a 60. Aqui incluímos a distribuição Normal Inversa e removemos a Poisson (pois não é contínua). A próxima figura apresenta a comparação das 4 distribuições:

end <- 80
mean <- 10
phi <- 2
a <- curve(dnorm(x, mean = mean, sd = phi), from = -3, end)
b <- curve(dtweedie(x, power = 1.5, mu = mean, phi = phi), from = 0,to = end, col = "purple", add = T)
c <- curve(dtweedie(x, power = 2, mu = mean, phi = phi), from = 0.01,to = end, col = "red", add = T)
d <- curve(dtweedie(x, power = 3, mu = mean, phi = phi), from = 0.01,to = end, col = "blue", add = T)
a$y <- a$y/max(a$y)
b$y <- b$y/max(b$y)
c$y <- c$y/max(c$y)
d$y <- d$y/max(d$y)

plot(a, type = "l", ylab = "Densidade", xlab = "Y")
legend("topright",
       col = c("black", "purple", "red", "blue"),
       legend = c("p = 0 (normal)", "p = 1.5 (gamma-Poisson)","p = 2 (gamma)", "p = 3 (normal inversa)"),
       lty = 1)
lines(b, col = "purple")
lines(c, col = "red")
lines(d, col = "blue")

# dtweedie(0, power = 1.5, mu = 2, phi = 5)

Uma vez que a distribuição normal é simétrica, a sua densidade está em torno da média, conforme esperado. A gamma-Poisson por sua vez chega a possuir um valor de densidade para Y = 0 (enquanto as outras não) e possui caráter assimétrico possuindo maior massa para valores ao redor de 5. A distribuição gamma tem maior valor de densidade para valores pequenos (até Y = 5 aproximadamente) e uma cauda pesada para os demais valores. Já a normal inversa assume maiores valores de densidade até Y = 10 em relação a distribuição gamma e cauda pesada para os demais valores de Y. Essa cauda é mais pesada que a da distribuição gamma, para isso, é necessário limitar os limites do eixo Y (isto é, dar um “zoom” no gráfico):

plot(a, type = "l", ylab = "Densidade", xlab = "Y", ylim = c(0, 0.1))
legend("topright",
       col = c("black", "purple", "red", "blue"),
       legend = c("p = 0 (normal)", "p = 1.5 (gamma-Poisson)","p = 2 (gamma)", "p = 3 (normal inversa)"),
       lty = 1)
lines(b, col = "purple")
lines(c, col = "red")
lines(d, col = "blue")

Estimação via Máxima Verossimilhança

Até aqui, foi apresentada a distribuição Tweedie a partir das suas propriedades e seus domínios, que são necessários para compreender como ela funciona na prática. Ao ajustar um modelo de regressão Tweedie no contexto de GLM, o primeiro passo é estimar o valor de \(p\), isto é, é necessário identificar qual forma da distribuição. Uma forma de se fazer isso é via máxima verossimilhança, e esse será o primeiro método abordado. Para isso, utilizaremos a função tweedie::tweedie.profile para traçar o perfil da verossimilhança para diferentes valores de \(p\) (parâmetro potência) e escolher o valor de \(p\) que maximize a função de log verossimilhança (\(l(x)\)). Em seguida, após identificarmos o valor de p que maximiza a função de log verossimilhança, iremos ajustar um GLM tweedie utilizando a família Tweedie, disponível no pacote statmod.

Para isso, foram simuladas 1000 amostras de \(Y \sim G(\alpha = 1, \beta = \frac{1}{0.5})\) sendo o primeiro o parâmetro de forma (\(\alpha = shape\)) e o segundo o parâmetro de taxa (\(\beta = rate = \frac{1}{escala }\)), cuja representação da distribuição gamma foi dada da seguinte forma: \[f(y) = \frac{\beta^{\alpha}y^{a-1}e^{-\beta y}}{\Gamma(\alpha)}\] e têm-se as seguintes propriedades:

\(E[Y] = \frac{\alpha}{\beta}\) \(V[Y] = \frac{\alpha}{\beta^2}\)

Nesse caso:

set.seed(323)
y <- rgamma(1000, shape = 1, rate = 1/0.5)
mean(y)

[1] 0.53276

var(y)

[1] 0.26624

ou seja:

\[E[Y] = \frac{\alpha}{\beta} = \frac{1}{1/0.5} = 0.5 = \widehat{\bar{Y}} = 0.53\] \[V[Y] = \frac{\alpha}{\beta^2} = \frac{1}{(1/0.5)^2} = \frac{1}{4} = 0.25 = \widehat{\sigma}^2_{Y} = 0.27\]

valores simulados muito próximos dos valores dos parâmetros.

Uma vez que os dados foram simulados, é necessário estimar o parâmetro de pontência (\(p\)), que deve ser próximo de 2, uma vez que p = 2 representa a distribuição gamma. Para realizar a estimação de p, é necessário fornecer os seguintes parâmetros:

• form = nesse caso é a fórmula da regressão. Como não temos nenhuma covariável nesse exemplo, apresentaremos apenas \(y \sim 1\), caso contrário deveríamos informar as covariáveis.

• p.vec = um vetor de possíveis valores para o \(p\), nesse caso, será fornecido o valor 0 (que representa uma distribuição Normal), e uma sequência de 20 valores equidistantes entre o valor 1.1 e 5 (para melhor estudar o grau de assimetria dos dados;

• link.power = será passado o valor 0, que representa o link log;

• do.plot = será passado o valor T, indicado a função para plotar o perfil de verossimilhança;

• do.smooth = será passado o valor F, pois não quero uma curva suave nos dados nesse primeiro momento.

Segue o gráfico da estimação de \(p\):

adj1 <- tweedie::tweedie.profile(y ~ 1,
                         p.vec = c(0, seq(1.1, 3, length.out = 15)),
                         link.power = 0,
                         do.plot = T,
                         do.smooth = F,
                         method = "series")

0 1.1 1.2357 1.3714 1.5071 1.6429 1.7786 1.9143 2.05 2.1857 2.3214 2.4571 2.5929 2.7286 2.8643 3 
................Done.

Note que para \(p=0\) temos um valor bem baixo de \(l(p)\), o que já era de se esperar, uma vez que os dados foram gerados de uma distribuição gamma. A partir do perfil da verossimilhança, nota-se que o valor de p que maximiza \(l(p)\) encontra-se próximo de 2. Ainda, o método não convergiu para valores de p superiores de 3, e ainda retornou avisos (warnings messages). Por um primeiro momento, tal aviso pode ser um pouco “intimidador”, no entanto, ao inspecionar o gráfico anterior, nota-se que não faz sentido valores de p > 3 para essa variável, e isso pode ter sido o motivo da não estimação de \(p\). A estimação do valor de \(p\) não é tão simples computacionalmente via MLE, e pode não ser estimada dependendo do valor de \(p\), e para isso, é necessário sempre checar o gráfico do perfil de \(l(p)\).

Nesse caso, faremos uma nova simulação da seguinte forma:

• Valores de \(p\) mais próximos de 2;
• Método de estimação series para que a estimação de p seja mais rápida;
• Solicitando o intervalo de confiança (para isso, é necessário solicitar uma curva smooth). Argumentos do.smooth = T e do.ci = T.

Segue resultado:

adj1 <- tweedie::tweedie.profile(y ~ 1,
                         p.vec = seq(1.5, 2.5, length.out = 10),
                         link.power = 0,
                         do.plot = T,
                         do.smooth = T,
                         do.ci = T,
                         method = "series")

1.5 1.6111 1.7222 1.8333 1.9444 2.0556 2.1667 2.2778 2.3889 2.5 
..........Done.

A linha horizontal pontilhada nesse caso apresenta o intervalo de confiança de 95% para o valor de \(p\). Os pontos do gráficos foram os valores de p utilizados, e a linha que conecta os pontos foi feita a partir da suavização da curva via splines.

O valor de p e do intervalo são encontrados respectivamente por:

# Estimativa pontual de p
adj1$p.max

[1] 2.0102

# Intervalo de confiança de p
adj1$ci

[1] 1.9281 2.0787

Embora o valor p = 2.010204 não tenha sido utilizado para estimar o valor de \(l(p)\), a partir da suavização foi o valor que maximiza \(l(p)\), e o mesmo será utilizado para estimar os parâmetros de regressão.

A estimação de um glm através da distribuição Tweedie é simples. Basta no argumento family especificar a função de ligação e o valor do parâmetro de potência:

m1 <- glm(y ~ 1, family = statmod::tweedie(var.power = adj1$p.max, link.power = 0))
summary(m1)

Call:
glm(formula = y ~ 1, family = statmod::tweedie(var.power = adj1$p.max, 
    link.power = 0))

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-3.743  -1.036  -0.315   0.333   2.798  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.6297     0.0306   -20.6   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Tweedie family taken to be 0.94407)

    Null deviance: 1116.3  on 999  degrees of freedom
Residual deviance: 1116.3  on 999  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

par(mfrow = c(1,3))
plot(m1, which = c(1:3))

É lógico que nesse modelo o ajuste seria feito sem mais complicações. A única “desvantagem” de usar essa metodologia, é que nem todos os demais pacotes do R tem funções implementadas para a família de distribuição Tweedie. Por exemplo, fazer half-normal plots não é possível. No entanto, como o valor de p foi praticamente da distribuição Gamma, basta ajustar um modelo gamma e realizar demais verificações do modelo:

m1 <- glm(y ~ 1, family = Gamma(link = "log"))
par(mfrow = c(1,1))
hnp::hnp(m1)

Gamma model 

y.norm <- rnorm(100, mean = 20)
adj2 <- tweedie::tweedie.profile(y.norm ~ 1,
                                 p.vec = c(0, seq(1.1, 5, length.out = 30)),
                                 link.power = 0,
                                 do.plot = T,
                                 do.smooth = T,
                                 do.ci = T,
                                 method = "series")

0 1.1 1.2345 1.369 1.5034 1.6379 1.7724 1.9069 2.0414 2.1759 2.3103 2.4448 2.5793 2.7138 2.8483 2.9828 3.1172 3.2517 3.3862 3.5207 3.6552 3.7897 3.9241 4.0586 4.1931 4.3276 4.4621 4.5966 4.731 4.8655 5 
...............................Done.

# Pacote CPLM
y.tw <- tweedie::rtweedie(100, power = 1.05, phi = 1, mu = 4)
# plot(density(y.tw))
adj2 <- tweedie::tweedie.profile(y.tw ~ 1,
                                 p.vec = c(0, seq(1.01, 2, length.out = 20)),
                                 link.power = 0,
                                 do.plot = T,
                                 do.smooth = T,
                                 do.ci = T,
                                 method = "series")

When the response variable contains exact zeros, all values of p must be between 1 and 2; other values have been removed.
1.01 1.0621 1.1142 1.1663 1.2184 1.2705 1.3226 1.3747 1.4268 1.4789 1.5311 1.5832 1.6353 1.6874 1.7395 1.7916 1.8437 1.8958 1.9479 
...................Done.

Estimação via quase-verossimilhança

A estimação do parâmetro de potência \(p\) agora será tratado num contexto de quase-verossimilhança. Aqui, não é necessário traçar o perfil da verossimilhança (ou da pseudo-verossimilhança no caso) em uma primeira etapa, para então fornecer o valor de \(p\) para o ajuste. Para isso, irei utilizar a função mcglm::mcglm que estima o valor de \(p\) automaticamente no processo de estimação dos parâmetros do modelo (BONAT, 2016). A grande versatilidade do pacote é que permite ajustar um modelo multivarido (múltiplas respostas) com diferentes estruturas de média (generalizado) para dados correlacionados (dados longitudinais/temporais/espaciais/medidas repetidas) (BONAT; JØRGENSEN, 2016). A única diferença, é que é necessário incluir uma covariável correlacionada com a variável resposta para estimar o parâmetro de potência, pois não sabe-se a forma da densidade, como no caso da estimação via verossimilhança.

Sendo assim, foi gerado a variável resposta Y com distribuição Gamma condicional a média X:

\[Y\mid x \sim G(\mu(x), \phi)\]

E a média foi gerada a partir de uma equação de regressão, isto é:

\[ln(\mu) = \beta_{0} + \beta_{1}x\]

em que X:

\[X\sim U (0, 1)\]

Segue o gráfico de dispersão entre as variáveis Y e X:

set.seed(323)
x <- runif(200)
b0 <- 2;b1 <- 0.5;phi <- 1.3
set.seed(888)
y <- rgamma(200, shape = exp(b0 + b1*x), rate = 1/phi)
# mean(y)
# mean(exp(b0 + b1*x))/(1/phi)
data <- data.frame(y = y,
                   x = x)
with(data, plot(y ~ x))

# abline(a = 2, b = 0.5, col = "blue")
# abline(coef = exp(c(2.28, 0.47)), col = "red")

Nota-se uma relação positiva entre ambas variáveis (uma vez que foi induzida pelo modelo). Agora, será ajustado um modelo de regressão de Y em função de X. Para tal, iremos especificar os seguintes parâmetros:

• matrix_pred: Matriz que identifica a estrutura de covariância nos dados. Nesse caso, é apenas uma matriz diagonal identidade (pois não tem estrutura de dependência dos dados).

• control_algorithm: Controla a estimação dos parâmetros do modelo de regressão através algoritimo. Tuning = 0.35 faz com que o algoritmo seja estudada de forma mais lenta, no entanto, permite o ajuste; Max_iter = 0.5 permite o algoritmo iterar mais vezes que o padrão (20);

• variance: Estrutura de variância para o modelo, nesse caso Tweedie;

• power_fixed: Para o modelo estimar o parâmetro de potência da Tweedie, é necessário deixar como F;

• link: Função de ligação utilizada para a média, nesse caso, link log.

Segue o summary do modelo ajustado:

z0 <- mcglm::mc_id(data)
m1 <- mcglm::mcglm(c(y ~ x),
            matrix_pred = list(z0),
            link = "log",
            variance = "tweedie",
            power_fixed = F,
            data = data,
            control_algorithm = list(tuning = 0.35, max_iter = 50))

Automatic initial values selected. 

summary(m1)

Call: y ~ x

Link function: log
Variance function: tweedie
Covariance function: identity
Regression:
            Estimates Std.error Z value
(Intercept)   2.28347  0.042772 53.3867
x             0.46811  0.072647  6.4436

Power:
  Estimates Std.error Z value
1    2.0549   0.66198  3.1041

Dispersion:
  Estimates Std.error Z value
1  0.081936   0.13502 0.60682

Algorithm: chaser
Correction: TRUE
Number iterations: 22

As estimativas dos parâmetros de regressão foram próximos aos induzidos pelo modelo (\(\beta_{0} = 2;\beta_{1} = 0.5\)). O parâmetro de potência foi igual a 2.05, valor esse próximo de 2. As interpretações dos parâmetros são feitos de forma análoga a da função glm utilizando family = Gamma(link = "log"). Para validar o modelo, segue análise de resíduos:

plot(m1)

Não nota-se relação média e variância nos resíduos em relação aos valores ajustados. O fato de uma pequena fuga de normalidade nesse gráfico não é um problema, uma vez que os resíduos não seguem uma distribuição normal (uma vez que a distribuição utilizada é a gamma, assimétrica).

Aplicação da Tweedie em modelos GAM

Contexto e metodologia

Uma outra utilidade de tal distribuição é na abordagem dos modelos aditivos generalizados (GAM). Imagine que tenhamos um problema em que eu tenho uma variável resposta positiva assimétrica que contém 0 e outra variável numérica qualquer e desejo verificar se as duas estão correlacionadas. Após fazer o gráfico de dispersão, nota-se que a relação entre as duas é não linear, inviabilizando assim, o cálculo do coeficiente de correlação de Pearson, que já foi apresentado nesse post. Como alternativa, é possível estimar uma curva suave para a relação entre y e x a partir da função mgcv::gam. Ela não retorna um coeficiente de correlação em si, no entanto permite visualizar graficamente tal relação. Os modelos GAM estendem os modelos GLM no sentido de ser possível estimar uma relação não linear entre uma covariável numérica e uma preditora numérica. Essa função suave, por exemplo, pode ser estimada tanto por polinômios locais, splines e densidade de kernel. A escolha da função suave normalmente é feita internamente pelo pacote através da função s(), entretanto, nesse exemplo foi necessário calibrar alguns parâmetros.

Aplicação

O conjunto de dados utilizado nesse exemplo foi de Swedish Motor Insurance que foi publicado em FREES (2009) e podemos baixá-lo aqui, cuja documentação está diponível aqui. Basicamente, estamos interessado em modelar o total de pagamentos feita pela companhia de seguros (Payment) em relação o número de anos do segurados (insured).