Regressão linear em Python para epidemiologistas em 6 passos

Neste tutorial, abordaremos as seguintes etapas:

1. Abre o conjunto de dados

2. Explorar os dados

3. Faça uma pergunta de investigação (que pode ser respondida usando um modelo de regressão linear)

4. Executa e interpreta uma regressão linear simples

5. Executa e interpreta uma regressão linear múltipla

6. Responde à pergunta de investigação

Vamos começar!

1. Abre o conjunto de dados

Vamos extrair uma versão disponível publicamente do estudo Framingham heart study https://www.kaggle.com/amanajmera1/framingham-heart-study-dataset/version/1

Neste site, você também podes consultar a codificação das diferentes variáveis. Por exemplo, em relação ao sexo (variável masculino; masculino=1 significa masculino e masculino=2 significa feminino)

Primeiro baixa o ficheiro para o teu computador.

import numpy as np 
import pandas as pd
import statsmodels.formula.api as smf

Se for tua primeira vez com o StatsModels, talvez seja necessário usar isso antes (instalar do Anaconda, distribuição Python).

conda install -c conda-forge statsmodels

Then open the Framingham database (or data frame). We will name it df to make things easier.

df = pd.read_csv('framingham.csv')

From now every time we want to “call” our data frame, we will use df in our code.

Vamos olhar para as primeiras 5 linhas e todas as nossas colunas.

df.head()

Vamos ver quantas colunas e linhas temos…

df.shape() 

(4240, 16)

O resultado significa que o quadro de dados tem 4.240 observações (linhas) e 16 variáveis (colunas). Lembra-te que .head() nos deu apenas um vislumbre das 4.240 observações/linhas.

Agora vamos descrever nossas variáveis

df.describe()

O parêntese vazio em df. describe() instrui o Python a mostrar todas as variáveis ​​disponíveis.

Aqui podemos ver o número de observações por coluna (contagem), seguido pela média e desvio padrão (dp). Observa que para variáveis ​​discretas binárias a interpretação é um pouco diferente e a média será igual à proporção relativa (0 a 1, em vez de 0–100%).

3. Faça uma pergunta de investigação

Existe associação entre pressão arterial sistólica e idade, após ajuste para fatores de confusão relevantes?

Como podemos responder a isso?

Aqui temos um resultado contínuo (pressão arterial sistólica em mmHg) e uma covariável contínua independente (idade em anos).

Como a pressão arterial sistólica (sBP) (variável chamada sysbp) tem uma distribuição normal e nós o que testar uma associação, devemos ir para regressão linear .

Consideraremos as seguintes variáveis ​​disponíveis como possíveis fatores de confusão relevantes para a associação entre sBP e idade; sexo (variável binária: masculino), escolaridade (variável categórica: educação), tabagismo (variável contínua: cigsPerDay para cigarros por dia), toma de medicamentos para hipertensão (variável binária: BPMeds) e colesterol total ( variável contínua: totChol).

Antes de correr os números, temos que definir uma hipótese nula e uma hipótese alternativa:

• Hipótese nula (H0): NÃO há associação entre sBP e idade (considerando os fatores de confusão escolhidos);
• Alternative hypothesis (H1): There is an association between sBP and age (considering the chosen confounders).

4. Executa e interpreta uma regressão linear simples

Vamos para uma regressão linear simples.

Descrevemos o modelo. Uma alternativa clássica em ciência de dados seria usar o pacote Scikilt-learn. Em alternativa, usaremos o pacote strong>StatsModels mais intuitivo para quem tem experiência em epidemiologia. Tem semelhanças com outros programas, como Stata e R /RStudio.

model1 = smf.ols(formula='sysBP ~ age', data=df).fit()

• smf chama o pacote StatsModels
• olhar para o Python que estamos com uma regressão mínimos quadrados ordinários (OLS) (um tipo de regressão linear)
• formula= usado para escrever as variáveis ​​dependentes e independentes
• primeira variável entre parênteses antes de “~”/variável dependente/resultado: A primeira variável é nossa única variável dependente. Este é o nosso resultado, a variável que determina qual tipo de regressão utilizar e aquela a ser associada a todas as outras covariáveis;
• ~ dentro de parênteses: Marca a fronteira entre o resultado (variável dependente) à esquerda e as covariáveis ​​(variáveis ​​independentes) à direita;
• covariáveis ​​independentes/variáveis ​​independentes: todas as outras variáveis ​​após o "~", entre parênteses;
• sinal + entre parênteses: o sinal + é usado para separar diferentes variáveis ​​independentes dentro do mesmo modelo (útil para modelos multivariáveis, também conhecidos como muitas variáveis ​​independentes)
• ,data=: This marks the name of the data frame
• .fit() diz ao Python que queremos preencher nossa função (“executar a função”)

