James Twose
1.7 MiB
1.7 MiB
Work Motivation Under Lockdown Exploratory Data Analysis (EDA)¶
Two datasets...
Post-workday evaluations of work motivation related stuff. Large spread across how many time points people collected; few decent-ish ones.
Traditional questionnaires from the same participants.
- Allocation indicates whether they took part in our motivation self-management training during the period or belonged to a waitlist control group.
- Pre-, mid- and post-questionnaires --> TIME == 1, 2 and 3 respectively.
- Everything from intrinsic1 to Psycap5 are raw answers. From Intrinsic to Optimism are sumscores.
- This was a feasibility & acceptability study so it's not powered for or aiming to look at outcomes. But if there was one, it could be the "Relative Autonomy Index" (RAI), a sumscore conventionally built as:
- Amotivation * -3 + ExternalTotal * -2 + Introjected * -1 + Identified * 2 + Intrinsic * 3
- If I remember correctly, there were some people whose RAI was boosted a lot
- In general:
- Good things to go up during the intervention (or life in general):
- Intrinsic, Identified
- Competence, Relatedness, Autonomy
- Bad things to go up during the intervention (or life in general):
- Amotivation, External (material & social), Introjected
- AutonomyThwarting, RelatednessThwarting, CompetenceThwarting
- Unsure if these were named in the data, or just some of the items named autonomyx/relatednessx/competencex
- Good things to go up during the intervention (or life in general):
This EDA only addresses the daily data, specifically data=="data/moti_feasibility_james_daily.csv"¶
In [1]:
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)
In [2]:
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
In [3]:
from sklearn.metrics import mean_squared_error, mean_absolute_error
Show the session information of the packages used in this analysis
In [4]:
session_info.show(write_req_file=False, req_file_name="work_motivation_pre_mid_post_EDA_requirements.txt",)
Out[4]:
Click to view session information
----- 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 -----
Click to view modules imported as dependencies
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
Read in the raw daily data¶
In [5]:
_df = pd.read_csv("data/moti_feasibility_james_pre-mid-post.csv")
In [6]:
_df = _df.assign(date=pd.to_datetime(_df["datestamp"], format="%Y-%m-%d %H:%M:%S"))
In [7]:
def convert_numeric_where_possible(x): try: return x.astype(np.number) except: return x
Reshape the data frame and convert data types to numeric where possible¶
In [8]:
demographics_columns = ['User', 'TIME', 'date', 'Allocation', 'gender', 'age', 'education']
In [9]:
target_columns = ["Amotivation", "ExternalTotal", "Introjected", "Identified", "Intrinsic"]
In [10]:
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)
Out[10]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| 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 |
2 rows × 31 columns
In [11]:
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
In [12]:
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 % =================
Count the amount of rows per user and plot¶
The aim of this is to assess whether there is enough data to do some of the more complext NLTSA methods
In [13]:
row_amount_df = (df .reset_index() .groupby("User") .count() .loc[:, ["date"]] .rename(columns={"date": "row_amount"}) .sort_values(by="row_amount") .reset_index())
In [14]:
row_amount_df.tail(1)
Out[14]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| User | row_amount | |
|---|---|---|
| 111 | Moti218 | 3 |
In [15]:
_ = 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()
Plot the heatmaps of the top 5 users with the highest amount of rows¶
The aim of this is to assess the amount of missingness time wise (i.e. the user could have a lot of data, spread over a long period with large gaps in the middle).
In [16]:
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")
In [17]:
plot_df = (df .set_index("User") .select_dtypes(np.number) .reset_index() .melt(id_vars=["User", "TIME"]))
In [18]:
_ = plt.figure(figsize=(20, 10)) _ = sns.lineplot(data=plot_df, x="TIME", y="value", hue="User", legend=False)
In [19]:
# _ = plt.figure(figsize=(20, 10)) # _ = sns.barplot(data=plot_df, x="variable", y="value") # _ = plt.xticks(rotation=90)
In [20]:
_ = plt.figure(figsize=(20, 10)) _ = sns.boxplot(data=plot_df, x="variable", y="value") _ = plt.xticks(rotation=90)
Visualize outliers¶
In [21]:
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')
In [22]:
O_summary.sort_values(by = "outliers-removed %", ascending=False).head(5)
Out[22]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| 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] |
Check the normality using the Kolmogrov-Smirnov test¶
In [23]:
# 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]}")
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| 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 |
3 rows × 28 columns
====== Normality test based on Kolmogorov-Smirnov test ===== - Size normally distributed feature: 9 - Size NOT normally distributed feature: 19
In [24]:
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']
In [25]:
len(X_normal.columns), len(use_transformation), len(use_nonparametric), (len(X_normal.columns) + len(use_transformation) + len(use_nonparametric))
Out[25]:
(9, 4, 15, 28)
In [26]:
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) )
In [27]:
trans_df["gender"].value_counts()
Out[27]:
0.0 75 1.0 23 3.0 3 Name: gender, dtype: int64
In [28]:
trans_df = df[df["gender"] != 3] # df["gender"] = df.gender.astype(int)
In [29]:
trans_df["gender"].value_counts()
Out[29]:
0.0 75 1.0 23 Name: gender, dtype: int64
In [30]:
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')
In [31]:
transformed_fitbit_cols = X_normal.columns.tolist() + use_transformation + use_nonparametric np.array(transformed_fitbit_cols)
Out[31]:
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')
In [32]:
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())
In [33]:
corrs[corrs["feature2"] == "RAI"].sort_values("p-value").head(20)
Out[33]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| 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 |
In [34]:
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())
In [35]:
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]
In [36]:
_ = 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()
In [37]:
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]
In [38]:
_ = 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()
Run the correlations as sample size increases for the highest correlations¶
In [39]:
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)
In [40]:
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() )
In [41]:
# # 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")
In [42]:
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()
Level of aggreement between two clinical measurements accounting for multiple data points¶
In [43]:
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(
In [44]:
target_name = "RAI"
In [45]:
feature_collection = " + ".join(features_list[3:-1]) feature_collection
Out[45]:
'ExternalSocial + ExternalMaterial + AutonomySat + CompetenceSat + RelatedSat + AutonomyFrust + CompetenceFrust + RelatedFrust + Vigor + Dedication + Absorption + WorkEngagement + SelfControl + Congruence + Interest + Control + AutonomousFunctioning + IncStructResources + DecHinderDemands + IncChallengeDemands + IncSocialResources + Resilience + Optimism'
In [46]:
trans_df[features_list[1:-1]].apply(lambda x: x.value_counts()) trans_df[features_list[1:-1]].describe()
Out[46]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| 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 |
8 rows × 25 columns
In [47]:
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]
Specify the model with the formula_string as the model syntax¶
In [48]:
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)
In [49]:
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} )
In [50]:
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)
In [51]:
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
/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 + 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/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.002142
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)
<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 + Control + IncStructResources + DecHinderDemands + 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 + Control + IncStructResources + DecHinderDemands + IncSocialResources + Resilience
/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 + Control + IncStructResources + DecHinderDemands + Resilience
/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 + IncStructResources + DecHinderDemands + Resilience
/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 + Dedication + Absorption + IncStructResources + DecHinderDemands + Resilience
/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 + Dedication + Absorption + IncStructResources + DecHinderDemands
/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 + AutonomyFrust + CompetenceFrust + Dedication + Absorption + IncStructResources + DecHinderDemands
/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 + AutonomyFrust + CompetenceFrust + Dedication + IncStructResources + DecHinderDemands
/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:
In [52]:
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
=============================================================
In [53]:
# top_feature = fixed_effect_output_no_empty.sort_values("P>|z|").drop("Intercept").head(1).index.tolist()[0] # top_feature
In [54]:
top_feature = "Dedication" #based on top correlation
In [55]:
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)
In [56]:
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
==========================================================
In [57]:
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")
In [59]:
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")
In [60]:
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")
In [61]:
trans_df.loc[mdf.fittedvalues.index, target_name].min(), trans_df.loc[mdf.fittedvalues.index, target_name].max()
Out[61]:
(0.0, 1.0)
In [62]:
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)
Out[62]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| RMSE | MAE | bias error | |
|---|---|---|---|
| RAI | 0.06 | 0.04 | 0.0 |
In [63]:
_ = 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")
In [64]:
_ = 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")
Create a train and test set, use a fitted model to predict previously unseen rows (this is not controlled for user)¶
In [65]:
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
In [66]:
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)
In [67]:
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
=========================================================
In [68]:
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()
In [69]:
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")
In [70]:
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)
Out[70]:
<style scoped="">
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
| RMSE | MAE | bias error | |
|---|---|---|---|
| RAI | 0.14 | 0.1 | 0.01 |
Users included in the training set¶
In [84]:
trans_df.loc[x_train_index, :].User.sort_values().unique()
Out[84]:
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)
Users included in the test set¶
In [85]:
trans_df.loc[x_test_index, :].User.sort_values().unique()
Out[85]:
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)
Users in the test set that were not included in the training set¶
In [89]:
# 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()))
Out[89]:
{'AAMoti?2',
'Moti145',
'Moti158',
'Moti196',
'Moti201',
'Moti204',
'Moti210',
'Moti211'}
In [90]:
0.14 - 0.06
Out[90]:
0.08000000000000002