A análise de dados binários é usualmente feita através da regressão logística, mas esse modelo possui limitações. Modificar a função de ligação da regressão permite maior flexibilidade na modelagem e diversas propostas já foram feitas nessa área.

Fenômenos que assumem apenas dois estados são chamados de binários e seu estudo sempre foi de grande interesse. Seja prever o resultado de uma eleição ou estudar a prevalência do uso de drogas ou remédios em uma certa população, é comum a presença de dados dicotômicos e o problema de estimar a proporção de ocorrência de um dos eventos faz-se importante. Além disso, pode ser interessante avaliar quanto algumas variáveis influenciam na proporção desse evento. No entanto, técnicas de regressão linear não podem ser utilizadas devido a violação de algumas suposições.

Diversas soluções para esse problema foram propostas e a regressão logística, um caso particular dos modelos lineares generalizados (Nelder e Wedderburn, 1972), é uma das mais utilizadas. Nesses, uma função de ligação é aplicada na combinação linear dos preditores de modo que a resposta seja um valor apropriado para o problema. Para dados binários, é natural utilizar funções quantil (inversa da função de distribuição) como é o caso da regressão logística que utiliza a distribuição logística como o nome sugere. Outras distribuições como a normal padrão (probito) e valor extremo (complementar log-log) também são opções comuns para a função de ligação (Agresti, 2002).

Diversas outras funções de ligações foram propostas na literatura para estender as mencionadas. Prentice (1975) apresentou uma função de ligação baseada no logaritmo da distribuição F-Snedecor. Aranda-Ordaz (1981) propôs uma função assimétrica que tem os modelos logístico e complementar log-log como casos particulares. Stukel (1988) introduziu uma generalização do modelo logístico adicionando dois parâmetros à função de ligação logística que controlam o comportamento das caudas da distribuição. Caron et al. (2009) e Caron (2010) sugerem um modelo baseado na distribuição Weibull que possui o modelo complementar log-log como caso especial e aproxima os modelos logístico e probito.

Apesar das qualidades destes modelos, até onde se sabe, não existe função para a estimação desses nos pacotes estatísticos a não ser pelos três primeiros casos mencionados. Todas as outras funções compartilham a qualidade de serem paramétricas, ou seja, há elementos extras a serem estimados além dos coeficientes de regressão. Prentice (1976) aponta dificuldades em estimar seu modelo com os dois parâmetros livres utilizando máxima verossimilhança e sugere fixar um deles e calcular o outro.

Dados

Banco de dados sobre mortalidade de besouros (Bliss, 1935). O objetivo original desse estudo era obter um inseticida eficaz contra besouros. Para isso, 481 animais foram expostos a diferentes concentrações de dissulfeto de carbono (CS2) durante cinco horas e contou-se o número de insetos mortos.

Esse conjunto de dados é conhecido por não ser bem ajustado por modelos simétricos, em particular logito e probito; por conta disso, é amplamente citado em trabalhos que buscam alternativas a esses modelos.

library(VGAM)
data("flourbeetle")
show(flourbeetle)
##    logdose CS2mgL exposed killed
## 1 1.690728  49.06      59      6
## 2 1.724194  52.99      60     13
## 3 1.755189  56.91      62     18
## 4 1.784189  60.84      56     28
## 5 1.811307  64.76      63     52
## 6 1.836894  68.69      59     53
## 7 1.860996  72.61      62     61
## 8 1.883888  76.54      60     60
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
grid()

Modelos 

m <- lapply(c("logit", "probit", "cloglog", "cauchit"), function(type)
    glm(cbind(killed,exposed-killed) ~ logdose, data = flourbeetle, family = binomial(link = type)))
