REGRESSÃO BETA E APLICAÇÕES

BETA REGRESSION AND APLICATIONS

Galarza, Christian

Mestrando em Estatística – IMECC / UNICAMP

RESUMO

Este trabalho apresenta uma introdução ao Modelo de Regressão Beta simples com dispersão variável calculando seus estimadores de Máxima Verossimilhança pelo Método Newton-Raphson e Scoring de Fisher. Além disso, tem-se a aplicação dos modelos nos dados Gasoline Yield Data de Prater (1956), onde se testam diferentes condições como modelos e funções de ligação, isto utilizando o pacote betareg do software estatístico R.

Palavras-chave: Modelo de regressão Beta simples, Modelo de regressão Beta com dispersão variável, Newton-Rapson, Scoring de Fisher, Aplicações.

ABSTRACT

This paper presents an introduction to the Simple and Varying dispersion Beta Regression Model calculating the Maximum likelihood estimators through Newton-Raphson and Fisher’s Scoring Methods. Futhermore, has an application of the models with the Gasoline Yield Data of Prater (1956), where we test different conditions like models and link functions, all using the betareg package of statistical R software.

Keywords: Simple Beta Regression, Varying dispersion beta regression model, Newton-Raphson, Fisher’s Scoring, Aplications.

2  

1. INTRODUÇÃO

Nos modelos de regressão, muitas vezes a variável y resposta é uma variável que apenas toma valores no intervalo contínuo [0,1], de modo que é necessário considerar algumas alternativas para ajustar o modelo de regressão sob esta condição. Alguns exemplos de variáveis são os seguintes:

Porcentagem do tempo dedicado a alguma atividade. Taxa de pobreza, taxa de desemprego. Pontuação de uma prova. Proporção de um composto químico numa mistura. Fração da superfície colhida de uma floreta.

Pode-se aplicar algumas transformações à variável de resposta, mas pelo geral se precisa supor normalidade. Aliás, às vezes os parâmetros não podem ser interpretados facilmente em termos da resposta original e não é viável arrumar problemas de heterocedasticidade ou assimetria. Uma melhor alternativa é supor que a variável de resposta y segue uma distribuição contínua com suporte (0,1), tal como a distribuição Beta, a qual é ligada às variáveis preditoras por uma função de ligação. Os modelos de regressão Beta usam um tipo de parametrização em termos de sua média e a precisão (dispersão) porque assim é muito mais fácil fazer as interpretações das estimativas. O modelo de regressão Beta fornece estimações precisas e seguras em relação aos coeficientes, independentemente da tendência dos dados (valores próximos de zero o de um) ou do tamanho da amostra. Vale mencionar que a Distribuição Beta não pertence à família exponencial de funções, pois sua distribuição não pode ser escrita da forma canônica; ela tem sua própria família que comtempla as distribuições Uniforme, Arco-Seno e Dirichlet segundo Sant’ Anna, Catten (2009). Em geral, o modelo de regressão Beta e muito útil para diversas aplicações práticas, além de que há um campo grande de pesquisa onde se tem desenvolvido muita literatura sobre este tema nos últimos anos. A Sua aplicação foi implementada no pacote betareg no software R (www.r-project.org).

3  

2. MODELO DE REGRESSÃO BETA

Uma variável aleatória Y segue uma Distribuição Beta com parâmetros , 0 denotado por ~ , se a distribuição de Y tem densidade:

; ,ΓΓ Γ

1 , 0 1

0,

sendo Γ ∙ a função gama e sua média e variância dadas por

1

Ferrari and Cribari-Neto (2004) apresentaram uma parametrização diferente para a função densidade dos modelos de regressão Beta, isto em termos da média e um parâmetro de precisão . Seja / e , i.e., e 1 ; sob a nova parametrização ~ , e sua densidade em (2.1) pode ser escrita como,

; ,Γ

Γ Γ 11 , 0 1

0,

com média e variância

11

Pode-se perceber que o parâmetro faz o papel de um parâmetro de precisão dado que para fixo, a maior valor de , menor variância da variável de resposta. Veja-se no gráfico 1

densidades betas com diferentes medias mas com mesma precisão, no quadro à esquerda com precisão de 5 e à direita com precisão de 15 unidades. Considere-se uma amostra independente , … . . , onde tem distribuição Beta com média

e precisão desconhecida , i.e., ~ , . Seja ∙ uma função 0,1 ⟼ , estritamente monótona e duas vezes diferençável. A função ∙ é chamada função de ligação dada por

(2.1)

(2.2)

(2.3)

(2.4)

(2.5)

4  

Gráfico 1. Beta densities for different values of( , ). Fonte: [4] Beta regression in R.

onde é o vetor de covariaveis de dimensão p e , …… . , é o vetor de parâmetros da regressão, i = 1, . . . , n. tem-se diferentes propostas de funções de ligação ∙ podendo ser escolhida a que leve para um bom ajuste. As funções de ligação mais comuns são as funções inversas acumuladas das distribuições Logística, Normal padrão, Valor mínimo extremo, Valor máximo extremo e Cauchy. Veja-se tabela 1 para mais detalhe.

DISTRIBUIÇÃO NOME FUNÇÃO DE LIGAÇÃO

Logística Logit log 1⁄

Normal padrão Probit Φ

Valor mínimo extremo Complementário log-log log log 1

Valor máximo extremo Log-log log log

Cauchy Cauchit tan 0.5

Ramalho, Ramalho & Murteira (2010) apresento sugestões para diferentes funções de ligação. A função de Log verossimilhança ℓ , | pode ser calculada como

ℓ , | ℓ , | ,

onde ℓ , | log Γ log Γ log Γ 1

(2.6)

5  

1 log 1 1 log 1 Note que é uma função de , o vetor de parâmetros de regressão. O vetor de parâmetros , é calculado por Máxima Verossimilhança (ML), usando o logit como função de ligação temos

logit log1

exp

1 exp

∗ logit log

1

obtendo as derivativas temos

ℓ

ℓ

onde,

ℓ

log Γ

log Γ 1 log log 1

1 log log 1

∗ ∗

1′

1logit′

1

sendo ∙ a função digama1, ∗ 1 e ∗ definido em (2.9), pelo que a derivada fica

ℓ

∗ ∗ 1

que pode ser escrito de maneira matricial como ∗ ∗ onde X é a matriz de desenho e T é uma matriz diagonal que contém os elementos 1 . Do mesmo jeito pode ser calculada a derivada

(2.7)

(2.8)

(2.9)

(2.10)

6  

ℓℓ

ℓ ∗ ∗ log 1 1

Os estimadores de máxima verossimilhança são obtidos pela solução do seguinte sistema

ℓ

0

ℓ0

A solução deste sistema não possui uma forma fechada, fazendo-se necessário o uso de algoritmos de otimização não-linear, como o algoritmo quasee-Newton BFGS. Para detalhes, ver Press et al. (1992). O algoritmo Newton-Raphson precisa da matriz de Informação Observada , sendo o negativo da matriz Hessiana ℓ dada por

ℓ

ℓ ℓ

ℓ

ℓ

e seus elementos são

ℓ

ℓ

ℓ

ℓ

ℓ

ℓℓ

onde calculando as derivadas parciais usando (2.10) e (2.11) tem se que

ℓ1 ∗ ∗ 1 2

ℓ

∗ ∗ 1

ℓ1 1 ′

(2.11)

1 Em geral, a função poligama está definida para m = 0, 1, . . . , como ⁄ log Γ , 0. Para maior detalhe ver Dishon & Weiss (1980). 

7  

onde 1 e 1 . Aplicando valor esperado se obtém a matriz de Informação Esperada de Fisher dada por

e seus elementos são

1

1

1 1 ′

elementos que podem ser expressados de forma matricial como

tr

onde diag , … , , diag , … , , , … , e diag , … , com elementos

