stress_at_work_analysis/exploration/ex_ml_pipeline.py

# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py:percent
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.13.0
#   kernelspec:
#     display_name: straw2analysis
#     language: python
#     name: straw2analysis
# ---

# %%
# %matplotlib inline
import datetime
import importlib
import os
import sys

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import yaml
from pyprojroot import here
from sklearn import linear_model
from sklearn.model_selection import LeaveOneGroupOut, cross_val_score
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.impute import SimpleImputer

nb_dir = os.path.split(os.getcwd())[0]
if nb_dir not in sys.path:
    sys.path.append(nb_dir)

import machine_learning.features_sensor
import machine_learning.labels
import machine_learning.model

# %%
import participants.query_db
from features import esm, helper, proximity

# %% [markdown] tags=[]
# # 1. Get the relevant data

# %%
participants_inactive_usernames = participants.query_db.get_usernames(
    collection_start=datetime.date.fromisoformat("2020-08-01")
)
# Consider only two participants to simplify.
ptcp_2 = participants_inactive_usernames[0:2]

# %% [markdown] jp-MarkdownHeadingCollapsed=true tags=[]
# ## 1.1 Labels

# %%
df_esm = esm.get_esm_data(ptcp_2)
df_esm_preprocessed = esm.preprocess_esm(df_esm)

# %%
df_esm_PANAS = df_esm_preprocessed[
    (df_esm_preprocessed["questionnaire_id"] == 8)
    | (df_esm_preprocessed["questionnaire_id"] == 9)
]
df_esm_PANAS_clean = esm.clean_up_esm(df_esm_PANAS)

# %% [markdown]
# ## 1.2 Sensor data

# %%
df_proximity = proximity.get_proximity_data(ptcp_2)
df_proximity = helper.get_date_from_timestamp(df_proximity)
df_proximity = proximity.recode_proximity(df_proximity)

# %% [markdown]
# ## 1.3 Standardization/personalization

# %% [markdown]
# # 2. Grouping/segmentation

# %%
df_esm_PANAS_daily_means = (
    df_esm_PANAS_clean.groupby(["participant_id", "date_lj", "questionnaire_id"])
    .esm_user_answer_numeric.agg("mean")
    .reset_index()
    .rename(columns={"esm_user_answer_numeric": "esm_numeric_mean"})
)

# %%
df_esm_PANAS_daily_means = (
    df_esm_PANAS_daily_means.pivot(
        index=["participant_id", "date_lj"],
        columns="questionnaire_id",
        values="esm_numeric_mean",
    )
    .reset_index(col_level=1)
    .rename(columns={8.0: "PA", 9.0: "NA"})
    .set_index(["participant_id", "date_lj"])
)


# %%
df_proximity_daily_counts = proximity.count_proximity(df_proximity, ["date_lj"])

# %%
df_proximity_daily_counts

# %% [markdown]
# # 3. Join features (and export to csv?)

# %%
df_full_data_daily_means = df_esm_PANAS_daily_means.join(
    df_proximity_daily_counts
).reset_index()

# %% [markdown]
# # 4. Machine learning model and parameters

# %%
lin_reg_proximity = linear_model.LinearRegression()

# %% [markdown]
# ## 4.1 Validation method

# %%
logo = LeaveOneGroupOut()
logo.get_n_splits(
    df_full_data_daily_means[["freq_prox_near", "prop_prox_near"]],
    df_full_data_daily_means["PA"],
    groups=df_full_data_daily_means["participant_id"],
)

# %% [markdown]
# ## 4.2 Fit results (export?)

# %%
cross_val_score(
    lin_reg_proximity,
    df_full_data_daily_means[["freq_prox_near", "prop_prox_near"]],
    df_full_data_daily_means["PA"],
    groups=df_full_data_daily_means["participant_id"],
    cv=logo,
    n_jobs=-1,
    scoring="r2",
)

# %%
lin_reg_proximity.fit(
    df_full_data_daily_means[["freq_prox_near", "prop_prox_near"]],
    df_full_data_daily_means["PA"],
)

# %%
lin_reg_proximity.score(
    df_full_data_daily_means[["freq_prox_near", "prop_prox_near"]],
    df_full_data_daily_means["PA"],
)

# %% [markdown]
# # Merging these into a pipeline

# %%
from machine_learning import features_sensor, labels, model, pipeline

# %%
importlib.reload(features_sensor)

# %%
with open("../machine_learning/config/minimal_features.yaml", "r") as file:
    sensor_features_params = yaml.safe_load(file)
print(sensor_features_params)

# %%
sensor_features = machine_learning.features_sensor.SensorFeatures(
    **sensor_features_params
)
sensor_features.data_types

# %%
sensor_features.set_participants_label("nokia_0000003")

# %%
sensor_features.data_types = ["proximity", "communication"]
sensor_features.participants_usernames = ptcp_2

# %%
sensor_features.get_sensor_data("proximity")

# %%
sensor_features.set_sensor_data()

# %%
sensor_features.get_sensor_data("proximity")

