Modelos de regressão para dados categóricos multinomiais

Prévia do material em texto

Aula 7: Modelos de regressão para dados categóricos multinomiais
	Introdução
	A distribuição multinomial
	O modelo de regressão logística multinomial para dados categóricos não ordenados
	Ilustração
	Modelos de regressão para dados categóricos ordinais
	Ilustração
Introdução
Em alguns estudos médicos e/ou epidemiológicos, especialmente aqueles baseados em questionários (inquéritos epidemiológicos, por exemplo), poderá haver interesse na avaliação do papel de determinantes de desfechos categóricos, e.g. auto-percepção de saúde, uso de contraceptivo, indicadores de qualidade de vida, etc.
O modelo probabilístico apropriado para descrever a variabilidade deste tipo de desfecho foi introduzido anteriormente, a saber a distribuição multinomial. Esta distribuição é uma extensão, ou generalização, da distribuição binomial, quando a variável aleatória de interesse tomar mais do que duas categorias.
Existem diferentes tipos de modelos de regressão para desfechos multicategóricos, três dos quais serão discutidos neste texto: (i) o modelo multinomial para desfechos categóricos sem ordenação, também conhecido como modelo “mlogito”; (ii) o modelo multinomial para desfechos categóricos com ordenação, também conhecido como “ologito”; e (iii) um modelo intermediário entre o “mlogito” e o “ologito”, conhecido como o modelo estereótipo.
 
Algumas ilustrações serão apresentadas ao longo do texto e os resultados dos ajustes das regressões comentados.
O modelo de regressão logística multinomial para dados categóricos não ordenados
Na introdução deste modelo de regressão, considera-se primeiro o modelo mais simples, aquele com apenas uma variável explicativa. Daqui pra frente, a seguinte notação será usada: k denota o número de categorias da variável desfecho e j será usado para indexar as categorias. A letra i como sempre, serve como índice para as observações da amostra, em geral indivíduos.
Supõe-se naturalmente, como em outros modelos introduzidos até aqui, que as observações da amostra não são correlacionadas. Este pressuposto pode ser relaxado, é claro, porém dependendo do tipo de estrutura de dependência imposto aos dados, diferentes classes de modelos vão ser mais ou menos apropriados. O modelo de regressão logística multinomial simples é descrito da forma abaixo.
Na equação deste modelo, assume-se que a k-ésima categoria é a referência para todas as outras. A escolha da referência é arbitrária, mas em geral esta escolha deve ser feita de forma a facilitar a interpretação dos resultados do ajuste do modelo. Por exemplo, em estudos epidemiológicos esta escolha geralmente incide no grupo de indivíduos para o qual o desfecho de saúde está ausente, ou o grupo com o menor grau de severidade do desfecho de saúde. Em cada situação haverá uma forma mais apropriada desta escolha, porém é importante notar que do ponto de vista estatístico modelos ajustados com diferentes grupos de referência são totalmente equivalentes. O que difere é a interpretação dos resultados.
Note que na equação do modelo de regressão logística multinomial para cada categoria do desfecho, menos a referência, dois parâmetros são necessários para descrever a associação da função dos valores esperados do desfecho com a covariável. Como em toda a classe de modelos lineares generalizados, estes parâmetros são estimados pelo critério da máxima verossimilhança e são interpretados de forma análoga ao modelo de regressão logística, i.e. o valor da função exponencial aplicada aos coeficientes são as razões de chances (odds ratios) relativas ao acréscimo de uma unidade nos níveis da covariável .
Considere agora o mesmo modelo com um número maior de covariáveis, i.e. o modelo de regressão logística multinomial multivariado1 No jargão estatístico, o termo multivariado se refere a um desfecho múltiplo, ou seja ao invés de apenas uma variável resposta sendo avaliada na amostra de indivíduos, várias respostas são avaliadas e modeladas simultaneamente. Por outro lado, quando o número de covariáveis no modelo de regressão é maior do que um usa-se o termo múltiplo. Finalmente o termo multinomial é usado para variáveis categóricas com mais de duas categorias. . Na versão mais geral, este modelo será descrito da forma abaixo.
Ou seja, na versão mais completa deste modelo, haveriam parâmetros para cada comparação (contraste) realizada com a categoria de referência, portanto haveriam parâmetros no modelo completo.
Por outro lado, as equações dos diferentes contrastes não necessariamente devem ser idênticas em sua composição em termos de covariáveis. Na prática, cada equação pode ter uma configuração específica, de acordo com o modelo teórico utilizado e/ou segundo a significância estatística dos efeitos das covariáveis, permitindo assim que qualquer covariável seja usada no j-ésimo contraste sem que seja incluída no preditor linear de outro contraste.
Alguns aplicativos de estatística permitem que a opção de seleção específica de covariáveis seja utilizada, mas normalmente a opção padrão considera a mesma especificação para todos os preditores lineares. Há vantagens e desvantagens em manter a mesma configuração do preditor linear em todos os contrastes.
De qualquer forma, com um pouco de álgebra é fácil derivar, a partir da equação do modelo de regressão logística multinomial, a expressão para a probabilidade condicional do i-ésimo indivíduo estar na j-ésima categoria.
Para completar o método de cálculo das probabilidades (prevalências) de categorias do desfecho a probabilidade de ocorrer a k-ésima categoria, para um dado perfil de covariável, é obtida pela propriedade de probabilidades complementares, i.e.
Note que omitir alguma covariável na composição de um preditor linear específico significa tornar o seu efeito nulo na equação da previsão das probabilidades.
Ilustração
O ajuste do modelo apresentado acima será ilustrado com um estudo sobre os determinantes da escolha do tipo de contraceptivo em uma região de Bangladesh, em 1988.
A variável desfecho, uso de contraceptivo (uso4), foi categorizada da seguinte maneira:
	esterilização (homem ou mulher);
	método reversível moderno (pílula anticoncepcional ou DIU);
	métodos tradicionais / ineficientes;
	não utilizava anticoncepcional.
As variáveis do banco de dados deste estudo armazenadas no arquivo “bangladesh.txt”, bem como suas descrições, são listadas abaixo.
	mulher
	ID para a mulher
	distrito
	ID para o distrito de residência
	uso
	Status de uso de contraceptivo: 1 = uso de algum método contraceptivo; 0 = não usa contraceptivo
	uso4
	Status de uso de contraceptivo, como descrito acima
	nfv
	Número de filhos vivos
	idade
	Idade da mulher, centrada na média amostral de 30 anos
	urbano
	Tipo de região da residência: 1 = urbano, 0 = rural
	educ
	Escolaridade: 0 = nenhuma; 1 = baixa; 2 = media; 3 = alta
	relig
	Religião: 1 = hindu; 0 = muçulmana
	prop_alfab
	Proporção de mulheres alfabetizadas no distrito de residência
	prop_rezar
	Proporção de mulheres no distrito de residência que rezam diariamente (medida de religiosidade)
As duas últimas variáveis são indicadores do contexto no qual os indivíduos estão inseridos e são usadas em modelos multiníveis.
Em alguns aplicativos existe um comando único para ajustar o modelo de regressão logístico multinomial para desfechos categóricos não ordenados. Contudo, na prática o ajuste deste modelo equivale a realizar regressões logísticas, sempre considerando a mesma categoria de referência. Assim na plataforma R, bem como em qualquer outro aplicativo de estatística, esta é uma forma alternativa para ajustar este modelo. 
Portanto, para analisar o perfil de uso de contraceptivo em Bangladesh, em 1988, os seguintes comandos podem ser usados para ajustar, por exemplo um modelo logístico multinomial para o desfecho tipo de anticontraceptivo, considerando apenas as covariáveis “idade” e “número de filhos vivos”, por enquanto.
bang <- read.table("bangladesh.txt")
names(bang) <- c("mul","dist","uso","uso4","nfv","idade","urb","edu","rel")names(bang)
attach(bang)
ccep1 <- ifelse(uso4==1,1,0)
ccep2 <- ifelse(uso4==2,1,0)
ccep3 <- ifelse(uso4==3,1,0)
bang <- cbind(bang,ccep1,ccep2,ccep3)
table(ccep1,uso4)
table(ccep2,uso4)
table(ccep3,uso4)
table(ccep1[uso4==4|uso4==1],nfv[uso4==4|uso4==1])
table(ccep2[uso4==4|uso4==2],nfv[uso4==4|uso4==2])
table(ccep3[uso4==4|uso4==3],nfv[uso4==4|uso4==3])
m1 <- glm(ccep1 ~ nfv + idade, subset(bang,uso4==1 | uso4==4), family=binomial(link="logit"))
m2 <- glm(ccep2 ~ nfv + idade, subset(bang,uso4==2 | uso4==4), family=binomial(link="logit"))
m3 <- glm(ccep3 ~ nfv + idade, subset(bang,uso4==3 | uso4==4), family=binomial(link="logit"))
Os resultados dos três ajustes estão relacionados a seguir.
 Call: glm(formula = ccep1 ~ nfv + idade, family = binomial(link = "logit"), data = subset(bang, uso4 == 1 | uso4 == 4)) Deviance Residuals: Min 1Q Median 3Q Max -0.8086 -0.6994 -0.4296 -0.3318 2.4332 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.678737 0.162894 -16.445 < 2e-16 *** nfv 0.465056 0.075729 6.141 8.2e-10 *** idade 0.017553 0.009088 1.931 0.0534 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1707.5 on 2029 degrees of freedom Residual deviance: 1597.9 on 2027 degrees of freedom AIC: 1603.9 Number of Fisher Scoring iterations: 5 
 Call: glm(formula = ccep2 ~ nfv + idade, family = binomial(link = "logit"), data = subset(bang, uso4 == 2 | uso4 == 4)) Deviance Residuals: Min 1Q Median 3Q Max -1.0906 -0.7959 -0.6958 -0.3956 2.3430 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.772399 0.115757 -15.311 < 2e-16 *** nfv 0.340367 0.057170 5.954 2.62e-09 *** idade -0.060409 0.008452 -7.147 8.86e-13 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 2532.4 on 2282 degrees of freedom Residual deviance: 2474.9 on 2280 degrees of freedom AIC: 2480.9 Number of Fisher Scoring iterations: 4 
 Call: glm(formula = ccep3 ~ nfv + idade, family = binomial(link = "logit"), data = subset(bang, uso4 == 3 | uso4 == 4)) Deviance Residuals: Min 1Q Median 3Q Max -0.6681 -0.6556 -0.4887 -0.4172 2.2442 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.413376 0.149916 -16.098 < 2e-16 *** nfv 0.337244 0.073471 4.590 4.43e-06 *** idade -0.001732 0.009335 -0.186 0.853 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1630.1 on 2009 degrees of freedom Residual deviance: 1590.6 on 2007 degrees of freedom AIC: 1596.6 Number of Fisher Scoring iterations: 5 