1 , ,

1 ,

1 1 .

As estimativas de Máxima Verossimilhança para podem ser calculadas iterativamente pelos algoritmos:

Newton-Raphson Escoring de Fisher

ℓ ℓ

até alcançar convergência, isto é

, 10 ,

(2.12)

8  

onde sugestões para os valores iniciais dos parâmetros foram apresentados por Ferrari & Cribari-Neto (2004). Para amostras grandes, e sob condições de regularidade, os estimadores de Máxima Verossimilhança e , têm aproximadamente distribuição de densidade conjunta normal (p+1) multivariada, i.e.,

~ ,

onde os erros padrões dos estimadores podem ser obtidos da matriz de variâncias e covariâncias

como

diag ⁄ . Uma extensão do modelo de regressão Beta simples foi apresentado formalmente por Simas et al. (2010) é o modelo de regressão Beta com dispersão variável, i.e., considerando um parâmetro de precisão não fixo para todas as variáveis, mas que pode ser modelado de um jeito parecido à média por covariaveis. Mais especificamente ~ , para cada i independentemente, com funções de ligação:

sendo e vetores dos coeficientes de regressão, e preditores lineares e e vetores de regressores. A função de ligação mais comum para é log, forçando que seja sempre positivo. As estimativas de ML são feitas do mesmo jeito que no caso do modelo simples, mas só se tem que substituir por na equação (2.6).

(2.13)

9  

3. IMPLEMENTAÇÃO E APLICAÇÕES

3.1. Modelo de Regressão Beta Simples

Com os dados Gasoline Yield Data tomados de Prater (1956) será aplicado um modelo básico de regressão Beta sugerido por Ferrari e Cribari-Neto (2004). tem-se 32 observações onde a variável de interesse yield, é a proporção de petróleo bruto convertida em gasolina depois dos processos de destilação e fracionamento, onde y pode ser estimado naturalmente pelo modelo de regressão Beta. As covariáveis são dois: temp é temperatura (em graus Fahrenheit) a qual toda a gasolina é evaporada, e batch que é um fator que contém dez níveis correspondentes a diferentes lotes, os quais foram sometidos a diferentes condições experimentais. Veja-se o gráfico 2. Inicialmente proponhamos um modelo simples com função de ligação logit para e parâmetro de precisão constante. > data("GasolineYield", package = "betareg")

> BR_logit = betareg(yield ~ batch + temp, data = GasolineYield)

> summary(BR_logit)