# %%
sensor_features.calculate_features(cached=False)
features_all_calculated = sensor_features.get_features("all", "all")

# %%
sensor_features.calculate_features(cached=True)
features_all_read = sensor_features.get_features("all", "all")

# %%
features_all_read = features_all_read.reset_index()
features_all_read["date_lj"] = features_all_read["date_lj"].dt.date
features_all_read.set_index(["participant_id", "date_lj"], inplace=True)
# date_lj column is parsed as a date and represented as Timestamp, when read from csv.
# When calculated, it is represented as date.

# %%
np.isclose(features_all_read, features_all_calculated).all()

# %%
with open("../machine_learning/config/minimal_labels.yaml", "r") as file:
    labels_params = yaml.safe_load(file)

# %%
labels = machine_learning.labels.Labels(**labels_params)
labels.participants_usernames = ptcp_2
labels.set_participants_label("nokia_0000003")
labels.questionnaires

# %%
labels.set_labels()

# %%
labels.get_labels("PANAS")

# %%
labels.aggregate_labels(cached=False)
labels_calculated = labels.get_aggregated_labels()

# %%
labels.aggregate_labels(cached=True)
labels_read = labels.get_aggregated_labels()
labels_read = labels_read.reset_index()
labels_read["date_lj"] = labels_read["date_lj"].dt.date
labels_read.set_index(["participant_id", "date_lj"], inplace=True)
# date_lj column is parsed as a date and represented as Timestamp, when read from csv.
# When calculated, it is represented as date.

# %%
np.isclose(labels_read, labels_calculated).all()

# %%
model_validation = machine_learning.model.ModelValidation(
    sensor_features.get_features("all", "all"),
    labels.get_aggregated_labels(),
    group_variable="participant_id",
    cv_name="loso",
)
model_validation.model = linear_model.LinearRegression()
model_validation.set_cv_method()

# %%
model_validation.cross_validate()

# %%
model_validation.groups

# %% [markdown]
# # Use RAPIDS

# %%
with open(here("rapids/config.yaml"), "r") as file:
    rapids_config = yaml.safe_load(file)

# %%
for key in rapids_config.keys():    
    if isinstance(rapids_config[key], dict):  # Remove top-level configs
        if ("PROVIDERS" in rapids_config[key]):  # Retain features (that have providers)
            if rapids_config[key]["PROVIDERS"]:  # Remove non-implemented features
                for provider in rapids_config[key]["PROVIDERS"]:
                    if rapids_config[key]["PROVIDERS"][provider]["COMPUTE"]:  # Check that the features were actually calculated
                        if "FEATURES" in rapids_config[key]["PROVIDERS"][provider]:
                            print(key)
                            print(provider)
                            print(rapids_config[key]["PROVIDERS"][provider]["FEATURES"])

# %%
features_rapids = pd.read_csv(here("rapids/data/processed/features/all_participants/all_sensor_features.csv"), parse_dates=["local_segment_start_datetime", "local_segment_end_datetime"])

# %%
features_rapids.columns

# %%
features_rapids = features_rapids.assign(date_lj=lambda x: x.local_segment_start_datetime.dt.date)

# %%
features_rapids["participant_id"] = features_rapids["pid"].str.extract("(\d+)")
features_rapids["participant_id"] = pd.to_numeric(features_rapids["participant_id"])
features_rapids.set_index(["participant_id", "date_lj"], inplace=True)

# %%
with open("../machine_learning/config/minimal_labels.yaml", "r") as file:
    labels_params = yaml.safe_load(file)

# %%
labels = machine_learning.labels.Labels(**labels_params)
labels.set_participants_label("all")

# %%
labels.aggregate_labels(cached=True)
labels_read = labels.get_aggregated_labels()
labels_read = labels_read.reset_index()
labels_read["date_lj"] = labels_read["date_lj"].dt.date
labels_read.set_index(["participant_id", "date_lj"], inplace=True)
# date_lj column is parsed as a date and represented as Timestamp, when read from csv.
# When calculated, it is represented as date.

# %%
features_rapids.shape

# %%
labels_read.shape

# %%
features_labels = features_rapids.join(labels_read, how="inner").reset_index()

# %%
features_labels.shape

# %%
features_labels.columns

# %%
imputer = SimpleImputer(missing_values=np.nan, strategy='mean')

# %%
feature_columns = features_labels.columns[6:-3]
label_column = "NA"
group_column = "pid"

# %%
lin_reg_rapids = linear_model.LinearRegression()
logo = LeaveOneGroupOut()
logo.get_n_splits(
    features_labels[feature_columns],
    features_labels[label_column],
    groups=features_labels[group_column],
)

# %%
cross_val_score(
    lin_reg_rapids,
    X=imputer.fit_transform(features_labels[feature_columns]),
    y=features_labels[label_column],
    groups=features_labels[group_column],
    cv=logo,
    n_jobs=-1,
    scoring="r2",
)