O conjunto de linhas de comandos abaixo calcula as razões de chances para as covariáveis do modelo ajustado. Por não haver interesse direto, em geral, na interpretação da razão de chances associado ao intercepto, estas serão omitidas dos resultados. No script abaixo, o objeto denominado “m”, disposto nas duas primeiras linhas de comandos, representam o modelo logístico ajustado.
 beta <- coef(m) se.beta <- sqrt(diag(summary(m)$cov.unscaled)) raz_z <- beta / se.beta pval <- pnorm(abs(raz_z),lower.tail=FALSE,log.p=FALSE) or <- exp(beta) liminf <- exp(beta-1.96*se.beta) limsup <- exp(beta+1.96*se.beta) data.frame("Odds_ratio"=or, "razão-z"=raz_z, "p-valor"=pval,"LI_IC95"=liminf,"LS_IC95"=limsup) 
As razões de chances para as covariáveis número de filhos vivos e idade estão listados abaixo, para os três contrastes disponíveis nesta análise.
 Odds_ratio razão.z p.valor LI_IC95 LS_IC95 nfv 1.5921031 6.141046 4.098992e-10 1.37249032 1.84685632 idade 1.0177076 1.931406 2.671645e-02 0.99974017 1.03599788 
 Odds_ratio razão.z p.valor LI_IC95 LS_IC95 nfv 1.4054629 5.953638 1.311229e-09 1.2564803 1.5721104 idade 0.9413799 -7.147215 4.427809e-13 0.9259135 0.9571047 
 Odds_ratio razão.z p.valor LI_IC95 LS_IC95 nfv 1.4010812 4.5901482 2.214657e-06 1.21317445 1.6180925 idade 0.9982693 -0.1855543 4.263971e-01 0.98017043 1.0167025 