Call: betareg(formula = yield ~ batch + temp, data = GasolineYield) Standardized weighted residuals 2: Min 1Q Median 3Q Max -2.8750 -0.8149 0.1601 0.8384 2.0483 Coefficients (mean model with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -6.1595710 0.1823247 -33.784 < 2e-16 *** batch1 1.7277289 0.1012294 17.067 < 2e-16 *** batch2 1.3225969 0.1179020 11.218 < 2e-16 *** batch3 1.5723099 0.1161045 13.542 < 2e-16 *** batch4 1.0597141 0.1023598 10.353 < 2e-16 *** batch5 1.1337518 0.1035232 10.952 < 2e-16 *** batch6 1.0401618 0.1060365 9.809 < 2e-16 *** batch7 0.5436922 0.1091275 4.982 6.29e-07 *** batch8 0.4959007 0.1089257 4.553 5.30e-06 *** batch9 0.3857930 0.1185933 3.253 0.00114 ** temp 0.0109669 0.0004126 26.577 < 2e-16 *** Phi coefficients (precision model with identity link): Estimate Std. Error z value Pr(>|z|) (phi) 440.3 110.0 4.002 6.29e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Type of estimator: ML (maximum likelihood) Log-likelihood: 84.8 on 12 Df Pseudo R-squared: 0.9617 Number of iterations: 51 (BFGS) + 3 (Fisher scoring)

10  

Observe-se que a covariável temp e o fator batch são significativos para explicar a variável de interesse yield, além disso, obtém-se um valor alto da precisão de 440,3, e o seu valor é determinante para um bom ajuste. 3.2. Modelo de Regressão Beta com precisão variável Embora o modelo de regressão Beta já incorpora informação da dispersão dos dados de um jeito natural (ver equação 2.4b), considerando o parâmetro de precisão como não fixo, podemos modelar a heterocedasticidade, isto utilizando uma função de ligação log e a variável temp como covariável como na equação (2.13). > BR_logitV = betareg(yield ~ batch + temp|temp, data = GasolineYield)

> summary(BR_logitV)

Call: betareg(formula = yield ~ batch + temp | temp, data = GasolineYield, link = "logit") Standardized weighted residuals 2: Min 1Q Median 3Q Max -2.5399 -0.7792 -0.1167 0.8621 2.9419 Coefficients (mean model with logit link):

Gráfico 2. Temperatura em graus Fahrenheit onde a gasolina é evaporada Vs. proporção de petróleo bruto convertida em gasolina depois dos processos de destilação e fracionamento. As observações ligadas por números iguais representam um nível do fator batch.

11  

Estimate Std. Error z value Pr(>|z|) (Intercept) -5.9232361 0.1835262 -32.275 < 2e-16 *** batch1 1.6019877 0.0638561 25.087 < 2e-16 *** batch2 1.2972663 0.0991001 13.090 < 2e-16 *** batch3 1.5653383 0.0997392 15.694 < 2e-16 *** batch4 1.0300720 0.0632882 16.276 < 2e-16 *** batch5 1.1541630 0.0656427 17.582 < 2e-16 *** batch6 1.0194446 0.0663510 15.364 < 2e-16 *** batch7 0.6222591 0.0656325 9.481 < 2e-16 *** batch8 0.5645830 0.0601846 9.381 < 2e-16 *** batch9 0.3594390 0.0671406 5.354 8.63e-08 *** temp 0.0103595 0.0004362 23.751 < 2e-16 *** Phi coefficients (precision model with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 1.364089 1.225781 1.113 0.266 temp 0.014570 0.003618 4.027 5.65e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Type of estimator: ML (maximum likelihood) Log-likelihood: 86.98 on 13 Df Pseudo R-squared: 0.9519 Number of iterations: 33 (BFGS) + 28 (Fisher scoring)

e com isto se tem evidencia de melhora ao incluir a temperatura de evaporação da gasolina como covariável do parâmetro de precisão . Agora os dois modelos (com precisão constante e não) são comparados pelo Teste de Ração de Verossimilitudes. > lrtest(BR_logit, BR_logitV)

Likelihood ratio test Model 1: yield ~ batch + temp Model 2: yield ~ batch + temp | temp #Df LogLik Df Chisq Pr(>Chisq) 1 12 84.798 2 13 86.977 1 4.359 0.03681 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

sendo evidencia de que o Modelo 2 (com precisão não constante) é melhor do que o Modelo 1. Note-se que é o mesmo que testar uma hipótese nula de igual dispersão vs. dispersão variável, i.e.,

: ~ , : ~ , cuja hipóteses nula é rejeitada a favor da hipóteses alternativa. O critério AIC também comprovou que o modelo 2 é melhor. > AIC(BR_logit,BR_logitV) df AIC BR_logit 12 -145.5951 BR_logitV 13 -147.9541

12  

3.3. Propostas de funções de ligação Ramalho & Murteira (2010) fala dos critérios de seleção de diferentes funções de ligação para o modelo de regressão Beta, onde uma simples escolha de uma delas pode melhorar significativamente o ajuste, especialmente quando se tem dados pertos aos valores extremos zero e um. Com o objetivo de comparar o ajuste, utilizaremos três modelos com igual precisão

e diferentes funções de ligação: logit, probit e loglog.

> BR_logit = betareg(yield ~ batch + temp, GasolineYield,link="logit")

> BR_probit = betareg(yield ~ batch + temp, GasolineYield,link="probit")

> BR_loglog = betareg(yield ~ batch + temp, GasolineYield,link="loglog")

> AIC(BR_logit, BR_probit, BR_loglog)

df AIC BR_logit 12 -145.5951 BR_probit 12 -155.6575 BR_loglog 12 -168.3101

Com o modelo com função de ligação loglog proposto também em Cribari-Neto and Lima (2007), é obtido o menor valor do critério, sendo não só melhor do que o modelo logit com

Gráfico 3. Curvas de predição da proporção de petróleo bruto convertida em gasolina depois dos processos dedestilação e fracionamento dada a temperatura em graus Fahrenheit onde a gasolina é evaporada. Cada curva comdiferente cor representa um nível do fator batch.

13  

precisão fixa, mas ainda melhor do que o modelo com modelação da dispersão proposto em linhas anteriores. Olhe-se o gráfico 4.

4. MODELOS DE REGRESSÃO BETA MAIS GERAIS

Varying dispersion beta regression models: Smithson & Verkuilen (2006). A general class of beta regression models: Simas, Barreto-Souza & Rocha (2010). Inflated beta regression models: Cook, Kieschnick, McCullough (2008), Ospina & Ferrari

(2010,2012a), Calabrese (2012). Truncated inflated beta regression models: Pereira, Botter & Sandoval (2011, 2013). Semi-parametric beta regression: Branscum, Jonhson & Thurmond (2007), Weihua et al

(2012). Time series: Rydlewski (2007), Rocha & Cribari–Neto (2009), Billio & Casarin (2011),

Casarin, Dalla Valle, Leisen (2012); da-Silva, Migon & Correia (2011), da-Silva & Migon (2012), Guolo & Varin (2012).

Multivariate beta regression: Souza & Moura (2012a, 2012b)

Gráfico 4. Curvas médias de predição da proporção de petróleo bruto convertida em gasolina depois dos processosde destilação e fracionamento dada a temperatura em graus Fahrenheit onde a gasolina é evaporada. Cada curvarepresenta a função média ajustando os dados com modelos com funções de ligação logit, probit e loglog.

14  

Mixed beta regression: Zimprich (2010), Verkuilen & Smithson (2012), Figueroa–Zúñiga, Arellano Valle & Ferrari (2013), Bonat, Ribeiro Jr & Zeviani (2013).

Errors-in-variables beta regression models: Carrasco, Ferrari, Arellano–Valle (2012) (more

later). Beta rectangular regression models: Bayes, Bazán & García (2012).

5. PESQUISAS

Além, são apresentados diferentes pesquisas e aplicações do Modelo Beta: Johnson et al. (1995, p. 235). “The beta distributions are among the most frequently

employed to model theoretical distributions”.

Bury (1999). “Applications of the beta distribution in engineering”.

Janardan and Padmanabhan (1986). “Modelling of hydrological variables using the beta distribution”.

McNally (1990). “Use of the beta distribution in the study of reproducibility of cows”.

Graham e Hollands (1990) e Milyutin e Yaromenko (1991). “Use the beta distribution in studies of indices related to the transmission of solar radiation”

Maffet and Wackerman (1991). “Power of radar signals is modeled by using the beta law.

Wiley et al. (1989). “Develop a beta model to estimate the probability of HIV transmission during sexual intercourse involving infected and non-infected individuals”.

15  

ANEXOS

Implementação do Scoring de Fisher para o Modelo de Regressão Beta Simples. library(betareg) library(lmtest) y = GasolineYield$yield X = model.matrix(gy_logit) EFBeta = function(x,y,th) { #Contas iniciais b = th[1:p] phi = exp(th[p+1]) eta = X%*%b y. = log(y/(1-y)) mu = exp(eta)/(1+exp(eta)) mu. = digamma(mu*phi)-digamma((1-mu)*phi) #Contas Adicionais t = as.numeric(mu*(1-mu)) om = psigamma(mu*phi,2)+psigamma((1-mu)*phi,2) c = phi*(mu*om-psigamma((1-mu)*phi,2)) w = as.numeric(phi*om*t^2) d = as.numeric(((1-mu)^2)*psigamma((1mu)*phi,2)+(mu^2) *psigamma(mu*phi,2)-psigamma(phi,2)) TT = diag(t) W = diag(w) D = diag(d) #Gradiente G_bet = phi*(t(X)%*%TT%*%(y.-mu.)) G_phi = sum(mu*(y.-mu.)+log(1-y)-digamma((1-mu)*phi)+digamma(phi)) GG = rbind(G_bet,G_phi) #Matriz Esperada de Fisher M_bb = phi*(t(X)%*%W%*%X) M_bp = t(X)%*%TT%*%c M_pb = t(M_bp) M_pp = sum(diag(D)) MIEF = cbind(rbind(M_bb,M_pb),rbind(M_bp,M_pp)) return(list(GG = GG, MIEF = MIEF)) } #Função de Log-Verossimilhança logverB<-function(param,y,X) { p = ncol(X) eta = X%*%param[1:p] mu = exp(eta)/(1+exp(eta)) vero = sum(lgamma(param[p+1])-lgamma(mu*param[p+1])-lgamma((1-mu) *param[p+1])+((mu*param[p+1])-1)*log(y)+(((1-mu)*param[p+1])-1) *log(1-y)) return(-vero) }

16  

#Regressão Beta RegBeta = function(X,y,b,phi,MaxIter=200,error=0.0001) { n = nrow(X) p = ncol(X) thv = c(b,phi) count = 0 criterio = 1 while(criterio > error) { count = count + 1 GM = EFBeta(x,y,thv) ite = solve(GM$MIEF)%*%GM$GG thn = thv + ite criterio = sqrt(t(thv-thn)%*%(thv-thn)) thv = thn npar = p+1 if (count==MaxIter) { break } } EPbeta = sqrt(diag(solve(GM$MIEF))) loglik = logverB(param,y,X) #Criterios AIC<- -2*loglik +2*npar BIC <- -2*loglik +log(n)*npar HQ <- -2*loglik +2*log(log(n))*npar return(list(iter=count, theta = thn, EPbeta = EPbeta ,AIC = AIC, BIC = BIC, HQ = HQ)) }

17  

REFERÊNCIAS

1. Bayer, F.M. (2011) Modelagem e Inferência em Regressão Beta, Tese de Doutorado,

Universidade Federal de Pernambuco, Recife, Brasil. 2. Biguelini, C.B.,(2009). Modelo de Regressão Beta para a Análise da Origem dos Problemas

de Sistemas Prediais. Monografia. Universidade Federal do Rio Grande do Sul. 3. Branscum, A.J., Johnson, W.O. & Thurmond, M.C. (2007). Bayesian beta regression:

applications to household expenditure data and genetic distance between foot-and-mouth disease viruses. Australian and New Zealand Journal of Statistics, 49, 287–301.

4. Cribari–Neto, F. & Zeiles, A. (2010). Beta regression in R. Journal of Statistical Software,

34, 1–24. 5. Espinheira, P.L., Ferrari, S.L.P. & Cribari–Neto, F. (2008a). Influence diagnostics in beta

regression. Computational Statistics and Data Analysis, 52, 4417–4431. 6. Ferrari, S.L.P. & Cribari–Neto, F. (2004). Beta regression for modelling rates and

proportions. Journal of Applied Statistics, 31, 799–815. 7. Ferrari, S.L.P., Espinheira, P.L. & Cribari–Neto, F. (2011). Diagnostic tools in beta

regression with varying dispersion. Statistica Sinica, 65, 337–351. 8. Ospina, R., Cribari–Neto, F. & Vasconcellos, K.L.P. (2006). Improved point and interval

estimation for a beta regression model. Computational Statistics and Data Analysis, 51, 960–981. Erratum at Computational Statistics and Data Analysis, 55, 2445.