# %%
sns.set(rc={"figure.figsize":(16, 8)})
sns.heatmap(features_labels[feature_columns].isna(), cbar=False)

# %% [markdown] tags=[]
# ```yaml
# ALL_CLEANING_INDIVIDUAL:
#   PROVIDERS:
#     RAPIDS:
#       COMPUTE: True
#       IMPUTE_SELECTED_EVENT_FEATURES: # Fill NAs with 0 only for event-based features, see table below
#         COMPUTE: True
#         MIN_DATA_YIELDED_MINUTES_TO_IMPUTE: 0.33 # Any feature value in a time segment instance with phone data yield > [MIN_DATA_YIELDED_MINUTES_TO_IMPUTE] will be replaced with a zero. 
#       COLS_NAN_THRESHOLD: 0.3 # Discard columns with missing value ratios higher than [COLS_NAN_THRESHOLD]. Set to 1 to disable
#       COLS_VAR_THRESHOLD: True # Set to True to discard columns with zero variance
#       ROWS_NAN_THRESHOLD: 1 # Discard rows with missing value ratios higher than [ROWS_NAN_THRESHOLD]. Set to 1 to disable
#       DATA_YIELD_FEATURE: RATIO_VALID_YIELDED_HOURS # RATIO_VALID_YIELDED_HOURS or RATIO_VALID_YIELDED_MINUTES
#       DATA_YIELD_RATIO_THRESHOLD: 0.3 # Discard rows with ratiovalidyieldedhours or ratiovalidyieldedminutes feature less than [DATA_YIELD_RATIO_THRESHOLD]. The feature name is determined by [DATA_YIELD_FEATURE] parameter. Set to 0 to disable
#       DROP_HIGHLY_CORRELATED_FEATURES:
#         COMPUTE: False
#         MIN_OVERLAP_FOR_CORR_THRESHOLD: 0.5
#         CORR_THRESHOLD: 0.95
#       SRC_SCRIPT: src/features/all_cleaning_individual/rapids/main.R
# ```

# %%
features_rapids_cleaned = pd.read_csv(here("rapids/data/processed/features/all_participants/all_sensor_features_cleaned_rapids.csv"), parse_dates=["local_segment_start_datetime", "local_segment_end_datetime"])
features_rapids_cleaned = features_rapids_cleaned.assign(date_lj=lambda x: x.local_segment_start_datetime.dt.date)
features_rapids_cleaned["participant_id"] = features_rapids_cleaned["pid"].str.extract("(\d+)")
features_rapids_cleaned["participant_id"] = pd.to_numeric(features_rapids_cleaned["participant_id"])
features_rapids_cleaned.set_index(["participant_id", "date_lj"], inplace=True)

# %%
features_cleaned_labels = features_rapids_cleaned.join(labels_read, how="inner").reset_index()
feature_clean_columns = features_cleaned_labels.columns[6:-3]

# %%
print(feature_columns.shape)
print(feature_clean_columns.shape)

# %%
sns.set(rc={"figure.figsize":(16, 8)})
sns.heatmap(features_cleaned_labels[feature_clean_columns].isna(), cbar=False)

# %%
lin_reg_rapids_clean = linear_model.LinearRegression()
logo = LeaveOneGroupOut()
logo.get_n_splits(
    features_cleaned_labels[feature_clean_columns],
    features_cleaned_labels[label_column],
    groups=features_cleaned_labels[group_column],
)

# %%
features_clean_imputed = imputer.fit_transform(features_cleaned_labels[feature_clean_columns])

# %%
cross_val_score(
    lin_reg_rapids_clean,
    X=features_clean_imputed,
    y=features_cleaned_labels[label_column],
    groups=features_cleaned_labels[group_column],
    cv=logo,
    n_jobs=-1,
    scoring="r2",
)

# %%
lin_reg_full = linear_model.LinearRegression()
lin_reg_full.fit(features_clean_imputed,features_cleaned_labels[label_column])

# %%
NA_pred = lin_reg_full.predict(features_clean_imputed)

# %%
# The coefficients
print("Coefficients: \n", lin_reg_full.coef_)
# The mean squared error
print("Mean squared error: %.2f" % mean_squared_error(features_cleaned_labels[label_column], NA_pred))
# The coefficient of determination: 1 is perfect prediction
print("Coefficient of determination: %.2f" % r2_score(features_cleaned_labels[label_column], NA_pred))

# %%
feature_clean_columns[np.argmax(lin_reg_full.coef_)]

# %% [markdown]
# Ratio between stationary time and total location sensed time. A lat/long coordinate pair is labeled as stationary if its speed (distance/time) to the next coordinate pair is less than 1km/hr. A higher value represents a more stationary routine.

# %%
plt.scatter(features_clean_imputed[:,np.argmax(lin_reg_full.coef_)], features_cleaned_labels[label_column], color="black")
plt.scatter(features_clean_imputed[:,np.argmax(lin_reg_full.coef_)], NA_pred, color="red", linewidth=3)

plt.xticks()
plt.yticks()

fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)

plt.show()