model1 = smf.ols(formula='sysBP ~ age', data=df).fit() 
print(model1.summary())

Na metade superior da saída (dentro do quadrado vermelho) podemos ver alguns dados gerais importantes sobre o modelo.

Na coluna da esquerda, vemos

• A variável dependente é sysBP
• Estamos executando um mínimo quadrado ordinário (OLS) modelo de regressão (uma forma de regressão linear)
• Incluímos 4.240 observações (pacientes)

Na coluna da direita, podemos ver algumas informações importantes sobre o diagnóstico do modelo:

• R-quadrado = 0,155, o que significa 16% da variabilidade do resultado é explicada pelo modelo. O R-quadrado varia de 0 a 1 e mede quão próximos os dados estão da linha de regressão ajustada. Em outras palavras, indica com % que a variabilidade do resultado (dados de resposta) é explicada pelo modelo;
• O valor da estatística F é usado para o cálculo do valor p do modelo, Prob (estatística F), que aqui é <0,05. Isso também nos diz que o Python está usando um teste ANOVA, o que implica uma distribuição F.

Como estamos realizando um estudo epidemiológico para observar/testar uma associação, em vez de construir um modelo preditivo, não dividiremos nossos dados em treino e teste forte>. Portanto, não te surpreendas se essa abordagem parecer diferente do que encontrarás na maioria dos tutoriais de ciência de dados.

Em relação ao output na metade inferior

Vamos ler o output, coluna por coluna. Aqui estamos a avaliar a associação entre o nosso resultado (sysbp) e a nossa variável independente, idade (primeira coluna), que tem uma associação positiva de 1,0128 (ou coeficiente, ou segunda coluna), com um erro padrão de 0,036 (terceira coluna), com um t-ratio (para um t-distribution) de 27,911 (quarta coluna) que fornece um valor p <0,001 (quinta coluna) , correspondendo a um intervalo de confiança de 95% para o coeficiente beta entre 0,94 2 (limite inferior, a sexta coluna) e 1,084 (limite superior, sétima coluna).

Para tornar as coisas mais simples…

Podemos rejeitar a hipótese nula de não associação entre idade e PAs. Existe uma associação positiva (coeficiente beta>1,0) entre idade e pressão arterial sistólica, que é estatisticamente significativa (p<0,01 e IC 95% não inclui coeficiente=0,0).

O coeficiente beta pode assumir um valor entre -∞ e +∞.

• Um coeficiente positivo representa uma associação linear positiva;
• Um coeficiente negativo representa uma associação linear negativa.

Em relação ao Intervalo de confiança de 95% (IC de 95%) para o coeficiente beta

• Se o IC 95% incluir 0, a associação não é estatisticamente significativa (p≥0,05), a hipótese nula de nenhuma associação não pode ser rejeitada;
• Se o IC 95% não incluir 0, a associação é estatisticamente significativa (p<0,05), a hipótese nula de associação pode ser rejeitada.

E a “(Interceção)”?

Neste caso (Intercept)= 82.14 representa a interceptação, ou o β0. Isso significa que, quando idade=0, sBP=82mmHg. Nesse contexto, o β0 é irrelevante para nossa questão de pesquisa, porque NÃO estamos a construit um modelo preditivo. Estamos a avaliar uma associação.

5. Executa e interpreta uma regressão linear múltipla

O processo aqui é bem simples, basta adicionar todas as covariáveis ​​(variáveis ​​independentes) após sua variável independente (neste caso, após sysbp).

model2 = smf.ols(formula='sysBP ~ age + male + education + cigsPerDay + BPMeds + totChol', data=df).fit()print(model2.summary())

Antes de ler o resultado, vemos que este é um modelo válido, mas… não é bem o que queríamos… Variáveis ​​discretas com múltiplas categorias (por exemplo, educação) não são estratificadas por grupos. Podemos especificar o tipo de covariáveis.

model3 = smf.ols(formula='sysBP ~ age + male + C(education) + cigsPerDay + BPMeds + totChol', data=df).fit()print(model3.summary())

Ótimo, agora podemos passar para a interpretação do nosso output.