Para o número de filhos vivos, a razão de chances relacionada ao desfecho homem / mulher esterilizada (1) vs. não usa contraceptivo (0), é igual a 1.56 para cada filho vivo a mais. Considerando os dois próximos contrastes, as razões de chances correspondentes são iguais a 1.41. Para a variável idade, a razão de chances do desfecho método reversível moderno (1) vs. nao usa contraceptivo (0) é igual a 0.94 e é estatisticamente significativa, ou seja para cada ano a mais na idade da mulher há uma diminuição na chance dela pertencer ao grupo que usa métodos reversíveis modernos com respeito ao grupo que não usa contraceptivos.
O pacote nnet possui uma função que ajusta o modelo multinomial para categorias não ordinais de forma direta. Para os dados deste exemplo, as seguintes linhas de comando devem ser executadas.
 library(nnet) us <- ifelse(uso4==4,0,uso4) m <- multinom(us ~ nfv + idade) summary(m) 
A segunda linha de comando serve apenas para modificar a codificação da variável desfecho de forma que a categoria “näo usa contraceptivo” seja colocada como referência, uma vez que o comando multinom considera automaticamente como referência o valor numérico mais baixo da variável desfecho. Os resultados do ajuste estão abaixo. Note que os valores numéricos são bastante próximos aos obtidos nas regressões logísticas separadas.
 Call: multinom(formula = us ~ nfv + idade) Coefficients: (Intercept) nfv idade 1 -2.694051 0.4709955 0.0197571341 2 -1.789346 0.3480600 -0.0631948465 3 -2.421547 0.3410286 0.0007136124 Std. Errors: (Intercept) nfv idade 1 0.1644997 0.07666769 0.009300547 2 0.1162456 0.05732006 0.008625155 3 0.1512711 0.07433410 0.009734456 Residual Deviance: 6028.138 AIC: 6046.138 
