Etude de la nocivité du tabac chez les femmes

1.7.1. Distribution des âges des participantes

2.1.1. Calcul des paramètres

rounding = 1

smokers = data[(data['Smoker'] == 'Yes')]
nonsmokers = data[(data['Smoker'] == 'No')]

dead_smokers = smokers[(smokers['Status'] == 'Dead')]
dead_nonsmokers = nonsmokers[(nonsmokers['Status'] == 'Dead')]

smokers_death_rate = round(dead_smokers.shape[0] / smokers.shape[0] * 100, rounding)
nonsmokers_death_rate = round(dead_nonsmokers.shape[0] / nonsmokers.shape[0] * 100, rounding)

def death_rate_ci95(sample_size: int,
                    deads: int,
                    per_people: int = 100) -> tuple:
    """Compute the 95% confidence interval upper and lower limits of a mortality rate in a sample.

    :param sample_size: number of people in the sample.
    :param deads: number of deads people in the sample.
    :param per_people: standardisation number: deads for per_people people. 
    """
    from math import sqrt

    
    upper = per_people / sample_size * (deads + 1.96 * sqrt(deads))
    lower = per_people / sample_size * (deads - 1.96 * sqrt(deads))

    return lower, upper

smokers_dr_95ic = death_rate_ci95(smokers.shape[0], dead_smokers.shape[0])
nonsmokers_dr_95ic = death_rate_ci95(nonsmokers.shape[0], dead_nonsmokers.shape[0])

2.1.2. Résultats

2.2.1. Calcul des paramètres et résultats

categories = [(18, 34), (34, 54), (54, 64), (64, 123)]

death_rates = []  # plotting purpose
dr_uncertainties = []  # plotting purpose

print(f"|Classe d'âge|Tabagisme|Effectif total|Individus décédés|\
        Taux de mortalité (%)|IC95%|\n"
      f"|-+-+-+-+-+-|")

for category in categories:
    for smoker in data['Smoker'].unique():
        
        alive, dead = nan, nan
        
        if smoker == 'Yes':
            smoke = 'Fumeuses'
        else:
            smoke = 'Non-fumeuses'
        
        for status in data['Status'].unique():
            
            subset = data[(data['Age'] >= category[0])
                          & (data['Age'] < category[1])
                          & (data['Smoker'] == smoker)
                          & (data['Status'] == status)]
            
            if status == 'Alive':
                alive = subset.shape[0]
            else:
                dead = subset.shape[0]

        death_rate = dead / (alive + dead) * 100
        death_rates.append(death_rate)

        dr_ci95 = death_rate_ci95(alive + dead, dead)
        dr_uncertainties.append(dr_ci95[1] - death_rate)

        rounded_dr_ci95 = list(dr_ci95)
        rounded_dr_ci95[0] = round(rounded_dr_ci95[0], rounding)
        rounded_dr_ci95[1] = round(rounded_dr_ci95[1], rounding)
        r_dr_ci95 = tuple(rounded_dr_ci95)
        
        print(f"|[{category[0]};{category[1]}[|{smoke}|{alive + dead}|\
                {dead}|{round(death_rate, rounding)}|{r_dr_ci95}|")

fig, ax = plt.subplots()

x_pos = arange(len(categories)) 
width = 0.35

xlabels = list(map(str, categories))
for x_l, xlabel in enumerate(xlabels):
    head, _, tail = xlabel.partition(',')
    xlabels[x_l] = '[' + head[1:] +  ';' + tail[:-1] + '['

plt.bar(x=x_pos - width / 2,
        height=death_rates[0:len(death_rates):2],
        yerr=dr_uncertainties[0:len(death_rates):2],
        label='Smoker', width=width, edgecolor='k',
        error_kw=dict(lolims=True))

plt.bar(x=x_pos + width / 2,
        height=death_rates[1:len(death_rates):2],
        yerr=dr_uncertainties[1:len(death_rates):2],
        label='Nonsmoker', width=width, edgecolor='k',
        error_kw=dict(lolims=True))

# In order to modify the default lolims error bars shape
for ch in ax.get_children():
    if str(ch).startswith('Line2D'):
        ch.set_marker('_')
        ch.set_markersize(6)

plt.ylabel('Death rate (%)')
plt.xlabel('Age category')
plt.xticks(x_pos)
ax.set_xticklabels(xlabels)
plt.legend()
plt.tight_layout()
plt.savefig("age_categories_death_rate.png")