names(m) <- c("logit", "probit", "cloglog", "cauchit")
m
## $logit
## 
## Call:  glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type), 
##     data = flourbeetle)
## 
## Coefficients:
## (Intercept)      logdose  
##      -60.72        34.27  
## 
## Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
## Null Deviance:       284.2 
## Residual Deviance: 11.22     AIC: 41.42
## 
## $probit
## 
## Call:  glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type), 
##     data = flourbeetle)
## 
## Coefficients:
## (Intercept)      logdose  
##      -34.94        19.73  
## 
## Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
## Null Deviance:       284.2 
## Residual Deviance: 10.11     AIC: 40.3
## 
## $cloglog
## 
## Call:  glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type), 
##     data = flourbeetle)
## 
## Coefficients:
## (Intercept)      logdose  
##      -39.57        22.04  
## 
## Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
## Null Deviance:       284.2 
## Residual Deviance: 3.44  AIC: 33.64
## 
## $cauchit
## 
## Call:  glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type), 
##     data = flourbeetle)
## 
## Coefficients:
## (Intercept)      logdose  
##      -77.32        43.52  
## 
## Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
## Null Deviance:       284.2 
## Residual Deviance: 20.15     AIC: 50.35

Observemos os resultados de cada ajuste.

plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
lines(fitted(m[[2]]) ~ logdose, data = flourbeetle, col = "green")
lines(fitted(m[[3]]) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[4]]) ~ logdose, data = flourbeetle, col = "blue")
grid()

Percebemos que algumas curvas subestimam algumas das probabilidades observadas de mortalidade e outras superestimam as probabilidades. Por exemplo, na logdose 1.724194 todas as curvas estimadas subestimam a probabilidade de mortalidade observada, sendo as que pior estimam essa probabilidade as obtidas com as ligações probit, logit e cauchit nessa ordem. Logo percebemos que a probabilidade estimada na logdose 1.784189 todas as curvas estimadas superestimam essa probabilidade observada.

Uma melhoria nas funções de ligação na regressão logística seria incorporar parámetros a estas funções com o intuito claro de diminuir os erros de estimação.

Funções de ligação paramétricas

Modelos lineares generalizados (GLM’s) permitem o tratamento de problemas de regressão nos quais a resposta pode ser distribuída de forma não normal. Mais especificamente, a distribuição da resposta pode ser qualquer distribuição numa família exponencial. Isso inclui gama, Poisson, binomial, normal e respostas gaussianas inversas.

Além disso, uma função de ligação conectando a resposta média ao preditor linear deve ser escolhida. GLM’s com ligações canônicas, como a ligação logit na regressão binomial, garantem o máximo de informação e uma interpretação simples dos parâmetros de regressão porque neste caso obtemos um modelo linear para o parâmetro natural da família exponencial subjacente.

Por exemplo, a ligação logit fornece uma representação simples das probabilidades que auxilia na interpretação dos resultados. Além disso, a concavidade da função de verossimilhança garante a unicidade dos estimadores de máxima verossimilhança (MLE).

No entanto, as ligações canônicas nem sempre fornecem o melhor ajuste disponível para um determinado conjunto de dados. Neste caso, a ligação pode ser mal especificada, o que pode produzir um viés substancial para as estimativas dos parâmetros de regressão, bem como para as estimativas da resposta média, por exemplo, ver Czado e Santner (1992) no caso da resposta binomial.

A abordagem mais comum para se proteger contra tal especificação incorreta é incorporar a função de ligação em uma ampla classe paramétrica de ligações \(\mathcal{F}=\{F(\cdot,\psi), \psi\in\Psi\}\), que inclui as ligações canônicas como um caso especial quando \(\psi = \psi_*\) digamos.

Muitas dessas classes paramétricas de funções de ligação para dados de regressão binária foram propostas na literatura. Por exemplo, Van Montford e Otten (1976), Copenhaver e Mielke (1977), Aranda-Ordaz (1981), Guerrero e Johnson (1982), Morgan (1983) e Whittmore (1983) propuseram famílias de um parâmetro enquanto Prentice (1976), Pregibon (1980), Stukel (1988) e Czado (1992) consideraram duas famílias de parâmetros. O caso geral de funções de ligação em GLMs foi estudado por Pregibon (1980) e Czado (1992). Czado (1995) desenvolveu critérios sobre como escolher tal família.