O modelo de regressão multinomial ordinal
Quando o desfecho categórico de interesse principal possui uma ordenação natural o modelo de regressão pode explorar tal ordenação na análise de dados, proporcionando um tipo diferente de interpretação de resultados. Isto pode ser feito de diferentes formas. A primeira versão, disponível na maior parte dos aplicativos estatísticos, é o chamado modelo de chances proporcionais. Neste modelo, a função de ligação logito também é utilizada, entretanto a forma de criar contrastes é modificada em relação ao modelo anterior. O modelo é descrito, de forma geral, como abaixo.
Seja a categoria correspondente ao i-ésimo indivíduo da amostra e o índice que representa as categorias da variável desfecho. Denota-se por o logaritmo da razão entre as seguintes probabilidades complementares2 Alguns autores parametrizam este modelo de outra forma, considerando no numerador a probabilidade do desfecho ser menor ou igual à j-ésima categoria e no denominador o oposto. Os dois modelos são equivalentes, é claro, mas a interpretação dos parâmetros deverá ser adaptada ao tipo de parametrização. .
O preditor linear deste modelo tem a seguinte configuração:
Note que o intercepto deste modelo varia de acordo com o contraste específico, mas o restante dos parâmetros do preditor linear é mantido fixo ao longo dos contrastes entre as probabilidades cumulativas. Este pressuposto é conhecido por chances proporcionais ou regressão paralela.
Alguns autores consideram este modelo mais apropriado para desfechos ordinais oriundos da categorização de uma variável contínua, e.g. a variável categórica status de índice de massa corporal descrita, por exemplo, por três categorias: eutrofia, sobrepeso e obesidade.
Ilustração
Os dados de um estudo amostral do tipo “survey”, realizado nos EUA e intitulado “Inventory of Positive Psychological Attitudes” estão no arquivo “survey_ippa.dta” O questionário consiste de 30 perguntas sobre as atitudesdas pessoas em situações corriqueiras, completadas por dados demográficos e perguntas adicionais sobre percepção de saúde e da forma física. As variáveis têm nomes óbvios ou estão categorizadas na própria base de dados com nomes bastante claros para as categorias. A tabela abaixo descreve o significado das variáveis.
	Variável
	Significado
	Gender
	Gênero
	Age
	Idade
	Marital
	Estado civil
	Educ
	Número de anos de educação
	Smoke
	Status de fumante
	Work
	Status de emprego
	Polaff
	Afiliação política
	Depress
	Freqüência do tempo em estado depressivo
	Satcurwt
	Satisfação com o peso atual
	Satwt18
	Satisfação com o peso aos 18 anos
	Health
	Estado geral de saúde
	Qolcur
	Qualidade de vida atual
	Qol18
	Qualidade de vida aos 18 anos
	Winter
	Tipo de férias de inverno que gostaria de realizar
	Ipa1
	Nível de energia individual corriqueiro
	Ipa2
	Forma de reação quando está sob pressão
	Ipa3
	Como você vê a vida como um todo?
	Ipa4
	Como você percebe o valor das atividades diárias?
	Confid
	Autoconfiança em situações de estresse
	Life
	Satisfação e objetivos de vida
	Ippatot
	Total do inventário de atitudes psicológicas positivas
 
As variáveis “Satcurwt”, “Satwt18”, “Health”, “Qolcur”, “Qol18”, “Ipa1”, “Ipa2”, “Ipa3” e “Ipa4” são escores que variam de 1 a 7 ou de 1 a 10. As escalas foram construídas de forma que respostas positivas aos itens correspondem a valores grandes da escala, respostas negativas a valores pequenos e respostas intermediárias a valores medianos da escala.
As variáveis “Ipa1” até “Ipa4” são relativas às quatro primeiras perguntas do inventário de 30 perguntas, enquanto as variáveis “confid”, “life” e “Ippatot” são somas dos escores das 30 perguntas. “Ippatot” é a soma de todos os escores, e “confid” e “life” são somas parciais dos escores de forma que os valores de “Ippatot” são iguais aos de “confid” mais os de “life”. As outras variáveis da base de dados são categóricas e suas categorias estão dispostas na base de dados, ou são contínuas (“age” e “educ”). 
Na ilustração sobre o modelo de regressão ordinal a variável desfecho considerada é a qualidade de vida atual, “Qolcur”, contudo as duas primeiras categorias foram agregadas. As covariáveis investigadas foram gênero, idade, estado civil (dicotomizado como morando junto ou morando sozinho) e anos de escolaridade.
Inicialmente, uma análise univariada será realizada. A função prop.trend.test implementa o teste qui-quadrado para tendência de proporções, ferramenta essencial para avaliar se uma covariável categórica ou discreta está associada com o desfecho categórico.