2.3.1. Manipulation des données

deaths = []

for status in data['Status']:
    if status == 'Alive':
        deaths.append(0)
    else:
        deaths.append(1)

data['Death'] =  deaths

smokers = data[(data['Smoker'] == 'Yes')]
nonsmokers = data[(data['Smoker'] == 'No')]
print(smokers.head(), '\n', '\n', nonsmokers.head())

  Smoker Status   Age  Death
0    Yes  Alive  21.0      0
1    Yes  Alive  19.3      0
4    Yes  Alive  81.4      0
7    Yes   Dead  57.5      1
8    Yes  Alive  24.8      0 
 
    Smoker Status   Age  Death
2      No   Dead  57.5      1
3      No  Alive  47.1      0
5      No  Alive  36.8      0
6      No  Alive  23.8      0
11     No   Dead  66.0      1

2.3.2. Régressions logistiques

1. Réglages et code de la modélisation
   

   Pour effectuer la régression logistique, nous utilisons la bibliothèque statsmodels. Nous laissons les paramètres par défaut : le modèle est le logit, la méthode est celle du maximum de vraisemblance (Maximum Likelyhood Estimation).

   def logistic_regression(independent_variable: list,
                           dependent_variable: list) -> object:
       """Logisitic regression with Logit model and MLE method.
   
       :param independent_variable: the independent variable
       :param dependent_variable: the dependent variable
       
       :return fitted_model: the logisitc regression fitted model.
       """
       import statsmodels.api as sm
   
       model = sm.Logit(dependent_variable, sm.add_constant(independent_variable))
       fitted_model = model.fit()
   
       print('\n', fitted_model.summary())
   
       return fitted_model
   

2. Cas des fumeuses
   

   proba_smokers, model_smokers = logistic_regression(smokers['Age'],
                                                      smokers['Death'])
   

   /usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
     return ptp(axis=axis, out=out, **kwargs)
   Optimization terminated successfully.
            Current function value: 0.412727
            Iterations 7
   
                               Logit Regression Results                           
   ==============================================================================
   Dep. Variable:                  Death   No. Observations:                  582
   Model:                          Logit   Df Residuals:                      580
   Method:                           MLE   Df Model:                            1
   Date:              jeu., 14 oct. 2021   Pseudo R-squ.:                  0.2492
   Time:                        10:05:54   Log-Likelihood:                -240.21
   converged:                       True   LL-Null:                       -319.94
   Covariance Type:            nonrobust   LLR p-value:                 1.477e-36
   ==============================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
   ------------------------------------------------------------------------------
   const         -5.5081      0.466    -11.814      0.000      -6.422      -4.594
   Age            0.0890      0.009     10.203      0.000       0.072       0.106
   ==============================================================================
   

3. Cas des non-fumeuses
   

   proba_nonsmokers, model_nonsmokers = logistic_regression(nonsmokers['Age'],
                                                            nonsmokers['Death'])
   

   /usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
     return ptp(axis=axis, out=out, **kwargs)
   Optimization terminated successfully.
            Current function value: 0.354560
            Iterations 7
   
                               Logit Regression Results                           
   ==============================================================================
   Dep. Variable:                  Death   No. Observations:                  732
   Model:                          Logit   Df Residuals:                      730
   Method:                           MLE   Df Model:                            1
   Date:              jeu., 14 oct. 2021   Pseudo R-squ.:                  0.4304
   Time:                        10:05:54   Log-Likelihood:                -259.54
   converged:                       True   LL-Null:                       -455.62
   Covariance Type:            nonrobust   LLR p-value:                 2.808e-87
   ==============================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
   ------------------------------------------------------------------------------
   const         -6.7955      0.479    -14.174      0.000      -7.735      -5.856
   Age            0.1073      0.008     13.742      0.000       0.092       0.123
   ==============================================================================
   

2.3.3. Visualisation des résultats

fig, ax = plt.subplots()
ax.set_xlim(15, 100)
ax.set_ylabel('Death probability')
ax.set_xlabel('Age')

lmplot(data=data, x='Age', y='Death', hue='Smoker',
       logistic=True, y_jitter=0.01, truncate=False,
       markers=['o', 'x'], scatter_kws=dict(alpha=0.3))

plt.tight_layout()

plt.savefig('death_probabilities.png')