Existem muitos casos em que obtemos uma especificação incorreta da função de ligação. O motivo é simples: temos que escolher a função do ligação antes de termos informações suficientes sobre a escolha da ligação. Assim, gostaríamos de descrever uma maneira de melhorar a qualidade do ajuste dos modelos lineares generalizados (GLM’s), reduzindo o desvio.

Transformações de potência gerias \(h(\cdot)\)

A seguir, vamos apresentar as funções de ligação transformada de potência \(h(\cdot)\) (Czado, 2007) para valores específicos de \(\psi\). Portanto, temos que passar à função um valor inicial \(\eta_0\in\mathbb{R}\) e \(\psi_1\) ou \(\psi_2\in \mathbb{R}\) para uma modificação de cauda única ou os valores do vetor \(\psi\in\mathbb{R}^2\) para uma modificação de ambas as caudas.

Para uma modificação da cauda esquerda, todo ponto \(<\eta_0\) será modificado, enquanto para uma modificação da cauda direita todo ponto é \(\geq\eta_0\). Uma modificação de ambas as caudas modifica ambas as caudas, ou seja, todos os pontos \(<\eta_0\) e \(\geq \eta_0\) são modificados.

Ao definirmos os parâmetros como 1, não obtemos nenhuma modificação, ou seja, uma linha reta. Se definirmos um parâmetro como 1 na modificação de ambas as caudas, obtemos uma única modificação da cauda, por exemplo, \(\psi= (1,\psi_2)\) modifica a cauda esquerda. Numa modificação de cauda direita, o parâmetro \(\psi_1 <1\) diminui a inclinação, enquanto a configuração \(\psi_1> 1\) a aumenta. Na modificação da cauda esquerda ocorre o contrário para \(\psi_2\).

Todas as funções de ligação paramétricas que apresentamo a seguir estão programadas no pacote glmx, o qual depende do pacote gld.

library(glmx); library(gld)

Pregibon

A família de funções de liagação proposta por Pregibon (1980) depende de dois parámetros e é definida como \[ g(\mu;\alpha,beta) = \dfrac{\mu^{\alpha-\beta}-1}{\alpha-\beta}-\dfrac{(1-\mu)^{\alpha+\beta}-1}{\alpha+\beta}\cdot \]

Quando \(\alpha=\beta=0\) obtemos a ligação logit. Para \(\beta=0\) esta função é simétrica e \(\beta\) controla a assimetria, o peso das caudas é controlado por \(\alpha\). A implementação usa a distribuição lambda generalizada.

## Pregibon model
pregibin <- function(shape) binomial(link = pregibon(shape[1], shape[2]))
ajuste1 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = pregibin, xstart = c(0, 0), xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste1) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Pregibon"), col = c("black","red"), fill = c(1,2))
grid()

Gosset

A família Gosset é dada pelo inverso da função de distribuição acumulada da distribuição \(t\)-Student. Os graus de liberdade \(\nu\) controlam o peso das caudas e está restrito a valores positivos. Para \(\nu=1\) se obtém a ligação Cauchy (cauchit) e para \(\nu\to\infty\) esta função de ligação converga para a probit.

## Gosset model
gossbin <- function(nu) binomial(link = gosset(nu))
ajuste2 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = gossbin, xstart = 0, xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste2) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Gosset"), col = c("black","red"), fill = c(1,2))
grid()

Aranda-Ordaz

A função Aranda-Ordaz (1981) simétrica é definida como \[ g(\mu;\phi) = \dfrac{2}{\phi}\dfrac{\mu^\phi-(1-\mu)^\phi}{\mu^\phi+(1-\mu)^\phi} \] e está implementada na função ao1. A função Aranda-Ordaz (1981) assimétrica é dada por \[ g(\mu;\phi) = \log\left( \dfrac{(1-\mu)^{-\phi}-1}{\phi} \right)\cdot \] Ambas contêm a função logit para \(\phi=0\) e \(\phi=1\) respectivamente, onde esta última também inclui o complementar log-log (cloglog) para \(\phi=0\).

