1

Olá!

Eu tenho um conjunto de dados, com o qual pretendo realizar um teste para analisar a influência de algumas variáveis preditoras sobre uma variável resposta; como existem muitos zeros na minha variável resposta (são 766 zeros de 2830 unidades amostrais), decidi utilizar a abordagem Hurdle Model. Lá no R, eu escrevi esses comandos:

fórmula <- dados$BC ~ dados$z_primeiro_artigo +
 dados$z_capacidade_científica + dados$z_tamanho_corporal +
 z_reproduções_por_ano + dados$Red_List_Status +
 dados$Tipo_de_desenvolvimento | dados$z_capacidade_científica +
 dados$z_tamanho_corporal + z_reproduções_por_ano +
 dados$Red_List_Status + dados$Tipo_de_desenvolvimento

resultado <- hurdle(formula = fórmula, dist = "negbin", data = dados, na.action = "na.fail")
summary(resultado)

Call:
hurdle(formula = fórmula, data = dados, na.action = "na.fail", dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.1840 -0.6896 -0.2369  0.1864 16.3096 

Count model coefficients (truncated negbin with log link):
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                           89.2998674  0.1065855 837.824  < 2e-16 ***
dados$z_primeiro_artigo               -0.0475314         NA      NA       NA    
dados$z_capacidade_científica          0.0751863  0.0048415  15.530  < 2e-16 ***
dados$z_tamanho_corporal               0.0020403  0.0006407   3.185  0.00145 ** 
z_reproduções_por_ano                  0.1797664  0.0761702   2.360  0.01827 *  
dados$Red_List_StatusEN               -0.4140505  0.1725280  -2.400  0.01640 *  
dados$Red_List_StatusLC                0.2434877  0.1372437   1.774  0.07604 .  
dados$Red_List_StatusNT               -0.2326801  0.1856711  -1.253  0.21014    
dados$Red_List_StatusVU                0.0002679  0.1702307   0.002  0.99874    
dados$Tipo_de_desenvolvimentoLarval    0.4254052  0.0928358   4.582  4.6e-06 ***
dados$Tipo_de_desenvolvimentoVivípara  0.0109588  0.3846127   0.028  0.97727    
Log(theta)                            -1.1538934  0.1093832 -10.549  < 2e-16 ***
Zero hurdle model coefficients (binomial with logit link):
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            1.3147054  0.2712539   4.847 1.25e-06 ***
dados$z_capacidade_científica          0.0682073  0.0100039   6.818 9.23e-12 ***
dados$z_tamanho_corporal               0.0015036  0.0008404   1.789   0.0736 .  
z_reproduções_por_ano                  0.3522174  0.2009335   1.753   0.0796 .  
dados$Red_List_StatusEN               -0.4264203  0.1776977  -2.400   0.0164 *  
dados$Red_List_StatusLC               -0.1618832  0.1555683  -1.041   0.2981    
dados$Red_List_StatusNT               -0.2458956  0.2064901  -1.191   0.2337    
dados$Red_List_StatusVU               -0.2674147  0.1880392  -1.422   0.1550    
dados$Tipo_de_desenvolvimentoLarval    0.0385487  0.0989498   0.390   0.6968    
dados$Tipo_de_desenvolvimentoVivípara  0.1392403  0.4588545   0.303   0.7615    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 0.3154
Number of iterations in BFGS optimization: 27 
Log-likelihood: -5853 on 22 Df
Warning message:
In sqrt(diag(object$vcov)): NaNs produzidos

Notem que os valores da variável "z_primeiro_artigo" aparecem como "NA", e eu não compreendi essa mensagem de aviso no final: "In sqrt(diag(object$vcov)): NaNs produzidos". Alguém saberia como me dar um help?

1 Resposta 1

2

Modelos lineares generalizados não fazem mágica. Não adianta ter dados, tentar ajustar um modelo a eles e acreditar que tudo vai dar certo. Além disso, é muito difícil (talvez impossível) te dar uma resposta definitiva sem trabalhar com os mesmos dados que tu está utilizando. Entretanto, é possível levantar algumas hipóteses sobre o que talvez esteja acontecendo.

0) Antes de jogar os dados num modelo, faça uma análise exploratória. Plote eles. Faça estatísticas simples, como média e desvios padrão para variáveis quantitativas e tabelas de frequências para variáveis categóricas. Isto vai ajudar a definir maneiras melhores de resolver problemas que podem surgir nas tuas análises.

1) Eu contei 6 variáveis para a modelagem das contagens e 5 para o excesso de zeros. Está correto isso? Tem algum motivo para excluir z_primeiro_artigo da modelagem extra de zeros? A parte de modelagem do excesso de zeros precisa ser tão complexa assim? De todo modo, com 6 covariáveis, é possível que entre estas 6 variáveis preditoras, algum par delas tenha correlação alta. Isto cria um problema chamado multicolinearidade. Pesquise sobre isso e veja como ela pode afetar a tua regressão.

2) z_primeiro_artigo tem erro padrão igual a NA. Isto significa que não foi possível calcular a variabilidade do erro de estimação deste parâmetro. Verifique se z_primeiro_artigo é constante. O fato de não haver variação nesta covariável pode ser o motivo para isto estar ocorrendo.

3) In sqrt(diag(object$vcov)): NaNs produzidos significa que alguns elementos da diagonal da matriz object$vcov são negativos. Verifique se diag(resultado$vcov) possui números negativos. Se possuir, isto significa que a matriz Hessiana do modelo não é positiva definida. Uma maneira de resolver este problema é verificar se os teus dados estão na mesma escala. Por exemplo, alguma covariável pode estar na ordem das unidades e, outra, na ordem de centenas. Isso quase sempre dá problema na hora de ajustar modelos lineares. Veja como transformar dados usando scaling. Apenas tome cuidado que as inferências feitas a partir de dados transformados são diferentes das inferências feitas nos dados originais.

Como pode ver, nenhuma resposta minha é definitiva. Este problema não é simples e fica impossível dar um diagnóstico preciso sem ter acesso aos dados. Por fim, não acho que o tamanho amostral seja problema. 2830 é um tamanho bastante razoável para este caso, com este número de covariáveis.

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .