
import os

tmp = os.getcwd()
_ = os.chdir(tmp.split("Coding/")[0] + "Coding/data-science-core-develop")

from neuropy.correlation import (correlation_analysis, 
                                 correlations_as_sample_increases, 
                                plot_correlogram)
from neuropy.utils import tiny_summary
from neuropy.anomaly_detection import univariate_outlier_removal
from neuropy.frequentist_statistics import (analysis_independent_t_test,
                                            analysis_independent_mannwhitneyu_test,
                                            normal_check, 
                                            find_optimal_transformation, 
                                            apply_power_transformations,
                                           multiple_ICCs,
                                           bland_altman_plot_repeated,
                                           diagnostic_plots)

_ = os.chdir(tmp)


import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import session_info

# from jmspack.NLTSA import (flatten, 
#                           ts_levels,
#                           fluctuation_intensity,
#                           distribution_uniformity,
#                           complexity_resonance,
#                           complexity_resonance_diagram)

from sklearn.preprocessing import MinMaxScaler, StandardScaler

from work_motivation_extras import *

import statsmodels.api as sm
import statsmodels.formula.api as smf

# Import the lmm model class
# from pymer4.models import Lmer
# from pymer4.stats import rsquared, rsquared_adj, vif

# To use this experimental feature, we need to explicitly ask for it:
# from sklearn.experimental import enable_iterative_imputer  # noqa
# from sklearn.impute import IterativeImputer
# from sklearn.tree import DecisionTreeRegressor
# from sklearn.ensemble import ExtraTreesRegressor
# from sklearn.neighbors import KNeighborsRegressor
# from sklearn.pipeline import make_pipeline

# from pyunicorn.timeseries import RecurrencePlot


from sklearn.metrics import mean_squared_error, mean_absolute_error


session_info.show(write_req_file=False,
                  req_file_name="work_motivation_pre_mid_post_EDA_requirements.txt",)

-----
matplotlib          3.3.4
neuropy             1.4.3
numpy               1.19.2
pandas              1.2.3
seaborn             0.11.1
session_info        1.0.0
sklearn             0.24.1
statsmodels         0.12.2
-----

PIL                         8.1.2
appnope                     0.1.2
backcall                    0.2.0
cffi                        1.14.5
colorama                    0.4.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.1
decorator                   4.4.2
ipykernel                   5.3.4
ipython_genutils            0.2.0
ipywidgets                  7.6.3
jedi                        0.17.2
joblib                      0.17.0
kiwisolver                  1.3.1
mpl_toolkits                NA
parso                       0.7.0
patsy                       0.5.1
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
prompt_toolkit              3.0.8
ptyprocess                  0.7.0
pyexpat                     NA
pygments                    2.8.1
pyparsing                   2.4.7
pytz                        2021.1
scipy                       1.5.3
six                         1.15.0
storemagic                  NA
tornado                     6.1
traitlets                   5.0.5
wcwidth                     0.2.5
work_motivation_extras      NA
zmq                         20.0.0

-----
IPython             7.21.0
jupyter_client      6.1.7
jupyter_core        4.7.1
jupyterlab          2.2.6
notebook            6.2.0
-----
Python 3.9.2 (default, Mar  3 2021, 11:58:52) [Clang 10.0.0 ]
macOS-10.16-x86_64-i386-64bit
-----
Session information updated at 2021-06-13 20:54


_df = pd.read_csv("data/moti_feasibility_james_pre-mid-post.csv")


_df = _df.assign(date=pd.to_datetime(_df["datestamp"], format="%Y-%m-%d %H:%M:%S"))


def convert_numeric_where_possible(x):
    try:
        return x.astype(np.number)
    except:
        return x


demographics_columns = ['User', 'TIME', 'date', 'Allocation', 'gender', 'age', 'education']


target_columns = ["Amotivation", "ExternalTotal", "Introjected", "Identified", "Intrinsic"]


df = (_df
      .drop("datestamp", axis=1)
      .apply(convert_numeric_where_possible)
      .loc[:, demographics_columns + sum_scores_columns]
      .assign(RAI = lambda x: x["Amotivation"] * -3 + x["ExternalTotal"] * -2 + x["Introjected"] * -1 + x["Identified"] * 2 + x["Intrinsic"] * 3)
     .drop(target_columns, axis=1))
df.head(2)


df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 234 entries, 0 to 233
Data columns (total 31 columns):
 #   Column                 Non-Null Count  Dtype         
---  ------                 --------------  -----         
 0   User                   234 non-null    object        
 1   TIME                   234 non-null    float64       
 2   date                   234 non-null    datetime64[ns]
 3   Allocation             228 non-null    float64       
 4   gender                 101 non-null    float64       
 5   age                    98 non-null     float64       
 6   education              101 non-null    object        
 7   ExternalSocial         214 non-null    float64       
 8   ExternalMaterial       214 non-null    float64       
 9   AutonomySat            211 non-null    float64       
 10  CompetenceSat          212 non-null    float64       
 11  RelatedSat             211 non-null    float64       
 12  AutonomyFrust          211 non-null    float64       
 13  CompetenceFrust        211 non-null    float64       
 14  RelatedFrust           211 non-null    float64       
 15  Vigor                  212 non-null    float64       
 16  Dedication             212 non-null    float64       
 17  Absorption             212 non-null    float64       
 18  WorkEngagement         212 non-null    float64       
 19  SelfControl            209 non-null    float64       
 20  Congruence             209 non-null    float64       
 21  Interest               209 non-null    float64       
 22  Control                209 non-null    float64       
 23  AutonomousFunctioning  209 non-null    float64       
 24  IncStructResources     183 non-null    float64       
 25  DecHinderDemands       183 non-null    float64       
 26  IncChallengeDemands    183 non-null    float64       
 27  IncSocialResources     183 non-null    float64       
 28  Resilience             183 non-null    float64       
 29  Optimism               183 non-null    float64       
 30  RAI                    214 non-null    float64       
dtypes: datetime64[ns](1), float64(28), object(2)
memory usage: 56.8+ KB


tiny_summary(df)

====== SUMMARY =====
- Number of columns = 31
- Number of samples = 234 

- Percentage % of missing values per feature
User                      0.00
TIME                      0.00
date                      0.00
Allocation                2.56
gender                   56.84
age                      58.12
education                56.84
ExternalSocial            8.55
ExternalMaterial          8.55
AutonomySat               9.83
CompetenceSat             9.40
RelatedSat                9.83
AutonomyFrust             9.83
CompetenceFrust           9.83
RelatedFrust              9.83
Vigor                     9.40
Dedication                9.40
Absorption                9.40
WorkEngagement            9.40
SelfControl              10.68
Congruence               10.68
Interest                 10.68
Control                  10.68
AutonomousFunctioning    10.68
IncStructResources       21.79
DecHinderDemands         21.79
IncChallengeDemands      21.79
IncSocialResources       21.79
Resilience               21.79
Optimism                 21.79
RAI                       8.55
dtype: float64

- Total number of missing values = 1124

- Overall percentage of missing values = 0.155 %
=================


row_amount_df = (df
.reset_index()
.groupby("User")
.count()
.loc[:, ["date"]]
.rename(columns={"date": "row_amount"})
.sort_values(by="row_amount")
.reset_index())


row_amount_df.tail(1)


_ = plt.figure(figsize=(20, 4))
_ = sns.barplot(data=row_amount_df, x="User", y="row_amount")
_ = plt.xticks(rotation=90)
_ = plt.axhline(3, c="red", ls="--", label="3 day cutoff")
_ = plt.title("Row amounts per user")
_ = plt.legend()


scale_data = True
plot_df=df.set_index("date").select_dtypes(np.number)
plot_df.index = plot_df.index.date
if scale_data:
    plot_df = pd.DataFrame(MinMaxScaler().fit_transform(plot_df), index=plot_df.index, columns=plot_df.columns)

_ = plt.figure(figsize=(20, 10))
_ = sns.heatmap(plot_df.T)
_ = plt.title(f"Heatmap of raw values - all users")


plot_df = (df
           .set_index("User")
           .select_dtypes(np.number)
           .reset_index()
           .melt(id_vars=["User", "TIME"]))


_ = plt.figure(figsize=(20, 10))
_ = sns.lineplot(data=plot_df, x="TIME", y="value", hue="User", legend=False)


# _ = plt.figure(figsize=(20, 10))
# _ = sns.barplot(data=plot_df, x="variable", y="value")
# _ = plt.xticks(rotation=90)


_ = plt.figure(figsize=(20, 10))
_ = sns.boxplot(data=plot_df, x="variable", y="value")
_ = plt.xticks(rotation=90)


plot_bool = True
df_RT_no_out, O_summary, _ = univariate_outlier_removal(df=df, 
                                                        nonparam_args={'factor': 3},
                                                        nonparam_name='inter-quartile range, k=3',
                                                        remove=False, #set remove to True to remove the outliers from the original data frame (as well as the returned data frame) 
                                                        plot=plot_bool,
                                                        figsize = (10, 10)
                                                       ) 
# _ = plt.savefig("images/Outliers_plot.png", dpi=400, format="png", bbox_inches='tight')


O_summary.sort_values(by = "outliers-removed %", ascending=False).head(5)


# subset data frame with only numerical variables
X = df.select_dtypes(np.number)
df_normal_check = normal_check(X)

# data frame containing variables with Gaussian distribution
X_normal    = X[df_normal_check.loc[df_normal_check['normality'] == True, 'feature']]

# data frame containing variables with non-Gaussian distribution
X_nonnormal = X[df_normal_check.loc[df_normal_check['normality'] == False, 'feature']]

# Print Info
display(df_normal_check.T)
print("====== Normality test based on Kolmogorov-Smirnov test =====\n")
print(f"- Size normally distributed feature: {X_normal.shape[1]}")
print(f"- Size NOT normally distributed feature: {X_nonnormal.shape[1]}")

====== Normality test based on Kolmogorov-Smirnov test =====

- Size normally distributed feature: 9
- Size NOT normally distributed feature: 19


transformed_X_nonnormal = apply_power_transformations(X[X_nonnormal.columns])

df_transformed_normal_check = normal_check(transformed_X_nonnormal)
use_transformation, use_nonparametric = find_optimal_transformation(df_transformed_normal_check, X_nonnormal.columns)

print(f"Normal features = {X_normal.columns.to_numpy()}")
print(f"Normal features after transformation = {np.array(use_transformation)}")
print(f"Non normal features after transformation = {np.array(use_nonparametric)}")

Normal features = ['age' 'ExternalSocial' 'ExternalMaterial' 'AutonomyFrust' 'Control'
 'AutonomousFunctioning' 'DecHinderDemands' 'IncSocialResources' 'RAI']
Normal features after transformation = ['box_square_AutonomySat' 'box_square_Interest'
 'box_square_IncStructResources' 'box_square_IncChallengeDemands']
Non normal features after transformation = ['TIME' 'Allocation' 'gender' 'CompetenceSat' 'RelatedSat'
 'CompetenceFrust' 'RelatedFrust' 'Vigor' 'Dedication' 'Absorption'
 'WorkEngagement' 'SelfControl' 'Congruence' 'Resilience' 'Optimism']


len(X_normal.columns), len(use_transformation), len(use_nonparametric), (len(X_normal.columns) + len(use_transformation) + len(use_nonparametric))

(9, 4, 15, 28)


trans_df = (pd.merge(X_normal, 
         transformed_X_nonnormal[use_transformation + use_nonparametric], 
         left_index=True, 
         right_index=True)
.merge(df[["date"]], left_index=True, right_index=True)
           )


trans_df["gender"].value_counts()

0.0    75
1.0    23
3.0     3
Name: gender, dtype: int64


trans_df = df[df["gender"] != 3]
# df["gender"] = df.gender.astype(int)


trans_df["gender"].value_counts()

0.0    75
1.0    23
Name: gender, dtype: int64


plot_densities = True
if plot_densities:
    # Consider only float features
    features_list = df.select_dtypes(float).drop("gender", axis=1).columns.tolist()
    # Consider the interested group to analyze
    grouping_var = "gender"

    # Plot
    set_plot_limits=False
    fig, axs = plt.subplots(figsize=(20,75), nrows=len(features_list), ncols=2, gridspec_kw={'width_ratios': [2, 1]})
    fig.subplots_adjust(hspace = 0.5, wspace=0.2)
    axs = axs.ravel()

    for i in range(0, len(features_list)*2, 2):
        try:
            _ = sns.histplot(data=trans_df, x=features_list[int(i/2)], hue=grouping_var, kde=False, bins=10, ax=axs[i])
            _ = sns.boxplot(data=trans_df, x=grouping_var, y=features_list[int(i/2)], ax=axs[i+1])
            _ = axs[i+1].axes.get_yaxis().get_label().set_visible(False)
        except:
            print(f"{features_list[int(i/2)]} - this had an error")

    # _ = plt.savefig(f"images/density_boxplots_all_fitbit_features_{grouping_var}.png", dpi=400, format="png", bbox_inches='tight')


transformed_fitbit_cols = X_normal.columns.tolist() + use_transformation + use_nonparametric
np.array(transformed_fitbit_cols)

array(['age', 'ExternalSocial', 'ExternalMaterial', 'AutonomyFrust',
       'Control', 'AutonomousFunctioning', 'DecHinderDemands',
       'IncSocialResources', 'RAI', 'box_square_AutonomySat',
       'box_square_Interest', 'box_square_IncStructResources',
       'box_square_IncChallengeDemands', 'TIME', 'Allocation', 'gender',
       'CompetenceSat', 'RelatedSat', 'CompetenceFrust', 'RelatedFrust',
       'Vigor', 'Dedication', 'Absorption', 'WorkEngagement',
       'SelfControl', 'Congruence', 'Resilience', 'Optimism'],
      dtype='<U30')


corr_out, _ = correlation_analysis(data=trans_df,
#                     col_list=["RAI", "Optimism"],
#                     row_list = transformed_fitbit_cols,
                                  check_norm=True)
corrs = corr_out["summary"]

/opt/miniconda3/envs/general/lib/python3.9/site-packages/scipy/stats/stats.py:4196: SpearmanRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.
  warnings.warn(SpearmanRConstantInputWarning())
/opt/miniconda3/envs/general/lib/python3.9/site-packages/scipy/stats/stats.py:4196: SpearmanRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.
  warnings.warn(SpearmanRConstantInputWarning())


corrs[corrs["feature2"] == "RAI"].sort_values("p-value").head(20)


fig = plot_correlogram(data=trans_df,
#                     col_list=transformed_fitbit_cols,
#                     row_list = ["RAI"],
                                  check_norm=True,
                    figsize=(20, 10),
                         font_scale=0.7)

/opt/miniconda3/envs/general/lib/python3.9/site-packages/scipy/stats/stats.py:4196: SpearmanRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.
  warnings.warn(SpearmanRConstantInputWarning())
/opt/miniconda3/envs/general/lib/python3.9/site-packages/scipy/stats/stats.py:4196: SpearmanRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.
  warnings.warn(SpearmanRConstantInputWarning())


top_feature1 = corrs[corrs["feature2"] == "RAI"].sort_values("r-value").head(1).feature1.tolist()[0]
top_feature2 = corrs[corrs["feature2"] == "RAI"].sort_values("r-value").head(1).feature2.tolist()[0]


_ = sns.regplot(x = top_feature2, y = top_feature1, data=trans_df)
_ = plt.title(f'{top_feature1} vs {top_feature2.replace("_", " ").title()} regression plot')
_ = sns.despine()


top_feature1 = corrs[corrs["feature2"] == "RAI"].sort_values("r-value").tail(1).feature1.tolist()[0]
top_feature2 = corrs[corrs["feature2"] == "RAI"].sort_values("r-value").tail(1).feature2.tolist()[0]


_ = sns.regplot(x = top_feature2, y = top_feature1, data=trans_df)
_ = plt.title(f'{top_feature1} vs {top_feature2.replace("_", " ").title()} regression plot')
_ = sns.despine()


amount_top_correlations = 2
top_corrs_sorted = corrs[corrs["feature2"] == "RAI"].sort_values("r-value").head(amount_top_correlations)
top_corrs_sorted_descending = corrs[corrs["feature2"] == "RAI"].sort_values("r-value", ascending=False).head(amount_top_correlations)


top_corrs_zip = zip(top_corrs_sorted.analysis.tolist() + top_corrs_sorted_descending.analysis.tolist(),
                    top_corrs_sorted.feature1.tolist() + top_corrs_sorted_descending.feature1.tolist(),
                    top_corrs_sorted.feature2.tolist() + top_corrs_sorted_descending.feature2.tolist()
                   )


# # Take all of the features from the feature2 column of the above data frame and show how the 
# # r-value and p-value changes as N increases when correlating them with petal_width
# # for feature in dict_results["summary"].feature2.tolist():
# for cor_type, feat_1, feat_2 in top_corrs_zip:
#     summary, fig = correlations_as_sample_increases(data = trans_df, 
#                                                     feature1=feat_1, 
#                                                     feature2=feat_2, 
#                                                     method=cor_type.split(" ")[0].lower(),
#                                                     bootstrap = False,
#                                                     bootstrap_per_N = 20,
#                                                     starting_N=20, 
#                                                     random_state=42069, 
#                                                     plot=True)
    
# #     fig.savefig(f"Images/{feat_1}_{feat_2}_correlation_as_sample_increases.png", dpi=400, format="png", bbox_inches = "tight")


hue_column = "gender"
hue_column = None
# POM
_ = plt.figure(figsize=(10,4))
_ = sns.boxplot(x = "TIME", y = top_feature1, hue=hue_column, data=trans_df)
_ = sns.stripplot(x = "TIME", y = top_feature1, hue=hue_column, data=trans_df, dodge=True, linewidth=1, 
                  edgecolor="white")
# _ = plt.legend(loc="lower center", ncol=6)
_ = sns.despine()

# Fitbit feature
_ = plt.figure(figsize=(10,4))
_ = sns.boxplot(x = "TIME", y = top_feature2, hue=hue_column, data=trans_df)
_ = sns.stripplot(x = "TIME", y = top_feature2, hue=hue_column, data=trans_df, dodge=True, linewidth=1, edgecolor="white")
# _ = plt.legend(loc="lower center", ncol=6)
_ = sns.despine()


fig, ax, means, diffs = bland_altman_plot_repeated(data=trans_df[[top_feature1, top_feature2, "TIME"]],
                                                   column1=top_feature1,
                                                   column2=top_feature2,
                                                   group_column="TIME", 
#                                                    hue_column="diagnosis_ms",
                                                   limits_approach="anova",
                                                   aggregate_scatter=False,
                                                   aggregate_scatter_method="median",
                                                   plt_title=f"Bland-Altman mean difference plot between {top_feature1} and\n {top_feature2}, taking into account multiple time points",
                                                   figsize=(9, 5))

/Users/jamestwose/Coding/data-science-core-develop/neuropy/frequentist_statistics.py:2576: UserWarning: The amount of data in the columns is not equal across columns, hence the amount of data shown in the plot is equal to 210, column == Dedication
  warnings.warn(


target_name = "RAI"


feature_collection = " + ".join(features_list[3:-1])
feature_collection

'ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + RelatedFrust + Vigor + Dedication + Absorption + WorkEngagement + SelfControl + Congruence + Interest + Control + AutonomousFunctioning + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism'


trans_df[features_list[1:-1]].apply(lambda x: x.value_counts())
trans_df[features_list[1:-1]].describe()


scale_bool = True
if scale_bool:
    trans_df[features_list[3:-1] + [target_name]] = pd.DataFrame(MinMaxScaler().fit_transform(trans_df[features_list[3:-1] + [target_name]]), index=trans_df.index, columns=features_list[3:-1] + [target_name])
#     trans_df[features_list[3:-1] + [target_name]] = pd.DataFrame(StandardScaler().fit_transform(trans_df[features_list[3:-1] + [target_name]]), index=trans_df.index, columns=features_list[3:-1] + [target_name])

/opt/miniconda3/envs/general/lib/python3.9/site-packages/pandas/core/frame.py:3191: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self[k1] = value[k2]


md = smf.mixedlm(f"{target_name} ~ {feature_collection}", data=trans_df, groups=trans_df["User"], re_formula="~TIME", missing="drop")
# md = smf.mixedlm(f"RAI ~ Allocation", data=trans_df, groups=trans_df["User"], missing="drop")
mdf = md.fit(method=["lbfgs"])
print(mdf.summary())

                    Mixed Linear Model Regression Results
==============================================================================
Model:                    MixedLM         Dependent Variable:         RAI     
No. Observations:         180             Method:                     REML    
No. Groups:               94              Scale:                      0.0027  
Min. group size:          1               Log-Likelihood:             204.1497
Max. group size:          3               Converged:                  No      
Mean group size:          1.9                                                 
------------------------------------------------------------------------------
                       Coef.   Std.Err.    z    P>|z|    [0.025       0.975]  
------------------------------------------------------------------------------
Intercept              12.758 178230.804  0.000 1.000  -349313.198  349338.714
ExternalSocial         -0.108      0.037 -2.907 0.004       -0.181      -0.035
ExternalMaterial       -0.114      0.034 -3.328 0.001       -0.181      -0.047
AutonomySat             0.231      0.054  4.241 0.000        0.124       0.337
CompetenceSat           0.104      0.051  2.048 0.041        0.004       0.204
RelatedSat              0.077      0.047  1.649 0.099       -0.015       0.170
AutonomyFrust          -0.143      0.043 -3.324 0.001       -0.227      -0.059
CompetenceFrust         0.124      0.042  2.924 0.003        0.041       0.206
RelatedFrust            0.007      0.034  0.191 0.849       -0.061       0.074
Vigor                  -3.905 100669.667 -0.000 1.000  -197312.827  197305.016
Dedication             -3.934 108413.487 -0.000 1.000  -212490.465  212482.596
Absorption             -5.103 139388.770 -0.000 1.000  -273202.071  273191.865
WorkEngagement         11.382 302009.001  0.000 1.000  -591915.382  591938.146
SelfControl             0.023      0.042  0.540 0.589       -0.060       0.105
Congruence             20.200 313374.475  0.000 1.000  -614182.485  614222.885
Interest               26.895 417832.634  0.000 1.000  -818910.018  818963.809
Control               -24.294 376049.370 -0.000 1.000  -737067.516  737018.928
AutonomousFunctioning -52.517 814773.636 -0.000 1.000 -1596979.498 1596874.464
IncStructResources      0.112      0.042  2.689 0.007        0.030       0.193
DecHinderDemands       -0.101      0.032 -3.187 0.001       -0.163      -0.039
IncChallengeDemands     0.031      0.045  0.686 0.493       -0.058       0.120
IncSocialResources      0.054      0.035  1.534 0.125       -0.015       0.122
Resilience              0.067      0.039  1.703 0.088       -0.010       0.144
Optimism               -0.044      0.037 -1.204 0.229       -0.116       0.028
Group Var               0.003      0.004                                      
Group x TIME Cov        0.000      0.001                                      
TIME Var                0.000      0.001                                      
==============================================================================

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/base/model.py:566: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2206: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2218: ConvergenceWarning: Gradient optimization failed, |grad| = 294.669673
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2261: ConvergenceWarning: The Hessian matrix at the estimated parameter values is not positive definite.
  warnings.warn(msg, ConvergenceWarning)


fig, axs = diagnostic_plots(model_fit=mdf,
                                             X=None,
                                             y=None,
                                             figsize = (8,8),
                                             limit_cooks_plot = False,
                                             subplot_adjust_args={"wspace": 0.3, "hspace": 0.3}
                                            )


cutoff = 0.06
drop_amount = 1
fixed_effect_output = mdf.summary().tables[1]
fixed_effect_output_no_empty = fixed_effect_output[fixed_effect_output["P>|z|"] != ""].astype(float)


while fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].shape[0] > 0:
    feature_to_drop = fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].sort_values(by = "P>|z|").tail(drop_amount).index.tolist()
    sig_predictors = fixed_effect_output_no_empty.drop(["Intercept", feature_to_drop[0]]).index.tolist()

    feature_collection = " + ".join(sig_predictors)

    print(feature_collection)
    
    md = smf.mixedlm(f"{target_name} ~ {feature_collection}", data=trans_df, groups=trans_df["User"], re_formula="~TIME", missing="drop")
    mdf = md.fit(method=["lbfgs"])
    
    fixed_effect_output = mdf.summary().tables[1]
    fixed_effect_output_no_empty = fixed_effect_output[fixed_effect_output["P>|z|"] != ""].astype(float)
#     print(mdf.summary())

<ipython-input-51-b9bb0bc68d69>:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  while fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].shape[0] > 0:
<ipython-input-51-b9bb0bc68d69>:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  feature_to_drop = fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].sort_values(by = "P>|z|").tail(drop_amount).index.tolist()

ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + RelatedFrust + Vigor + Dedication + Absorption + WorkEngagement + SelfControl + Congruence + Interest + Control + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/base/model.py:566: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2206: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2218: ConvergenceWarning: Gradient optimization failed, |grad| = 29.675986
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2261: ConvergenceWarning: The Hessian matrix at the estimated parameter values is not positive definite.
  warnings.warn(msg, ConvergenceWarning)
<ipython-input-51-b9bb0bc68d69>:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  while fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].shape[0] > 0:
<ipython-input-51-b9bb0bc68d69>:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  feature_to_drop = fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].sort_values(by = "P>|z|").tail(drop_amount).index.tolist()

ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + RelatedFrust + Vigor + Dedication + Absorption + SelfControl + Congruence + Interest + Control + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:1634: UserWarning: Random effects covariance is singular
  warnings.warn(msg)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)
<ipython-input-51-b9bb0bc68d69>:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  while fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].shape[0] > 0:
<ipython-input-51-b9bb0bc68d69>:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  feature_to_drop = fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].sort_values(by = "P>|z|").tail(drop_amount).index.tolist()

ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + RelatedFrust + Vigor + Dedication + Absorption + SelfControl + Interest + Control + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)
<ipython-input-51-b9bb0bc68d69>:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  while fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].shape[0] > 0:
<ipython-input-51-b9bb0bc68d69>:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  feature_to_drop = fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].sort_values(by = "P>|z|").tail(drop_amount).index.tolist()

ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + Vigor + Dedication + Absorption + SelfControl + Interest + Control + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)
<ipython-input-51-b9bb0bc68d69>:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  while fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].shape[0] > 0:
<ipython-input-51-b9bb0bc68d69>:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  feature_to_drop = fixed_effect_output_no_empty.drop("Intercept")[fixed_effect_output_no_empty["P>|z|"] > cutoff].sort_values(by = "P>|z|").tail(drop_amount).index.tolist()

ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + Vigor + Dedication + Absorption + Interest + Control + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism


print(mdf.summary())

            Mixed Linear Model Regression Results
=============================================================
Model:               MixedLM   Dependent Variable:   RAI     
No. Observations:    180       Method:               REML    
No. Groups:          94        Scale:                0.0029  
Min. group size:     1         Log-Likelihood:       192.6853
Max. group size:     3         Converged:            Yes     
Mean group size:     1.9                                     
-------------------------------------------------------------
                   Coef.  Std.Err.   z    P>|z| [0.025 0.975]
-------------------------------------------------------------
Intercept           0.269    0.044  6.180 0.000  0.184  0.355
ExternalSocial     -0.103    0.033 -3.066 0.002 -0.168 -0.037
ExternalMaterial   -0.115    0.032 -3.608 0.000 -0.177 -0.053
AutonomySat         0.246    0.050  4.948 0.000  0.148  0.343
CompetenceSat       0.124    0.030  4.143 0.000  0.066  0.183
AutonomyFrust      -0.127    0.036 -3.498 0.000 -0.198 -0.056
CompetenceFrust     0.116    0.034  3.412 0.001  0.049  0.182
Dedication          0.194    0.029  6.767 0.000  0.138  0.250
IncStructResources  0.129    0.015  8.408 0.000  0.099  0.159
DecHinderDemands   -0.090    0.015 -6.204 0.000 -0.119 -0.062
Group Var           0.002                                    
Group x TIME Cov    0.000                                    
TIME Var            0.000                                    
=============================================================


# top_feature = fixed_effect_output_no_empty.sort_values("P>|z|").drop("Intercept").head(1).index.tolist()[0]
# top_feature


top_feature = "Dedication" #based on top correlation


md = smf.mixedlm(f"{target_name} ~ {top_feature}", data=trans_df, groups=trans_df["User"], re_formula="~TIME", missing="drop")
mdf = md.fit(method=["lbfgs"])

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)


print(mdf.summary())

          Mixed Linear Model Regression Results
==========================================================
Model:              MixedLM  Dependent Variable:  RAI     
No. Observations:   210      Method:              REML    
No. Groups:         105      Scale:               0.0051  
Min. group size:    1        Log-Likelihood:      172.0057
Max. group size:    3        Converged:           Yes     
Mean group size:    2.0                                   
----------------------------------------------------------
                 Coef. Std.Err.   z    P>|z| [0.025 0.975]
----------------------------------------------------------
Intercept        0.285    0.033  8.773 0.000  0.222  0.349
Dedication       0.432    0.040 10.786 0.000  0.354  0.511
Group Var        0.007    0.123                           
Group x TIME Cov 0.001    0.042                           
TIME Var         0.000    0.023                           
==========================================================


fig, axs = diagnostic_plots(model_fit=mdf,
                                             X=None,
                                             y=None,
                                             figsize = (8,8),
                                             limit_cooks_plot = False,
                                             subplot_adjust_args={"wspace": 0.3, "hspace": 0.3}
                                            )
fig.savefig(f"images/{target_name}_{top_feature}_model_diagnostics.png", dpi=400, format="png")


df_test = pd.DataFrame({"User": trans_df.loc[mdf.fittedvalues.index, "User"],
    "fitted": mdf.fittedvalues,
              target_name: trans_df.loc[mdf.fittedvalues.index, target_name]}).set_index("User")


plt.figure(figsize=(30,6))
plt.title(f"Linear Mixed Model | Outcome Measure == {target_name} | RMSE = {round(np.sqrt(mean_squared_error(df_test['fitted'], df_test[target_name])),2)} | bias Error = { round(np.mean(df_test['fitted'] - df_test[target_name]), 2)} ")
markerline, stemlines, baseline = plt.stem(df_test.index, df_test['fitted'] - df_test[target_name], use_line_collection=True, linefmt="black", markerfmt='D')
plt.setp(markerline, "alpha", 0.5)

plt.hlines(y=round(np.sqrt(mean_squared_error(df_test['fitted'], df_test[target_name])),2), colors="grey", linestyles='-.', label='+ RMSE', 
           xmin = df_test.index.min(), xmax = "Moti201"
          ) 
plt.hlines(y=round(-np.sqrt(mean_squared_error(df_test['fitted'], df_test[target_name])),2), colors="grey", linestyles='-.', label='- RMSE', 
            xmin = df_test.index.min(), xmax = "Moti201"
          ) 
plt.xticks(rotation=90, ticks=df_test.index)
plt.ylabel(f"'Error = fitted values - real {target_name}'")
# plt.ylim([-.5,.6])
_ = plt.legend()
plt.savefig(f"images/{target_name}_fitted_RMSE_plot.png", dpi=400, format="png")


trans_df.loc[mdf.fittedvalues.index, target_name].min(), trans_df.loc[mdf.fittedvalues.index, target_name].max()

(0.0, 1.0)


fit_errors_df = pd.DataFrame({"RMSE": np.sqrt(mean_squared_error(df_test['fitted'], df_test[target_name])), 
              "MAE": mean_absolute_error(df_test['fitted'], df_test[target_name]),
              "bias error": np.mean(df_test['fitted'] - df_test[target_name])},
              index=[target_name])
fit_errors_df.round(2)


_ = sns.lmplot(x=top_feature, 
               y=target_name, 
               hue="TIME", 
               data=trans_df,
               ci=None,
               scatter_kws={'s':100}, 
               height=6, 
               aspect=1.5,
              legend=True)
_ = sns.regplot(x=top_feature, 
                y=target_name, 
                data=trans_df,
                color="black", 
                scatter_kws={'s': 20}, 
                line_kws={"linestyle": "--"})
_ = plt.title(f"{target_name} vs {top_feature.replace('_', ' ').title()}\nblack line is the regression line of the group | coloured lines are regression lines within TIME")


_ = sns.lmplot(x=top_feature, 
               y=target_name, 
               hue="User", 
               data=trans_df,
               ci=None,
               scatter_kws={'s':100}, 
               height=6, 
               aspect=1.5,
              legend=False)
_ = sns.regplot(x=top_feature, 
                y=target_name, 
                data=trans_df,
                color="black", 
                scatter_kws={'s': 20}, 
                line_kws={"linestyle": "--"})
_ = plt.title(f"{target_name} vs {top_feature.replace('_', ' ').title()}\nblack line is the regression line of the group | coloured lines are regression lines within User")


top_feature = "Dedication" #based on top correlation

x_train_index = trans_df.sample(frac=0.8, random_state=42).index

x_test_index = trans_df.drop(x_train_index, axis=0).index


md = smf.mixedlm(f"{target_name} ~ {top_feature}", data=trans_df.loc[x_train_index, :], groups=trans_df.loc[x_train_index, "User"], re_formula="~TIME", missing="drop")
mdf = md.fit(method=["lbfgs"])

/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/base/model.py:566: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2206: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2218: ConvergenceWarning: Gradient optimization failed, |grad| = 0.933821
  warnings.warn(msg, ConvergenceWarning)
/opt/miniconda3/envs/general/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.
  warnings.warn(msg, ConvergenceWarning)


print(mdf.summary())

          Mixed Linear Model Regression Results
=========================================================
Model:             MixedLM  Dependent Variable:  RAI     
No. Observations:  167      Method:              REML    
No. Groups:        97       Scale:               0.0043  
Min. group size:   1        Log-Likelihood:      140.8881
Max. group size:   3        Converged:           No      
Mean group size:   1.7                                   
---------------------------------------------------------
                 Coef. Std.Err.   z   P>|z| [0.025 0.975]
---------------------------------------------------------
Intercept        0.318    0.036 8.835 0.000  0.248  0.389
Dedication       0.397    0.044 9.099 0.000  0.312  0.483
Group Var        0.007    0.095                          
Group x TIME Cov 0.001    0.027                          
TIME Var         0.000    0.014                          
=========================================================


df_test = pd.DataFrame({"User": trans_df.loc[x_test_index, "User"],
    "predicted": mdf.predict(trans_df.loc[x_test_index, "Dedication"]),
              target_name: trans_df.loc[x_test_index, target_name]}).set_index("User").dropna()


plt.figure(figsize=(30,6))
plt.title(f"Linear Mixed Model | Outcome Measure == {target_name} | RMSE = {round(np.sqrt(mean_squared_error(df_test['predicted'], df_test[target_name])),2)} | bias Error = { round(np.mean(df_test['predicted'] - df_test[target_name]), 2)} ")
markerline, stemlines, baseline = plt.stem(df_test.index, df_test['predicted'] - df_test[target_name], use_line_collection=True, linefmt="black", markerfmt='D')
plt.setp(markerline, "alpha", 0.5)

plt.hlines(y=round(np.sqrt(mean_squared_error(df_test['predicted'], df_test[target_name])),2), colors="grey", linestyles='-.', label='+ RMSE', 
           xmin = df_test.index.min(), xmax = "Moti201"
          ) 
plt.hlines(y=round(-np.sqrt(mean_squared_error(df_test['predicted'], df_test[target_name])),2), colors="grey", linestyles='-.', label='- RMSE', 
            xmin = df_test.index.min(), xmax = "Moti201"
          ) 
plt.xticks(rotation=90, ticks=df_test.index)
plt.ylabel(f"'Error = predicted values - real {target_name}'")
# plt.ylim([-.5,.6])
_ = plt.legend()
plt.savefig(f"images/{target_name}_pred_test_RMSE_plot.png", dpi=400, format="png")


fit_errors_df = pd.DataFrame({"RMSE": np.sqrt(mean_squared_error(df_test['predicted'], df_test[target_name])), 
              "MAE": mean_absolute_error(df_test['predicted'], df_test[target_name]),
              "bias error": np.mean(df_test['predicted'] - df_test[target_name])},
              index=[target_name])
fit_errors_df.round(2)


trans_df.loc[x_train_index, :].User.sort_values().unique()

array(['AAMoti?1', 'AAMoti?4', 'AAMoti?5', 'AAMoti?6', 'Moti105',
       'Moti106', 'Moti108', 'Moti109', 'Moti110', 'Moti111', 'Moti112',
       'Moti113', 'Moti114', 'Moti115', 'Moti116', 'Moti117', 'Moti118',
       'Moti119', 'Moti120', 'Moti121', 'Moti122', 'Moti123', 'Moti125',
       'Moti126', 'Moti127', 'Moti128', 'Moti129', 'Moti130', 'Moti131',
       'Moti132', 'Moti134', 'Moti135', 'Moti137', 'Moti138', 'Moti140',
       'Moti141', 'Moti143', 'Moti144', 'Moti146', 'Moti147', 'Moti148',
       'Moti149', 'Moti150', 'Moti151', 'Moti152', 'Moti153', 'Moti154',
       'Moti155', 'Moti156', 'Moti157', 'Moti159', 'Moti160', 'Moti161',
       'Moti162', 'Moti163', 'Moti164', 'Moti165', 'Moti166', 'Moti167',
       'Moti168', 'Moti169', 'Moti170', 'Moti171', 'Moti172', 'Moti173',
       'Moti175', 'Moti176', 'Moti177', 'Moti178', 'Moti179', 'Moti180',
       'Moti181', 'Moti182', 'Moti183', 'Moti184', 'Moti185', 'Moti186',
       'Moti187', 'Moti188', 'Moti189', 'Moti190', 'Moti191', 'Moti192',
       'Moti193', 'Moti194', 'Moti195', 'Moti197', 'Moti198', 'Moti199',
       'Moti200', 'Moti202', 'Moti203', 'Moti205', 'Moti206', 'Moti207',
       'Moti208', 'Moti209', 'Moti212', 'Moti213', 'Moti215', 'Moti216',
       'Moti217', 'Moti218'], dtype=object)


trans_df.loc[x_test_index, :].User.sort_values().unique()

array(['AAMoti?2', 'Moti109', 'Moti110', 'Moti114', 'Moti115', 'Moti116',
       'Moti117', 'Moti118', 'Moti121', 'Moti122', 'Moti123', 'Moti127',
       'Moti143', 'Moti144', 'Moti145', 'Moti147', 'Moti150', 'Moti156',
       'Moti158', 'Moti160', 'Moti163', 'Moti164', 'Moti165', 'Moti167',
       'Moti168', 'Moti173', 'Moti178', 'Moti183', 'Moti187', 'Moti188',
       'Moti190', 'Moti191', 'Moti196', 'Moti201', 'Moti203', 'Moti204',
       'Moti207', 'Moti210', 'Moti211', 'Moti218'], dtype=object)


# set(trans_df.loc[x_train_index, :].User.unique().tolist()) - set(trans_df.loc[x_test_index, :].User.unique().tolist())
set(trans_df.loc[x_test_index, :].User.unique().tolist()) - set(trans_df.loc[x_train_index, :].User.unique().tolist()).intersection(set(trans_df.loc[x_test_index, :].User.unique().tolist()))

{'AAMoti?2',
 'Moti145',
 'Moti158',
 'Moti196',
 'Moti201',
 'Moti204',
 'Moti210',
 'Moti211'}


0.14 - 0.06

0.08000000000000002

	Allocation	age	ExternalSocial	ExternalMaterial	AutonomySat	CompetenceSat	RelatedSat	AutonomyFrust	CompetenceFrust	RelatedFrust	...	Congruence	Interest	Control	AutonomousFunctioning	IncStructResources	DecHinderDemands	IncChallengeDemands	IncSocialResources	Resilience	Optimism
count	226.000000	95.000000	212.000000	212.000000	209.000000	210.000000	209.000000	209.000000	209.000000	209.000000	...	207.000000	207.000000	207.000000	207.000000	181.000000	181.000000	181.000000	181.000000	181.000000	181.000000
mean	0.650442	46.115789	3.701258	3.166667	5.129984	5.654762	5.455742	2.801435	2.330542	2.262360	...	3.850483	3.421256	2.201208	5.070531	5.692818	3.046961	4.907182	4.020718	4.550645	5.298343
std	0.477889	9.934748	1.410806	1.336686	1.016055	0.922137	1.006954	1.273122	1.126968	0.965801	...	0.489228	0.897546	0.805208	1.245400	0.828069	1.211355	1.217741	1.286867	0.921309	1.374145
min	0.000000	23.000000	1.000000	1.000000	1.250000	2.750000	2.000000	1.000000	1.000000	1.000000	...	2.000000	1.000000	1.000000	0.400000	2.800000	1.000000	1.000000	1.400000	1.666667	1.000000
25%	0.000000	40.000000	2.666667	2.000000	4.500000	5.250000	5.000000	1.750000	1.500000	1.500000	...	3.600000	2.800000	1.600000	4.400000	5.200000	2.166667	4.200000	3.000000	4.000000	4.500000
50%	1.000000	47.000000	3.666667	3.333333	5.250000	6.000000	5.750000	2.750000	2.000000	2.000000	...	4.000000	3.600000	2.200000	5.000000	5.800000	3.000000	5.000000	4.000000	4.666667	5.500000
75%	1.000000	55.000000	4.666667	4.000000	6.000000	6.250000	6.000000	3.750000	3.000000	2.750000	...	4.000000	4.000000	2.600000	5.900000	6.400000	3.833333	5.800000	5.000000	5.000000	6.000000
max	1.000000	63.000000	7.000000	6.333333	7.000000	7.000000	7.000000	7.000000	5.500000	5.000000	...	5.000000	5.000000	4.600000	8.200000	7.000000	6.166667	7.000000	7.000000	7.000000	7.000000

Work Motivation Under Lockdown Exploratory Data Analysis (EDA)¶

This EDA only addresses the daily data, specifically `data=="data/moti_feasibility_james_daily.csv"`¶

Read in the raw daily data¶

Reshape the data frame and convert data types to numeric where possible¶

Count the amount of rows per user and plot¶

Plot the heatmaps of the top 5 users with the highest amount of rows¶

Visualize outliers¶

Check the normality using the Kolmogrov-Smirnov test¶

Run the correlations as sample size increases for the highest correlations¶

Level of aggreement between two clinical measurements accounting for multiple data points¶

Specify the model with the formula_string as the model syntax¶

Create a train and test set, use a fitted model to predict previously unseen rows (this is not controlled for user)¶

Users included in the training set¶

Users included in the test set¶

Users in the test set that were not included in the training set¶

	User	TIME	date	Allocation	gender	age	education	ExternalSocial	ExternalMaterial	AutonomySat	...	Interest	Control	AutonomousFunctioning	IncStructResources	DecHinderDemands	IncChallengeDemands	IncSocialResources	Resilience	Optimism	RAI
0	AAMoti?1	1.0	2020-03-27 09:01:15	NaN	0.0	54.0	A4	3.0	1.333333	4.50	...	2.8	1.6	4.8	6.2	3.0	4.6	4.2	4.666667	5.5	10.75
1	AAMoti?2	1.0	2020-04-23 13:16:36	NaN	0.0	41.0	A5	5.0	2.000000	3.75	...	4.6	3.6	5.0	4.0	2.0	2.8	2.2	4.666667	1.0	1.25

	var	method	outliers-removed	outliers-removed %	index
17	Congruence	inter-quartile range, k=3	3	1.4	[18, 105, 175]
20	AutonomousFunctioning	standard deviation, k=3.0	2	1.0	[105, 175]
9	AutonomyFrust	standard deviation, k=3.0	2	0.9	[17, 226]
14	Absorption	inter-quartile range, k=3	1	0.5	[167]
15	WorkEngagement	inter-quartile range, k=3	1	0.5	[80]

	0	1	2	3	4	5	6	7	8	9	...	18	19	20	21	22	23	24	25	26	27
feature	TIME	Allocation	gender	age	ExternalSocial	ExternalMaterial	AutonomySat	CompetenceSat	RelatedSat	AutonomyFrust	...	Interest	Control	AutonomousFunctioning	IncStructResources	DecHinderDemands	IncChallengeDemands	IncSocialResources	Resilience	Optimism	RAI
p-value	0.0	0.0	0.0	0.275094	0.103033	0.10707	0.032465	0.000052	0.000212	0.07133	...	0.046558	0.082886	0.276185	0.029861	0.435267	0.047239	0.218575	0.017274	0.000186	0.609656
normality	False	False	False	True	True	True	False	False	False	True	...	False	True	True	False	True	False	True	False	False	True

	analysis	feature1	feature2	r-value	p-value	stat-sign	N
286	Spearman Rank	Dedication	RAI	0.693083	2.148071e-31	True	210
167	Spearman Rank	AutonomySat	RAI	0.676722	2.418086e-29	True	209
187	Spearman Rank	CompetenceSat	RAI	0.628806	1.635579e-24	True	210
311	Spearman Rank	WorkEngagement	RAI	0.625542	3.316445e-24	True	210
224	Pearson	AutonomyFrust	RAI	-0.619377	1.573851e-23	True	209
272	Spearman Rank	Vigor	RAI	0.571688	1.277075e-19	True	210
241	Spearman Rank	CompetenceFrust	RAI	-0.514979	1.494000e-15	True	209
349	Pearson	Control	RAI	-0.496807	2.663397e-14	True	207
362	Spearman Rank	IncStructResources	RAI	0.526186	2.802045e-14	True	181
377	Spearman Rank	Optimism	RAI	0.518632	7.506659e-14	True	181
371	Pearson	IncChallengeDemands	RAI	0.497916	9.916126e-13	True	181
322	Spearman Rank	SelfControl	RAI	0.452299	7.829674e-12	True	207
299	Spearman Rank	Absorption	RAI	0.411519	5.469874e-10	True	210
124	Pearson	ExternalSocial	RAI	-0.378030	1.315857e-08	True	212
332	Spearman Rank	Congruence	RAI	0.374809	2.635422e-08	True	207
367	Pearson	DecHinderDemands	RAI	-0.393722	4.172293e-08	True	181
206	Spearman Rank	RelatedSat	RAI	0.354944	1.339931e-07	True	209
356	Pearson	AutonomousFunctioning	RAI	0.313403	4.262429e-06	True	207
146	Pearson	ExternalMaterial	RAI	-0.296031	1.168136e-05	True	212
257	Spearman Rank	RelatedFrust	RAI	-0.297875	1.184926e-05	True	209

Work Motivation Under Lockdown Exploratory Data Analysis (EDA)¶

This EDA only addresses the daily data, specifically data=="data/moti_feasibility_james_daily.csv"¶

Read in the raw daily data¶

Reshape the data frame and convert data types to numeric where possible¶

Count the amount of rows per user and plot¶

Plot the heatmaps of the top 5 users with the highest amount of rows¶

Visualize outliers¶

Check the normality using the Kolmogrov-Smirnov test¶

Run the correlations as sample size increases for the highest correlations¶

Level of aggreement between two clinical measurements accounting for multiple data points¶

Specify the model with the formula_string as the model syntax¶

Create a train and test set, use a fitted model to predict previously unseen rows (this is not controlled for user)¶

Users included in the training set¶

Users included in the test set¶

Users in the test set that were not included in the training set¶

This EDA only addresses the daily data, specifically `data=="data/moti_feasibility_james_daily.csv"`¶