fibro-height-weight-bmi-skinarea.ipynb

{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Download and Preprocess the Dataset",
   "id": "7466459b63682600"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:39:09.088834Z",
     "start_time": "2024-07-20T21:39:09.080608Z"
    }
   },
   "cell_type": "code",
   "source": [
    "### IMPORTS\n",
    "import pyreadstat\n",
    "import pandas as pd\n",
    "import pandas as pd\n",
    "from scipy.stats import pearsonr, spearmanr, pointbiserialr\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "from scipy import stats"
   ],
   "id": "a3f589a0a3dea36a",
   "outputs": [],
   "execution_count": 3
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Download the dataset",
   "id": "cc753af73f026bdb"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:39:28.194997Z",
     "start_time": "2024-07-20T21:39:09.090371Z"
    }
   },
   "cell_type": "code",
   "source": "original_df, meta = pyreadstat.read_sas7bdat(r'data_donnee\\hs.sas7bdat')",
   "id": "9e053b73b285fffa",
   "outputs": [],
   "execution_count": 4
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Filter and Preprocess the dataset",
   "id": "71f82a505f6aafd0"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:39:29.906549Z",
     "start_time": "2024-07-20T21:39:28.196002Z"
    }
   },
   "cell_type": "code",
   "source": [
    "#'HWTDGHTM' = height column\n",
    "#'HWTDGWTK' = weight column \n",
    "#'HWTDGBMI' = BMI column\n",
    "\n",
    "df = original_df \n",
    "# The initial size of the dataframe\n",
    "print(df.shape)\n",
    "\n",
    "# Filter the dataframe\n",
    "df['HWTDGHTM'] = pd.to_numeric(df['HWTDGHTM'], errors='coerce')\n",
    "df['HWTDGBMI'] = pd.to_numeric(df['HWTDGBMI'], errors='coerce')\n",
    "df['HWTDGWTK'] = pd.to_numeric(df['HWTDGWTK'], errors='coerce')\n",
    "\n",
    "# Filter the dataframe\n",
    "# Get rid of people who haven't answered\n",
    "df = df[(df['HWTDGHTM'] < 2.5) & (df['HWTDGBMI'] < 60) & (df['HWTDGWTK'] < 150)]\n",
    "\n",
    "# Filter so that you only have women in the dataframe\n",
    "df = df[df['DHH_SEX'] == 2]\n",
    "\n",
    "# Delete patients who are missing the chronic illness survey done \n",
    "df = df[df['DOCCC'] != 2]\n",
    "\n",
    "# Delete patients who are younger than 18 years old\n",
    "df = df[df['DHHGAGE'] > 2]\n",
    "\n",
    "# Now, we will further filter the DataFrame to remove rows where 'CCC_045' is 7, 8, or 9.\n",
    "# We only include rows that answered: \"I have fibromyalgia\" or \"I don't have fibromyalgia\"\n",
    "df = df[df['CCC_045'].isin([1, 2])]\n",
    "\n",
    "print(df.shape)\n",
    "unique_values = df['CCC_045'].value_counts()\n",
    "print(unique_values)"
   ],
   "id": "46dcc4dcd4104637",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(109659, 1284)\n",
      "(50028, 1284)\n",
      "CCC_045\n",
      "2.0    48133\n",
      "1.0     1895\n",
      "Name: count, dtype: int64\n"
     ]
    }
   ],
   "execution_count": 5
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Our dataset has 109659 rows (patients) and 1284 columns (variables) before any filtering. After filtering men and subjects with unacceptable answers out, we are left with 53596. Out of 1895 of them have fibromyalgia.\n",
    "\n",
    "Now, we will randomly select 1895 out of non-fibromyalgia subjects so that our statistical tests doesn't get skewed by the sheer size of the non-fibro sample."
   ],
   "id": "68a87d0145e6882f"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:39:30.206177Z",
     "start_time": "2024-07-20T21:39:29.910577Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# Filter the DataFrame to get rows where CCC_045 is 1.0 (yes fibromyalgia)\n",
    "df_1 = df[df['CCC_045'] == 1.0]\n",
    "\n",
    "# Filter the DataFrame to get rows where CCC_045 is 2.0 (no fibromyalgia)\n",
    "df_2 = df[df['CCC_045'] == 2.0]\n",
    "\n",
    "# Take a random sample of size 1897 from the 2.0 category\n",
    "sample_df_2 = df_2.sample(n=1895, random_state=42)\n",
    "\n",
    "# Combine the two DataFrames\n",
    "filtered_df = pd.concat([df_1, sample_df_2])\n",
    "\n",
    "# Replace values, now CCC_045 means that the patient doesn't have fibromyalgia\n",
    "filtered_df['CCC_045'] = filtered_df['CCC_045'].replace(2, 0)\n",
    "\n",
    "df = filtered_df\n",
    "\n",
    "unique_values = filtered_df['CCC_045'].value_counts()\n",
    "print(unique_values)\n",
    "\n",
    "unique_values = filtered_df['DHHGAGE'].value_counts()\n",
    "print(unique_values)"
   ],
   "id": "bde7e57d5cc52374",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CCC_045\n",
      "1.0    1895\n",
      "0.0    1895\n",
      "Name: count, dtype: int64\n",
      "DHHGAGE\n",
      "13.0    449\n",
      "11.0    444\n",
      "12.0    431\n",
      "10.0    383\n",
      "14.0    379\n",
      "16.0    280\n",
      "9.0     261\n",
      "8.0     236\n",
      "15.0    234\n",
      "7.0     204\n",
      "6.0     191\n",
      "5.0     151\n",
      "4.0     106\n",
      "3.0      41\n",
      "Name: count, dtype: int64\n"
     ]
    }
   ],
   "execution_count": 6
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Now, we have 1895 patients from each category (people who have or people who doesn't have fibromyalgia)\n",
    "\n",
    "Let's add an estimated skin area column to these subjects by using Du Bois formula: BSA(m^2)=0.007184 × ((Height(cm))^0.725) × ((Weight(kg))^0.425)\n",
    "\n",
    "BSA = Body Surface Area"
   ],
   "id": "1e75e6bb98977073"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:39:30.216487Z",
     "start_time": "2024-07-20T21:39:30.206529Z"
    }
   },
   "cell_type": "code",
   "source": "df['BSA'] = 0.007184 * (df['HWTDGHTM'] * 100) ** 0.725 * df['HWTDGWTK'] ** 0.425",
   "id": "494dd06d8f5b08fb",
   "outputs": [],
   "execution_count": 7
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Check correlations of height, weight, bmi and skin area (bsa) with fibromyalgia",
   "id": "788f935cec9e406f"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T12:39:04.241183Z",
     "start_time": "2024-07-15T12:39:04.195365Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# Columns\n",
    "height = df['HWTDGHTM']\n",
    "weight = df['HWTDGWTK']\n",
    "bmi = df['HWTDGBMI']\n",
    "bsa = df['BSA']\n",
    "fibromyalgia_status = df['CCC_045']\n",
    "\n",
    "###########################    HEIGHT    #################################\n",
    "# Pearson correlation\n",
    "pearson_corr, pearson_p = pearsonr(height, fibromyalgia_status)\n",
    "print(f'Height: Pearson correlation: {pearson_corr}, p-value: {pearson_p}')\n",
    "\n",
    "# Spearman correlation\n",
    "spearman_corr, spearman_p = spearmanr(height, fibromyalgia_status)\n",
    "print(f'Spearman correlation: {spearman_corr}, p-value: {spearman_p}')\n",
    "\n",
    "# Point-Biserial correlation (only if fibromyalgia_status is binary)\n",
    "if df['CCC_045'].nunique() == 2:\n",
    "    pointbiserial_corr, pointbiserial_p = pointbiserialr(fibromyalgia_status, height)\n",
    "    print(f'Point-Biserial correlation: {pointbiserial_corr}, p-value: {pointbiserial_p} \\n')\n",
    "else:\n",
    "    print('Point-Biserial correlation is not applicable as the fibromyalgia status variable is not binary.')\n",
    "\n",
    "###########################    WEIGHT    #################################\n",
    "# Pearson correlation\n",
    "pearson_corr, pearson_p = pearsonr(weight, fibromyalgia_status)\n",
    "print(f'WEIGHT: Pearson correlation: {pearson_corr}, p-value: {pearson_p}')\n",
    "\n",
    "# Spearman correlation\n",
    "spearman_corr, spearman_p = spearmanr(weight, fibromyalgia_status)\n",
    "print(f'Spearman correlation: {spearman_corr}, p-value: {spearman_p}')\n",
    "\n",
    "# Point-Biserial correlation (only if fibromyalgia_status is binary)\n",
    "if df['CCC_045'].nunique() == 2:\n",
    "    pointbiserial_corr, pointbiserial_p = pointbiserialr(fibromyalgia_status, weight)\n",
    "    print(f'Point-Biserial correlation: {pointbiserial_corr}, p-value: {pointbiserial_p} \\n')\n",
    "else:\n",
    "    print('Point-Biserial correlation is not applicable as the fibromyalgia status variable is not binary.')\n",
    "\n",
    "###########################    BMI    #################################\n",
    "# Pearson correlation\n",
    "pearson_corr, pearson_p = pearsonr(bmi, fibromyalgia_status)\n",
    "print(f'BMI: Pearson correlation: {pearson_corr}, p-value: {pearson_p}')\n",
    "\n",
    "# Spearman correlation\n",
    "spearman_corr, spearman_p = spearmanr(bmi, fibromyalgia_status)\n",
    "print(f'Spearman correlation: {spearman_corr}, p-value: {spearman_p}')\n",
    "\n",
    "# Point-Biserial correlation (only if fibromyalgia_status is binary)\n",
    "if df['CCC_045'].nunique() == 2:\n",
    "    pointbiserial_corr, pointbiserial_p = pointbiserialr(fibromyalgia_status, bmi)\n",
    "    print(f'Point-Biserial correlation: {pointbiserial_corr}, p-value: {pointbiserial_p} \\n')\n",
    "else:\n",
    "    print('Point-Biserial correlation is not applicable as the fibromyalgia status variable is not binary.')\n",
    "\n",
    "\n",
    "###########################    SURFACE AREA    #################################\n",
    "# Pearson correlation\n",
    "pearson_corr, pearson_p = pearsonr(bsa, fibromyalgia_status)\n",
    "print(f'BSA - ESTIMATED SKIN SURFACE AREA:\\nPearson correlation: {pearson_corr}, p-value: {pearson_p}')\n",
    "\n",
    "# Spearman correlation\n",
    "spearman_corr, spearman_p = spearmanr(bsa, fibromyalgia_status)\n",
    "print(f'Spearman correlation: {spearman_corr}, p-value: {spearman_p}')\n",
    "\n",
    "# Point-Biserial correlation (only if fibromyalgia_status is binary)\n",
    "if df['CCC_045'].nunique() == 2:\n",
    "    pointbiserial_corr, pointbiserial_p = pointbiserialr(fibromyalgia_status, bsa)\n",
    "    print(f'Point-Biserial correlation: {pointbiserial_corr}, p-value: {pointbiserial_p}')\n",
    "else:\n",
    "    print('Point-Biserial correlation is not applicable as the fibromyalgia status variable is not binary.')"
   ],
   "id": "b9805d82fd85376f",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Height: Pearson correlation: -0.048124182988062064, p-value: 0.00304253391776734\n",
      "Spearman correlation: -0.04408040980270724, p-value: 0.006644702049321087\n",
      "Point-Biserial correlation: -0.048124182988062064, p-value: 0.00304253391776734 \n",
      "\n",
      "WEIGHT: Pearson correlation: 0.1659135101184542, p-value: 8.486000575701152e-25\n",
      "Spearman correlation: 0.1671962750664145, p-value: 3.663633623736912e-25\n",
      "Point-Biserial correlation: 0.1659135101184542, p-value: 8.486000575701152e-25 \n",
      "\n",
      "BMI: Pearson correlation: 0.1875147513092818, p-value: 2.4740052245336616e-31\n",
      "Spearman correlation: 0.1898266581088423, p-value: 4.4064604150301994e-32\n",
      "Point-Biserial correlation: 0.1875147513092818, p-value: 2.4740052245336616e-31 \n",
      "\n",
      "BSA - ESTIMATED SKIN SURFACE AREA:\n",
      "Pearson correlation: 0.13278892713913593, p-value: 2.238592544270953e-16\n",
      "Spearman correlation: 0.13419980291509626, p-value: 1.0716276526950227e-16\n",
      "Point-Biserial correlation: 0.13278892713913593, p-value: 2.238592544270953e-16\n"
     ]
    }
   ],
   "execution_count": 17
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "##### All of the variables have weak (height being VERY weak) but statistically significant correlation to fibromyalgia (bmi - weight - bsa - height in order of strength high to low) \n",
    "\n",
    "Let's check the mean height/weight/bmi/bsa in fibro and non-fibro subjects"
   ],
   "id": "88da05c29ef2c999"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "#### Partial Correlations controlled for other Independent Variables",
   "id": "674adb9d19d8c67"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:40:23.396572Z",
     "start_time": "2024-07-20T21:40:23.193064Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import pandas as pd\n",
    "from pingouin import partial_corr\n",
    "\n",
    "# Load your data into a DataFrame (assuming the DataFrame is already loaded as df)\n",
    "# df = pd.read_csv('your_data.csv')  # Modify this line as needed to load your data\n",
    "\n",
    "# Assign columns to variables\n",
    "height = df['HWTDGHTM']\n",
    "weight = df['HWTDGWTK']\n",
    "bmi = df['HWTDGBMI']\n",
    "bsa = df['BSA']\n",
    "fibromyalgia_status = df['CCC_045']\n",
    "\n",
    "# Partial correlation between height and fibromyalgia status, controlling for weight, skin area, and BMI\n",
    "partial_corr_height = partial_corr(data=df, x='HWTDGHTM', y='CCC_045', covar=['HWTDGWTK', 'BSA', 'HWTDGBMI'])\n",
    "print(\"Partial correlation (height ~ fibromyalgia status | weight, skin area, BMI):\")\n",
    "print(partial_corr_height)\n",
    "\n",
    "# Partial correlation between weight and fibromyalgia status, controlling for height, skin area, and BMI\n",
    "partial_corr_weight = partial_corr(data=df, x='HWTDGWTK', y='CCC_045', covar=['HWTDGHTM', 'BSA', 'HWTDGBMI'])\n",
    "print(\"Partial correlation (weight ~ fibromyalgia status | height, skin area, BMI):\")\n",
    "print(partial_corr_weight)\n",
    "\n",
    "# Partial correlation between skin area (BSA) and fibromyalgia status, controlling for height, weight, and BMI\n",
    "partial_corr_bsa = partial_corr(data=df, x='BSA', y='CCC_045', covar=['HWTDGHTM', 'HWTDGWTK', 'HWTDGBMI'])\n",
    "print(\"Partial correlation (skin area ~ fibromyalgia status | height, weight, BMI):\")\n",
    "print(partial_corr_bsa)\n",
    "\n",
    "# Partial correlation between BMI and fibromyalgia status, controlling for height, weight, and skin area\n",
    "partial_corr_bmi = partial_corr(data=df, x='HWTDGBMI', y='CCC_045', covar=['HWTDGHTM', 'HWTDGWTK', 'BSA'])\n",
    "print(\"Partial correlation (BMI ~ fibromyalgia status | height, weight, skin area):\")\n",
    "print(partial_corr_bmi)\n"
   ],
   "id": "9b82899cead87bb9",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Partial correlation (height ~ fibromyalgia status | weight, skin area, BMI):\n",
      "            n         r          CI95%     p-val\n",
      "pearson  3790 -0.033423  [-0.07, -0.0]  0.039717\n",
      "Partial correlation (weight ~ fibromyalgia status | height, skin area, BMI):\n",
      "            n         r          CI95%     p-val\n",
      "pearson  3790 -0.006459  [-0.04, 0.03]  0.691119\n",
      "Partial correlation (skin area ~ fibromyalgia status | height, weight, BMI):\n",
      "            n         r          CI95%     p-val\n",
      "pearson  3790  0.023955  [-0.01, 0.06]  0.140512\n",
      "Partial correlation (BMI ~ fibromyalgia status | height, weight, skin area):\n",
      "            n         r          CI95%     p-val\n",
      "pearson  3790 -0.001249  [-0.03, 0.03]  0.938751\n"
     ]
    }
   ],
   "execution_count": 9
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "None of the independent variables (height, weight, skin area, BMI) show a strong or statistically significant direct relationship with fibromyalgia status after controlling for the other variables. The only significant result is the partial correlation between height and fibromyalgia status, but the effect size is very small and might not be practically significant.",
   "id": "248e1415428a1e40"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "#### Compare the means",
   "id": "40807a0583ca9287"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T12:44:00.773886Z",
     "start_time": "2024-07-15T12:44:00.741571Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# Define the columns\n",
    "height = df['HWTDGHTM']\n",
    "weight = df['HWTDGWTK']\n",
    "bmi = df['HWTDGBMI']\n",
    "bsa = df['BSA']\n",
    "age = df['DHHGAGE']\n",
    "fibromyalgia_status = df['CCC_045']\n",
    "\n",
    "# Separate the data into two groups\n",
    "fibro_group = df[df['CCC_045'] == 1]\n",
    "non_fibro_group = df[df['CCC_045'] == 0]\n",
    "\n",
    "# Calculate means\n",
    "mean_height_fibro = fibro_group['HWTDGHTM'].mean()\n",
    "mean_weight_fibro = fibro_group['HWTDGWTK'].mean()\n",
    "mean_BMI_fibro = fibro_group['HWTDGBMI'].mean()\n",
    "mean_BSA_fibro = fibro_group['BSA'].mean()\n",
    "mean_age_fibro = fibro_group['DHHGAGE'].mean()\n",
    "\n",
    "mean_height_non_fibro = non_fibro_group['HWTDGHTM'].mean()\n",
    "mean_weight_non_fibro = non_fibro_group['HWTDGWTK'].mean()\n",
    "mean_BMI_non_fibro = non_fibro_group['HWTDGBMI'].mean()\n",
    "mean_BSA_non_fibro = non_fibro_group['BSA'].mean()\n",
    "mean_age_non_fibro = non_fibro_group['DHHGAGE'].mean()\n",
    "\n",
    "\n",
    "# Perform t-tests\n",
    "t_height, p_height = stats.ttest_ind(fibro_group['HWTDGHTM'], non_fibro_group['HWTDGHTM'], equal_var=False)\n",
    "t_weight, p_weight = stats.ttest_ind(fibro_group['HWTDGWTK'], non_fibro_group['HWTDGWTK'], equal_var=False)\n",
    "t_BMI, p_BMI = stats.ttest_ind(fibro_group['HWTDGBMI'], non_fibro_group['HWTDGBMI'], equal_var=False)\n",
    "t_BSA, p_BSA = stats.ttest_ind(fibro_group['BSA'], non_fibro_group['BSA'], equal_var=False)\n",
    "t_age, p_age = stats.ttest_ind(fibro_group['DHHGAGE'], non_fibro_group['DHHGAGE'], equal_var=False)\n",
    "\n",
    "# Print the results\n",
    "print(\"Mean height of people with fibromyalgia:\", mean_height_fibro)\n",
    "print(\"Mean weight of people with fibromyalgia:\", mean_weight_fibro)\n",
    "print(\"Mean BMI of people with fibromyalgia:\", mean_BMI_fibro)\n",
    "print(\"Mean BSA of people with fibromyalgia:\", mean_BSA_fibro,)\n",
    "print(\"Mean age of people with fibromyalgia:\", mean_age_fibro, \"\\n\")\n",
    "\n",
    "print(\"Mean height of people without fibromyalgia:\", mean_height_non_fibro)\n",
    "print(\"Mean weight of people without fibromyalgia:\", mean_weight_non_fibro)\n",
    "print(\"Mean BMI of people without fibromyalgia:\", mean_BMI_non_fibro)\n",
    "print(\"Mean BSA of people without fibromyalgia:\", mean_BSA_non_fibro)\n",
    "print(\"Mean age of people without fibromyalgia:\", mean_age_non_fibro, \"\\n\")\n",
    "\n",
    "\n",
    "# Print the results\n",
    "print(f\"Height: t-statistic = {t_height:.4f}, p-value = {p_height:.4f}\")\n",
    "print(f\"Weight: t-statistic = {t_weight:.4f}, p-value = {p_weight:.4f}\")\n",
    "print(f\"BMI: t-statistic = {t_BMI:.4f}, p-value = {p_BMI:.4f}\")\n",
    "print(f\"BSA: t-statistic = {t_BSA:.4f}, p-value = {p_BSA:.4f}\")\n",
    "print(f\"Age: t-statistic = {t_age:.4f}, p-value = {p_age:.4f}\")\n"
   ],
   "id": "988b411cc488199c",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean height of people with fibromyalgia: 1.6163018469656991\n",
      "Mean weight of people with fibromyalgia: 73.93981530343008\n",
      "Mean BMI of people with fibromyalgia: 28.315540897097623\n",
      "Mean BSA of people with fibromyalgia: 1.7764724587004308\n",
      "Mean age of people with fibromyalgia: 11.624274406332454 \n",
      "\n",
      "Mean height of people without fibromyalgia: 1.6228559366754618\n",
      "Mean weight of people without fibromyalgia: 68.7397889182058\n",
      "Mean BMI of people without fibromyalgia: 26.117356200527706\n",
      "Mean BSA of people without fibromyalgia: 1.7283371064681714\n",
      "Mean age of people without fibromyalgia: 10.081266490765172 \n",
      "\n",
      "Height: t-statistic = -2.9653, p-value = 0.0030\n",
      "Weight: t-statistic = 10.3550, p-value = 0.0000\n",
      "BMI: t-statistic = 11.7493, p-value = 0.0000\n",
      "BSA: t-statistic = 8.2457, p-value = 0.0000\n",
      "Age: t-statistic = 14.8236, p-value = 0.0000\n"
     ]
    }
   ],
   "execution_count": 20
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "People who have fibromyalgia are, in average, older, shorter, heavier, has a higher BMI and skin area than those who don't have fibromyalgia. These differences are statistically significant as shown in the t-test.\n",
    "\n",
    "Age is divided into categories in our dataset: category 10 -> 50-54 , category 11 -> 55-59"
   ],
   "id": "6afc86d9128408e8"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "It is established that BMI and age positively correlate with fibromyalgia, in the following section we will take our new variables like heihgt, weight, skin area and see if they have an effect on the condition independent from BMI and Age",
   "id": "9a00c69f6100f254"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "#### Height and BMI as independent predictors of fibromylagia\n",
    "We will use logistic regression to see if height and bmi are independently predicting fibromyalgia, suggesting if each has its own effect on the disease"
   ],
   "id": "baff62e18209b5b2"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-14T23:02:25.075996Z",
     "start_time": "2024-07-14T23:02:25.046299Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# Only height and bmi model\n",
    "X = df[['HWTDGHTM', 'HWTDGBMI']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)\n"
   ],
   "id": "eae6750942b7973e",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.665031\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3794\n",
      "Model:                          Logit   Df Residuals:                     3791\n",
      "Method:                           MLE   Df Model:                            2\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.04056\n",
      "Time:                        02:02:25   Log-Likelihood:                -2523.1\n",
      "converged:                       True   LL-Null:                       -2629.8\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.697e-47\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.2758      0.821     -0.336      0.737      -1.885       1.333\n",
      "HWTDGHTM      -1.1710      0.487     -2.404      0.016      -2.126      -0.216\n",
      "HWTDGBMI       0.0809      0.006     13.457      0.000       0.069       0.093\n",
      "==============================================================================\n",
      "    Feature         VIF\n",
      "0     const  605.269428\n",
      "1  HWTDGHTM    1.014349\n",
      "2  HWTDGBMI    1.014349\n"
     ]
    }
   ],
   "execution_count": 33
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "(?) Coefficient of height shows that it is a better predictor than bmi, while both being statistically significant. \n",
    "\n",
    "The odds ratio represents the change in odds of the dependent variable (the likelihood of having fibromyalgia) for a one-unit increase in the predictor variable, holding all other variables constant.\n",
    "\n",
    "Height (HWTDGHTM) (per meters):\n",
    "Coefficient per meter: -1.1710\n",
    "Odds Ratio: \n",
    "e^−1.1710 ≈ 0.310\n",
    "\n",
    "Coefficient per cm = −1.1710/100 = −0.01171\n",
    "e^-0.01171 ≈ 0.988\n",
    "\n",
    "BMI (HWTDGBMI):\n",
    "Coefficient: 0.0809\n",
    "Odds Ratio: \n",
    "e^0.0809 ≈ 1.084\n",
    "\n",
    "Interpretation:\n",
    "\n",
    "For height, the odds ratio of 0.310 indicates that for each unit (per meters) increase in height (in meters), the odds of having fibromyalgia decrease by approximately 69% (since 1 - 0.310 = 0.69).\n",
    "For each centimeter increase in height, the odds of having fibromyalgia decrease by about 1.2% (since 1 - 0.988 = 0.012 or 1.2%).\n",
    "\n",
    "For BMI, the odds ratio of 1.084 suggests that for each unit increase in BMI, the odds of having fibromyalgia increase by about 8.4%.\n",
    "\n",
    "The VIF (Variance Inflation Factor) values for both height and BMI are around 1.014, indicating that multicollinearity is not a concern in this model.\n",
    "\n",
    "(?) It is suprising how height has weaker correlation but is a better predictor of fibromyalgia"
   ],
   "id": "31659a7b4b1236df"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T12:48:12.959963Z",
     "start_time": "2024-07-15T12:48:12.931872Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# height, bmi, age\n",
    "X = df[['HWTDGHTM', 'HWTDGBMI', 'DHHGAGE']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "d45f0079632b3cd4",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.648490\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3786\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.06443\n",
      "Time:                        15:48:12   Log-Likelihood:                -2457.8\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.612e-73\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -5.0899      0.914     -5.566      0.000      -6.882      -3.298\n",
      "HWTDGHTM       1.0147      0.519      1.954      0.051      -0.003       2.032\n",
      "HWTDGBMI       0.0660      0.006     10.918      0.000       0.054       0.078\n",
      "DHHGAGE        0.1516      0.011     13.670      0.000       0.130       0.173\n",
      "==============================================================================\n",
      "    Feature         VIF\n",
      "0     const  712.024753\n",
      "1  HWTDGHTM    1.078595\n",
      "2  HWTDGBMI    1.019577\n",
      "3   DHHGAGE    1.062905\n"
     ]
    }
   ],
   "execution_count": 24
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "##### It seems that hight looses it statistical significany (even though by a very close margin) when we include age. This might be due to the fact that people tends to lose some height as they grow older after a certain age because of the spine compression and posture problems.",
   "id": "365eaca9ebbf3d86"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Weight and BMI and Age as independent predictors of fibromyalgia",
   "id": "4a529c462d6f9642"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-14T23:28:20.621730Z",
     "start_time": "2024-07-14T23:28:20.598765Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['HWTDGWTK', 'HWTDGBMI']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)\n"
   ],
   "id": "3d9cd7ab1e762e5f",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.665019\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3794\n",
      "Model:                          Logit   Df Residuals:                     3791\n",
      "Method:                           MLE   Df Model:                            2\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.04058\n",
      "Time:                        02:28:20   Log-Likelihood:                -2523.1\n",
      "converged:                       True   LL-Null:                       -2629.8\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.494e-47\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -2.1714      0.164    -13.217      0.000      -2.493      -1.849\n",
      "HWTDGWTK      -0.0135      0.006     -2.422      0.015      -0.024      -0.003\n",
      "HWTDGBMI       0.1162      0.015      7.619      0.000       0.086       0.146\n",
      "==============================================================================\n",
      "    Feature        VIF\n",
      "0     const  21.846335\n",
      "1  HWTDGWTK   6.970374\n",
      "2  HWTDGBMI   6.970374\n"
     ]
    }
   ],
   "execution_count": 39
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Wight (HWTDGWTK) (per kilograms):\n",
    "Coefficient per meter: -0.0135\n",
    "Odds Ratio: \n",
    "e^−0.0135 ≈ 0.987\n",
    "\n",
    "\n",
    "BMI (HWTDGBMI):\n",
    "Coefficient:  0.1162\n",
    "Odds Ratio: \n",
    "e^ 0.1162 ≈ 1.123\n",
    "\n",
    "Interpretation:\n",
    "\n",
    "For weight, the odds ratio of 0.987 indicates that for each unit increase in weight (in kilograms), the odds of having fibromyalgia decrease by approximately 1.3% (since 1 - 0.987 = 0.013 or 1.3%).\n",
    "For BMI, the odds ratio of 1.123 suggests that for each unit increase in BMI, the odds of having fibromyalgia increase by about 12.3%.\n",
    "The VIF (Variance Inflation Factor) values for both weight and BMI are around 6.970, indicating some level of multicollinearity, but it is generally considered acceptable as it is below 10."
   ],
   "id": "8daa8304d2fe650d"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T12:55:53.902865Z",
     "start_time": "2024-07-15T12:55:53.851702Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['HWTDGWTK', 'HWTDGBMI', 'DHHGAGE']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "457571107a7b8c3e",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.648461\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3786\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.06447\n",
      "Time:                        15:55:53   Log-Likelihood:                -2457.7\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.122e-73\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -3.4490      0.211    -16.319      0.000      -3.863      -3.035\n",
      "HWTDGWTK       0.0117      0.006      2.011      0.044       0.000       0.023\n",
      "HWTDGBMI       0.0356      0.016      2.283      0.022       0.005       0.066\n",
      "DHHGAGE        0.1515      0.011     13.700      0.000       0.130       0.173\n",
      "==============================================================================\n",
      "    Feature        VIF\n",
      "0     const  33.586159\n",
      "1  HWTDGWTK   7.077416\n",
      "2  HWTDGBMI   7.098891\n",
      "3   DHHGAGE   1.055436\n"
     ]
    }
   ],
   "execution_count": 25
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "",
   "id": "e21b4f3800b487ed"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Weight (HWTDGWTK):\n",
    "Coefficient: 0.0117\n",
    "Odds Ratio: \n",
    "e^0.0117 ≈ 1.012\n",
    "\n",
    "BMI (HWTDGBMI):\n",
    "Coefficient: 0.0356\n",
    "Odds Ratio:\n",
    "e^0.0356 ≈ 1.036\n",
    " \n",
    "Age (DHHGAGE):\n",
    "Coefficient: 0.1515\n",
    "Odds Ratio: \n",
    "e^0.1515 ≈ 1.163\n",
    "\n",
    "\n",
    "Interpretation:\n",
    "\n",
    "For weight, the odds ratio of 1.012 indicates that for each unit increase in weight (in kilograms), the odds of having fibromyalgia increase by approximately 1.2%.\n",
    "For BMI, the odds ratio of 1.036 suggests that for each unit increase in BMI, the odds of having fibromyalgia increase by about 3.6%.\n",
    "For age, the odds ratio of 1.163 indicates that for each unit increase in age (in years), the odds of having fibromyalgia increase by approximately 16.3%.\n",
    "The VIF (Variance Inflation Factor) values indicate that there is some multicollinearity between weight and BMI (with VIFs around 7), but it is generally considered acceptable as it is below 10. Age shows a low VIF, indicating low multicollinearity with other variables.\n",
    "\n",
    "In our case a unit increase in age is as follows: \n",
    "Age between 18 and 19 - 03 \n",
    "Age between 20 and 24 - 04 \n",
    "Age between 25 and 29 - 05\n",
    "Age between 30 and 34 - 06\n",
    "Age between 35 and 39 - 07 \n",
    "Age between 40 and 44 - 08 \n",
    "Age between 45 and 49 - 09 \n",
    "Age between 50 and 54 - 10 \n",
    "Age between 55 and 59 - 11 \n",
    "Age between 60 and 64 - 12 \n",
    "Age between 65 and 69 - 13 \n",
    "Age between 70 and 74 - 14 \n",
    "Age between 75 and 79 - 15 \n",
    "Age 80 and older      - 16\n",
    "\n",
    "##### The fact that predictive power of weight stays stable but BMI lessens when we include age might be due to the fact that people grow shorter after a certain age which affects BMI."
   ],
   "id": "eef49ee347d8334d"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Weight and Skin Area and Age as independent predictors of fibromyalgia",
   "id": "35e021cd697ef0be"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T15:05:54.481996Z",
     "start_time": "2024-07-15T15:05:54.399149Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['HWTDGWTK', 'BSA']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "d4b9a2b6587aac56",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.675323\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3787\n",
      "Method:                           MLE   Df Model:                            2\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.02571\n",
      "Time:                        18:05:54   Log-Likelihood:                -2559.5\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.597e-30\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.6185      0.611      2.651      0.008       0.422       2.815\n",
      "HWTDGWTK       0.0588      0.007      8.078      0.000       0.045       0.073\n",
      "BSA           -3.3153      0.619     -5.359      0.000      -4.528      -2.103\n",
      "==============================================================================\n",
      "    Feature         VIF\n",
      "0     const  337.109308\n",
      "1  HWTDGWTK   11.119501\n",
      "2       BSA   11.119501\n"
     ]
    }
   ],
   "execution_count": 26
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Since there is significant multicollinearity, our coefficient estimates are not reliable.",
   "id": "b7605c7e7c6b7cfc"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T15:10:49.076039Z",
     "start_time": "2024-07-15T15:10:48.982831Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['HWTDGWTK', 'BSA', \"DHHGAGE\"]]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "489512a5829291c3",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.648472\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3786\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.06445\n",
      "Time:                        18:10:49   Log-Likelihood:                -2457.7\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.303e-73\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.9527      0.677     -2.883      0.004      -3.280      -0.625\n",
      "HWTDGWTK       0.0401      0.007      5.349      0.000       0.025       0.055\n",
      "BSA           -1.4603      0.645     -2.265      0.024      -2.724      -0.196\n",
      "DHHGAGE        0.1518      0.011     13.777      0.000       0.130       0.173\n",
      "==============================================================================\n",
      "    Feature         VIF\n",
      "0     const  390.481259\n",
      "1  HWTDGWTK   11.545639\n",
      "2       BSA   11.629896\n",
      "3   DHHGAGE    1.046727\n"
     ]
    }
   ],
   "execution_count": 30
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T15:09:09.322605Z",
     "start_time": "2024-07-15T15:09:09.284079Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['HWTDGWTK', 'BSA']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "# Adding the interaction term to the model\n",
    "X['weight_BSA_interaction'] = X['HWTDGWTK'] * X['BSA']\n",
    "\n",
    "interaction_model = sm.Logit(y, X).fit()\n",
    "print(interaction_model.summary())\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "29bfb7dab6a34bfd",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.674861\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3786\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.02638\n",
      "Time:                        18:09:09   Log-Likelihood:                -2557.7\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 7.541e-30\n",
      "==========================================================================================\n",
      "                             coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------------\n",
      "const                     -0.4567      1.264     -0.361      0.718      -2.934       2.020\n",
      "HWTDGWTK                   0.0956      0.021      4.572      0.000       0.055       0.137\n",
      "BSA                       -2.3105      0.818     -2.826      0.005      -3.913      -0.708\n",
      "weight_BSA_interaction    -0.0181      0.010     -1.879      0.060      -0.037       0.001\n",
      "==========================================================================================\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.674861\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3786\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.02638\n",
      "Time:                        18:09:09   Log-Likelihood:                -2557.7\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 7.541e-30\n",
      "==========================================================================================\n",
      "                             coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------------\n",
      "const                     -0.4567      1.264     -0.361      0.718      -2.934       2.020\n",
      "HWTDGWTK                   0.0956      0.021      4.572      0.000       0.055       0.137\n",
      "BSA                       -2.3105      0.818     -2.826      0.005      -3.913      -0.708\n",
      "weight_BSA_interaction    -0.0181      0.010     -1.879      0.060      -0.037       0.001\n",
      "==========================================================================================\n",
      "                  Feature          VIF\n",
      "0                   const  1404.213797\n",
      "1                HWTDGWTK    94.422511\n",
      "2                     BSA    19.165059\n",
      "3  weight_BSA_interaction   141.742694\n"
     ]
    }
   ],
   "execution_count": 29
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Also didn't work\n",
    "\n",
    "##### Didn't achieve any significant result from weight x skin area x age"
   ],
   "id": "83959bbd289fcc04"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Skin Area and BMI and Age as independent predictors of fibromyalgia",
   "id": "ce6ebe68e06ced16"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T15:11:20.754320Z",
     "start_time": "2024-07-15T15:11:20.722702Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['BSA', 'HWTDGBMI']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "17d8d553ae77257e",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.674980\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3787\n",
      "Method:                           MLE   Df Model:                            2\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.02621\n",
      "Time:                        18:11:20   Log-Likelihood:                -2558.2\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 1.250e-30\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.4373      0.354     -4.065      0.000      -2.130      -0.744\n",
      "BSA           -0.3466      0.287     -1.207      0.227      -0.909       0.216\n",
      "HWTDGBMI       0.0753      0.009      8.224      0.000       0.057       0.093\n",
      "==============================================================================\n",
      "    Feature         VIF\n",
      "0     const  113.717011\n",
      "1       BSA    2.454501\n",
      "2  HWTDGBMI    2.454501\n"
     ]
    }
   ],
   "execution_count": 31
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "Effect of BSA is not statistically significant\n",
    "\n",
    "This might be a counter-argument for the small fiber neuropathy since individual with larger skin areas have more small fibers. - > Not really: The hypothesis suggests that fibromyalgia patients may have a reduced density of small nerve fibers or dysfunctional small fibers, leading to pain and other sensory symptoms. This reduction in density or dysfunction can occur regardless of the total skin area.\n",
    "\n"
   ],
   "id": "a806fba701b5751d"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-15T15:13:23.048310Z",
     "start_time": "2024-07-15T15:13:23.016841Z"
    }
   },
   "cell_type": "code",
   "source": [
    "X = df[['BSA', 'HWTDGBMI','DHHGAGE']]\n",
    "y = df['CCC_045']\n",
    "\n",
    "# LOGISTIC REGRESSION\n",
    "# Adding a constant for the intercept\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "model = sm.Logit(y, X).fit()\n",
    "print(model.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "25e88cf182cc304a",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.648492\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3786\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Mon, 15 Jul 2024   Pseudo R-squ.:                 0.06442\n",
      "Time:                        18:13:23   Log-Likelihood:                -2457.8\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.640e-73\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -4.0484      0.413     -9.792      0.000      -4.859      -3.238\n",
      "BSA            0.5894      0.302      1.951      0.051      -0.003       1.181\n",
      "HWTDGBMI       0.0503      0.009      5.334      0.000       0.032       0.069\n",
      "DHHGAGE        0.1512      0.011     13.692      0.000       0.130       0.173\n",
      "==============================================================================\n",
      "    Feature         VIF\n",
      "0     const  140.650282\n",
      "1       BSA    2.573752\n",
      "2  HWTDGBMI    2.562858\n",
      "3   DHHGAGE    1.052597\n"
     ]
    }
   ],
   "execution_count": 34
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Effect of BSA is not statistically significant",
   "id": "c05229a4537ec37d"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-07-20T21:47:49.684132Z",
     "start_time": "2024-07-20T21:47:49.632296Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# Define the independent variables\n",
    "X = df[['HWTDGHTM', 'HWTDGWTK', 'HWTDGBMI', 'BSA', 'DHHGAGE']]\n",
    "X = sm.add_constant(X)  # Adds a constant term to the model\n",
    "\n",
    "# Define the dependent variable\n",
    "y = df['CCC_045']\n",
    "\n",
    "# Fit the logistic regression model\n",
    "logit_model = sm.Logit(y, X)\n",
    "result = logit_model.fit()\n",
    "\n",
    "# Print the summary of the logistic regression model\n",
    "print(result.summary())\n",
    "\n",
    "# Calculating VIF for each predictor\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data['Feature'] = X.columns\n",
    "vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]\n",
    "print(vif_data)"
   ],
   "id": "9c7c7de652f7fe1b",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.648444\n",
      "         Iterations 5\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                CCC_045   No. Observations:                 3790\n",
      "Model:                          Logit   Df Residuals:                     3784\n",
      "Method:                           MLE   Df Model:                            5\n",
      "Date:                Sun, 21 Jul 2024   Pseudo R-squ.:                 0.06449\n",
      "Time:                        00:47:49   Log-Likelihood:                -2457.6\n",
      "converged:                       True   LL-Null:                       -2627.0\n",
      "Covariance Type:            nonrobust   LLR p-value:                 4.391e-71\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -3.5522      4.022     -0.883      0.377     -11.435       4.331\n",
      "HWTDGHTM       0.7506      3.430      0.219      0.827      -5.972       7.473\n",
      "HWTDGWTK       0.0242      0.042      0.581      0.561      -0.057       0.106\n",
      "HWTDGBMI       0.0319      0.071      0.449      0.653      -0.107       0.171\n",
      "BSA           -1.0886      3.053     -0.357      0.721      -7.071       4.894\n",
      "DHHGAGE        0.1518      0.011     13.605      0.000       0.130       0.174\n",
      "==============================================================================\n",
      "    Feature           VIF\n",
      "0     const  12917.835008\n",
      "1  HWTDGHTM     42.800316\n",
      "2  HWTDGWTK    356.385239\n",
      "3  HWTDGBMI    138.933572\n",
      "4       BSA    249.173880\n",
      "5   DHHGAGE      1.082879\n"
     ]
    }
   ],
   "execution_count": 11
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}