Vamos dar uma olhada nas covariáveis ​​contínuas

Em primeiro lugar, observa que há um valor ligeiramente diferente para coeficiente e IC 95% para idade, quando comparado à regressão linear simples anterior. Este é um fenômeno normal e esperado. Conforme observado anteriormente, há uma regressão linear positiva estatisticamente significativa.

Repara que temos apenas 4007 observações. Isso ocorre porque 233 pacientes (4240–4007) tinham pelo menos uma observação ausente para qualquer uma das variáveis ​​do modelo.

Lembra-te que o OLS clássico só pode ser executado com dados completos, o que significa que sempre que uma observação (linha) tiver dados ausentes para pelo menos uma covariável, essa linha não entrará no modelo.

How do we read the output; taking as example age:

• The null hypothesis (H0) is: In the population, there is no linear association between age and sBP, when controlling for gender, education, smoking, anti-hypertensors intake, and total cholesterol (in other words, when controlling for all other covariates present in the model);
• Since p-value<0.05, we reject H0, meaning there is a linear association between age and sBP, when controlling for gender, education, smoking, anti-hypertensors intake, and total cholesterol;
• Interpretation of the slope: For each 1 unit increase in age, mean sBP increases +0.84mmHg.

Notice an important pattern regarding p-values and 95%CI

• When p<0.05, the confidence interval does NOT include 0.0 (eg: age, beta coefficient (95%CI): 0.84 [0.76;0.91]);
• When p≥0.05 the confidence interval includes 0.0 (eg; male, beta coefficient (95% CI): -0.087 [-1.40;+1.22]).

Vamos avaliar as covariáveis ​​binárias discretas

Vejamos o gênero (male = 0 significa paciente do sexo feminino; male= 1 significa paciente do sexo masculino):

• A hipótese nula (H0) é: Na população, não há diferença entre homens e mulheres no valor médio da sBP, quando controlados por idade, escolaridade, tabagismo, ingestão de anti-hipertensores , e colesterol total;
• Aqui coeficiente beta (IC 95%) = -0,09 (-1,4;+1,23). Como o IC de 95% inclui 0,0 (e p≥0,05), não rejeitamos H0.

Vamos avaliar as covariáveis ​​não binárias discretas

Agora chegamos à parte mais intrigante do modelo... A covariável não-binária discreta; Education (Educação).

Mas, o que está em falta? Category one (education=1)

Por quê? Como todas as outras categorias são comparadas em relação à categoria de referência (educação=1)

df.education.value_counts()

Como vemos, existem 4 categorias (1, 2, 3, 4). Na base de dados original

• 1 = Algum ensino secundário
• 2= ​​ensino secundário ou equivalente
• 3= Algum tempo de faculdade ou escola profissional
• 4= Ensino Superior

Por convenção, o modelo usa a categoria com o menor índice (neste caso, 1, algum ensino médio) como referência. Todas as categorias restantes no modelo são comparadas com a categoria de referência.

Mas, como faço para ler os coeficientes para a educação?

No caso de educação, olhando primeiro para a primeira linha, quando educação=2

• A hipótese nula (H0) é: Na população, não há diferença entre aqueles com “algum ensino médio” (educação=1) e aqueles com “ensino médio” (educação=2 ) em seu valor médio para sBP, quando controlados por idade, sexo, tabagismo, ingestão de anti-hipertensores e colesterol total;
• Como valor-p≥0,05, rejeitamos H0, não rejeitamos a hipótese nula para essa categoria.

E as diferentes categorias de educação?

Vemos que para a 3ª categoria para Educação (“Algum tempo de ensino superior”), o coeficiente beta (IC 95%): -3,57 (-5,35;-1,79), significando que; "Em comparação com pacientes com algum ensino secundário (categoria de educação =1), os pacientes com algum ensino superior (educação =3) têm uma sBP média mais baixa, e essa diferença é estatisticamente significativa".

Consideraque neste modelo não podemos comparar a associação entre pacientes na categoria Education=2 e pacientes na categoria Education=3. Para isso, teríamos que mudar a referência.

6. Responde à pergunta de investigação

Existe associação entre pressão arterial sistólica e idade?

Sim, existe uma associação positiva entre idade e pressão arterial sistólica (coeficiente Beta [95%IC]: 0,84[0,76;0,91]), o que é estatisticamente significativo (p<0,05 e intervalo de confiança não include 0), mesmo considerando o efeito de possíveis confundidores relevantes (covariáveis ​​no modelo).