## Aranda-Ordaz model asymmetric
aranda.ordaz <- function(phi) binomial(link = ao2(phi))
ajuste3 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = aranda.ordaz, xstart = 0, xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste3) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Aranda-Ordaz"), col = c("black","red"), fill = c(1,2))
grid()

Guerrero-Johnson

A seguinte função de ligação foi proposta em Guerrero, V. and Johnson, R. (1982). “Use of the Box-Cox Transformation with Binary Response Models.” Biometrika, 69, 309–314. A família Guerrero-Johnson é definida como \[ g(\mu;\phi) = \frac{1}{\phi}\left( \Big(\dfrac{\mu}{1-\mu}\Big)^\phi-1\right), \] sendo simétrica e contendo a função logística quando \(\phi=0\).

## Guerrero-Johnson model
guerrero.johnson = function(phi) binomial(link = gj(phi))
ajuste4 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = guerrero.johnson, xstart = 0, xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste3) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Guerrero-Johnson"), col = c("black","red"), fill = c(1,2))
grid()

Stukel

A família de funções de ligação Stukel, proposta por Thérèse A. Stukel em 1988 também generaliza a ligação logito, mas preservando a propriedade de que se \(\eta=0\) então \(\mu=0.5\). Define-se como \[ \log\left(\dfrac{\mu}{1-\mu}\right) = h_\alpha(\eta), \] onde \(\alpha=(\alpha_1,\alpha_2)\) é o vetor de parâmetros que definem a forma da função de ligação. Desta forma, o modelo proposto é \[ \mu_\alpha(\eta)=\dfrac{\exp\big(h_\alpha(\eta)\big)}{1+\exp\big(h_\alpha(\eta)\big)}\cdot \]

A função \(h_\alpha(\eta)\) é definida como \[ h_\alpha(\eta)=\left\{ \begin{array}{ccl} \alpha_1^{-1}\Big(\exp(\alpha_1|\eta|)-1\Big) & \mbox{se} & \alpha_1>0 \\ \eta & \mbox{se} & \alpha_1=0 \\ -\alpha_1^{-1}\log\big(1-\alpha_1|\eta|\big) & \mbox{se} & \alpha_1<0 \end{array}\right. \] se \(\eta\geq 0\) ou, equivalentemente, se \(\mu\geq 0.5\) e \[ h_\alpha(\eta)=\left\{ \begin{array}{ccl} -\alpha_2^{-1}\Big(\exp(\alpha_2|\eta|)-1\Big) & \mbox{se} & \alpha_2>0 \\ \eta & \mbox{se} & \alpha_2=0 \\ \alpha_2^{-1}\log\big(1-\alpha_2|\eta|\big) & \mbox{se} & \alpha_2<0 \end{array}\right. \] se \(\eta\leq 0\) ou, equivalentemente, se \(\mu\leq 0.5\).

Quando \(\alpha_1 = \alpha_2 = 0\) obtemos o modelo logístico (logit). Caso contrário cada parâmetro governa o comportamento da curva de maneira diferente. Se \(\alpha_1 = \αlpha_2\) a curva obtida é simétrica, se \(\alpha_1\neq \alpha_2\) a função de ligação é assimetrica. O modelo probito é obtido quando \(\alpha_1 = \αlpha_2 ≈ 0.165\), quando \(\alpha_1 = −0.037\) e \(\alpha_2 = 0.62\) se obtém a chamada função de ligação log-log e quando \(\alpha_1 = 0.62\) e \(\alpha_2 = −0.037\) obtem-se a função de ligação complemetar log-log. No caso de \(\alpha_1 = \alpha_2 ≈ −0.077\) a função de ligação chama-se Laplace.

Podemos utilizar também uma aproximação à função de ligação paramétrica Stukel, também proposta por Stukel (1988). Suponhamos ajustamos um modelo com a função de ligação \(g_1(\mu)=g(\mu;\alpha_1,\delta_1)\) quando a função de ligação correta, mas desconhecida, é \(g_0(\mu)=g(\mu;\alpha_1,\delta_1)\).

Utilizando a expansão em série de Taylor de primeira órdem, no ponto \(g_1(\mu)\), temos a seguinte relação aproximada \[ g_0(\mu) \approx g_1(\mu)+(\alpha_0-\alpha_1)D_\alpha(g_1)+(\delta_0-\delta_1) D_\delta(g_1), \] onde \[ D_\alpha(g_1)=\left. \dfrac{\partial}{\partial\alpha} g(\mu;\alpha,\delta)\right|_{\alpha=\alpha_1,\delta=\delta_1} \] e similarmente para \(D_\delta(g_1)\).

Observemos que agora podemos aproximar a função de ligação correta \(g_0(\mu) = X\beta\) por \[ g_1(\mu)=X\beta+Z\gamma, \] onde \(Z=\big(D_\alpha(g_1),D_\delta(g_1)\big)\) e \(\gamma=-(\alpha_0-\alpha_1,\delta_0-\delta_1)\).

A forma de aproximar as estimativas dos parâmetros do modelo ampliado é escrevendo como preditor linear \[ \eta = \beta_0+\beta_1 \mbox{logdose}+\dfrac{1}{2}\alpha_1\widehat{\eta}^2\mbox{I}_{\widehat{\eta}>0}-\dfrac{1}{2}\alpha_2\widehat{\eta}^2\mbox{I}_{\widehat{\eta}<0}, \] para o caso do conjunto de dados considerado.

Neste modelo ampliado as estimativas dos coeficientes \(\alpha_1\) e \(\alpha_2\) representam aproximações de primeira ordem das estimativas destes mesmos parâmetros se utilizados diretamente na definição do modelo de regressão. Este modelo permite além de encontrar aproximações dos parâmetros que definem a nova função de ligação permite também testar se devemos ou não utilizar funções de ligação paramétricas

ajuste5 = update(m[[1]], formula = . ~ . + ifelse(predict(m[[1]])>0, 0.5*predict(m[[1]])^2,0) + 
                   ifelse(predict(m[[1]])<0, -0.5*predict(m[[1]])^2,0), family = binomial(link="logit"))
summary(ajuste5)
## 
## Call:
## glm(formula = cbind(killed, exposed - killed) ~ logdose + ifelse(predict(m[[1]]) > 
##     0, 0.5 * predict(m[[1]])^2, 0) + ifelse(predict(m[[1]]) < 
##     0, -0.5 * predict(m[[1]])^2, 0), family = binomial(link = "logit"), 
##     data = flourbeetle)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
## -0.26801   0.67253  -0.44440  -0.40939   1.02043  -0.85234   0.03192   0.61935  
## 
## Coefficients:
##                                                          Estimate Std. Error
## (Intercept)                                              -53.7809    16.6477
## logdose                                                   30.1860     9.3888
## ifelse(predict(m[[1]]) > 0, 0.5 * predict(m[[1]])^2, 0)    0.3597     0.2666
## ifelse(predict(m[[1]]) < 0, -0.5 * predict(m[[1]])^2, 0)  -0.1765     0.2534
##                                                          z value Pr(>|z|)   
## (Intercept)                                               -3.231  0.00124 **
## logdose                                                    3.215  0.00130 **
## ifelse(predict(m[[1]]) > 0, 0.5 * predict(m[[1]])^2, 0)    1.349  0.17734   
## ifelse(predict(m[[1]]) < 0, -0.5 * predict(m[[1]])^2, 0)  -0.696  0.48613   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.2024  on 7  degrees of freedom
## Residual deviance:   3.0416  on 4  degrees of freedom
## AIC: 37.24
## 
## Number of Fisher Scoring iterations: 4
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste5) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Ligação Stukel"), col = c("black","red"), fill = c(1,2))
grid()

Percebemos claramente a melhoria na qualidade na estimação das probabilidades de mortalidade.

Selecionando a função de ligação

Um problema comdiversas destas funções de ligação é a falta de um critério para escolher qual à mais adequada. Estamos dizendo que na maioria destes ajustes o AIC, AICc e BIC não estão disponíveis nos objeto de ajuste, nem está claro se estes critérios devem ser utilizados nestas situações. Devemos lembrar que estamos adicionando parámetros aos modelos, os quais devem ser considerados no cálculo desses critérios.

Como alternativa propomos utilizar diversas medidas de precisão, sendo estas:

Um erro de estimação é a diferença entre um valor observado e seu estimado. Pode ser escrito como \[ e_i = y_i-\widehat{y}_i \] sendo \(y_i\) o valor observado da resposta e \(\widehat{y}_i\) sua estimativa pelo modelo.

Observe que os erros de estimação são diferentes dos resíduos de duas maneiras. Podemos medir a precisão da estimação utilizando os erros dependentes de escala.

Os erros de estimação estão na mesma escala dos dados. As medidas de precisão baseadas apenas em \(e_i\) são, portanto, dependentes da escala e não podem ser usadas para fazer comparações entre modelos que envolvem unidades diferentes.

As duas medidas dependentes de escala mais comumente usadas são baseadas em erros absolutos ou erros quadrados:

Erro médio absoluto: MAE

\[ MAE = \dfrac{1}{n}\sum_{i=1}^n |y_i-\widehat{y}_i| \]

Erro absoluto médio: RMSE

\[ RMSE = \sqrt{\dfrac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2} \]

Ao comparar modelos o MAE é popular porque é fácil de entender e calcular. Um método de precisão que minimize o MAE levará a precisões medianas, enquanto a minimização do RMSE levará a previsões médias. Consequentemente, o RMSE também é amplamente utilizado, apesar de ser mais difícil de interpretar.

MAE = function(obs,adj) mean(abs((obs-adj)))
RMSE = function(obs,adj) sqrt(mean((obs-adj)^2))
obs = flourbeetle$killed/flourbeetle$exposed
erro1 = MAE(obs,ajuste1$fitted.values)
erro2 = MAE(obs,ajuste2$fitted.values)
erro3 = MAE(obs,ajuste3$fitted.values)
erro4 = MAE(obs,ajuste4$fitted.values)
erro5 = MAE(obs,ajuste5$fitted.values)
erro6 = RMSE(obs,ajuste1$fitted.values)
erro7 = RMSE(obs,ajuste2$fitted.values)
erro8 = RMSE(obs,ajuste3$fitted.values)
erro9 = RMSE(obs,ajuste4$fitted.values)
erro0 = RMSE(obs,ajuste5$fitted.values)
## adicionando espaço extra à margem direita do gráfico dentro do quadro
par(mar=c(4, 6, 1, 1) + 0.1, cex.axis=0.6)
## Plotando o primeiro conjunto de dados e desenhando seu eixo
plot(1:5,cbind(erro1,erro2,erro3,erro4,erro5), xlab= "Código do modelo",
     pch=16, axes=FALSE, ylim=c(0.02,0.06), ylab="", type="b", col="black", main="")
axis(2, ylim=c(0.02,0.06), col="black", las=1)  ## las=1 faz etiquetas horizontais
mtext("MAE", side=2, line=2.5)
box()
par(new=TRUE)
## Plotando o segundo gráfico e colocando a escala do eixo à direita
plot(1:5,cbind(erro6,erro7,erro8,erro9,erro0), xlab= "Código do modelo",
     pch=15, ylab="", ylim=c(0.02,0.06), axes=FALSE, type="b", col="red")
axis(2, ylim=c(0.02,0.06), lwd=2, line=3.5, col="red", col.axis = "red")
axis(1, at = 1:5, labels = c("Pregibon","Gosset","Aranda-Ordaz","Guerrero-Johnson","Stukel"), tick = FALSE)
mtext("RMSE", side=2, line=5.3, col="red")
grid()

Concluímos então que, utilizando ambas as medida de precisão, o modelo no qual minimizamos o erro de estimação é àquele no qual utilizamos a função de ligação Stukel, ou seja, o modelo no objeto R ajuste5. Observemos que neste caso é onde foi necessário utilizar mais parámetros.

Decidimos utilizar somente algumas das funções de ligação disponíveis no pacote de funções glmx.