James
6.1 MiB
6.1 MiB
Windowed PCA and fluctuation intensity in COVID19 data¶
In [111]:
import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from jmspack.utils import apply_scaling, JmsColors from jmspack.NLTSA import (fluctuation_intensity, distribution_uniformity, complexity_resonance, complexity_resonance_diagram, cumulative_complexity_peaks, cumulative_complexity_peaks_plot) from sklearn import decomposition from sklearn.linear_model import LinearRegression from scipy.stats import norm import time
In [112]:
start = time.time()
In [113]:
if "jms_style_sheet" in plt.style.available: plt.style.use("jms_style_sheet")
In [114]:
# df = pd.read_csv("../shield-complexity/data/Finland_worldsurvey_nonmissing_since_2021-06-08_to_2022-01-20.csv") df = pd.read_csv("../shield-complexity/data/Finland_worldsurvey_nonmissing_c_of_v_since_2021-06-08_to_2022-01-20.csv")
In [115]:
df["date"] = pd.to_datetime(df["date"]).dt.date
In [116]:
display(df.head()); df.shape
<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>
| date | pct_covid | pct_flu | percent_mc | percent_hf | percent_anos | pct_twodoses | pct_community_cli | pct_activity_work_outside_home | pct_activity_shop | ... | pct_want_info_covid_variants | pct_want_info_children_education | pct_want_info_economic_impact | pct_want_info_mental_health | pct_want_info_relationships | pct_want_info_employment | pct_want_info_none | pct_delayed_care_cost | pct_mask_shop_1d | pct_mask_spent_time_1d | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2021-06-09 | 9.039486 | 9.039486 | 0.612679 | 2.813836 | 6.958875 | 2.033529 | 4.412780 | 1.701158 | 0.408198 | ... | 1.796540 | 3.263608 | 3.377915 | 2.276444 | 2.927897 | 3.459269 | 1.471821 | 2.490908 | 0.552925 | 1.417980 |
| 1 | 2021-06-10 | 11.782998 | 11.782998 | 0.512932 | 2.422617 | 8.330548 | 2.013005 | 5.110684 | 1.874262 | 0.391850 | ... | 1.617688 | 3.455750 | 4.811816 | 2.227195 | 2.701875 | 3.445995 | 1.478371 | 2.153090 | 0.386451 | 1.578405 |
| 2 | 2021-06-11 | 19.488523 | 15.235526 | 2.173259 | 2.121825 | 7.898911 | 3.045827 | 4.672163 | 2.991276 | 2.263277 | ... | 1.611505 | 3.279128 | 5.133363 | 2.492608 | 2.942794 | 4.426361 | 1.449034 | 2.251560 | 0.236805 | 1.628893 |
| 3 | 2021-06-12 | 13.118015 | 13.118015 | 1.162009 | 2.543860 | 8.630301 | 1.826629 | 5.658765 | 2.148999 | 0.995366 | ... | 1.701548 | 3.803658 | 3.688413 | 2.070748 | 2.775816 | 3.189387 | 1.242110 | 2.605254 | 0.580646 | 2.399832 |
| 4 | 2021-06-13 | 10.705081 | 12.118662 | 0.782932 | 2.329260 | 8.237337 | 2.684727 | 5.455478 | 2.414263 | 0.566504 | ... | 1.922228 | 3.769134 | 3.916612 | 2.752952 | 3.658005 | 3.275252 | 1.477509 | 2.832298 | 0.449128 | 2.521729 |
5 rows × 64 columns
Out[116]:
(206, 64)
In [117]:
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 206 entries, 0 to 205 Data columns (total 64 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 date 206 non-null object 1 pct_covid 206 non-null float64 2 pct_flu 206 non-null float64 3 percent_mc 206 non-null float64 4 percent_hf 206 non-null float64 5 percent_anos 206 non-null float64 6 pct_twodoses 206 non-null float64 7 pct_community_cli 206 non-null float64 8 pct_activity_work_outside_home 206 non-null float64 9 pct_activity_shop 206 non-null float64 10 pct_activity_restaurant_bar 206 non-null float64 11 pct_activity_spent_time 206 non-null float64 12 pct_activity_large_event 206 non-null float64 13 pct_activity_public_transit 206 non-null float64 14 pct_food_security 206 non-null float64 15 pct_anxious_7d 206 non-null float64 16 pct_depressed_7d 206 non-null float64 17 pct_symp_fever 206 non-null float64 18 pct_symp_cough 206 non-null float64 19 pct_symp_diff_breathing 206 non-null float64 20 pct_symp_fatigue 206 non-null float64 21 pct_symp_stuffy_nose 206 non-null float64 22 pct_symp_aches 206 non-null float64 23 pct_symp_sore_throat 206 non-null float64 24 pct_symp_chest_pain 206 non-null float64 25 pct_symp_nausea 206 non-null float64 26 pct_symp_headache 206 non-null float64 27 pct_symp_chills 206 non-null float64 28 pct_testing_rate 206 non-null float64 29 pct_avoid_contact 206 non-null float64 30 pct_worried_catch_covid 206 non-null float64 31 pct_belief_distancing_effective 206 non-null float64 32 pct_belief_masking_effective 206 non-null float64 33 pct_others_distanced_public 206 non-null float64 34 pct_others_masked_public 206 non-null float64 35 pct_belief_children_immune 206 non-null float64 36 pct_belief_no_spread_hot_humid 206 non-null float64 37 pct_received_news_local_health 206 non-null float64 38 pct_received_news_experts 206 non-null float64 39 pct_received_news_who 206 non-null float64 40 pct_received_news_govt_health 206 non-null float64 41 pct_received_news_politicians 206 non-null float64 42 pct_received_news_journalists 206 non-null float64 43 pct_received_news_friends 206 non-null float64 44 pct_received_news_none 206 non-null float64 45 pct_trust_covid_info_local_health 206 non-null float64 46 pct_trust_covid_info_experts 206 non-null float64 47 pct_trust_covid_info_who 206 non-null float64 48 pct_trust_covid_info_govt_health 206 non-null float64 49 pct_trust_covid_info_politicians 206 non-null float64 50 pct_trust_covid_info_journalists 206 non-null float64 51 pct_trust_covid_info_friends 206 non-null float64 52 pct_trust_covid_info_religious 206 non-null float64 53 pct_want_info_covid_treatment 206 non-null float64 54 pct_want_info_covid_variants 206 non-null float64 55 pct_want_info_children_education 206 non-null float64 56 pct_want_info_economic_impact 206 non-null float64 57 pct_want_info_mental_health 206 non-null float64 58 pct_want_info_relationships 206 non-null float64 59 pct_want_info_employment 206 non-null float64 60 pct_want_info_none 206 non-null float64 61 pct_delayed_care_cost 206 non-null float64 62 pct_mask_shop_1d 206 non-null float64 63 pct_mask_spent_time_1d 206 non-null float64 dtypes: float64(63), object(1) memory usage: 103.1+ KB
In [118]:
_ = plt.figure(figsize=(20, 10)) _ = sns.heatmap(df .set_index("date") .T )
In [119]:
_ = plt.figure(figsize=(20, 10)) _ = sns.heatmap(df .set_index("date") .pipe(apply_scaling) .T )
In [120]:
df_melt = (df .set_index("date") .pipe(apply_scaling) .reset_index() .melt(id_vars="date"))
In [121]:
_ = plt.figure(figsize=(30, 7)) _ = sns.lineplot(data=df_melt, x="date", y="value", hue="variable", legend=False)
In [122]:
g = sns.FacetGrid(df_melt, col="variable", col_wrap=3, aspect=2) _ = g.map(sns.lineplot, "date", "value")
In [123]:
def summary_window_FUN(x: pd.DataFrame, window_size: int = 7, user_func=decomposition.PCA, kwargs: dict = {}): window_range = np.arange(0, len(x)-window_size, window_size) cp_df = pd.DataFrame() for window_begin in window_range: current_cp_df = pd.DataFrame(user_func(n_components=None, **kwargs) .fit_transform(x.iloc[window_begin: window_begin + window_size])).iloc[:, 0] cp_df = pd.concat([cp_df, current_cp_df]) return cp_df.rename(columns={0: f"windowed_{user_func.__name__}"}).reset_index(drop=True)
In [124]:
def window_PCA(x: pd.DataFrame, window_size: int = 7, user_func=decomposition.PCA, kwargs: dict = {}): window_range = np.arange(0, len(x)-window_size) pca_list = list() for window_begin in window_range: current_pca = user_func(n_components=None, **kwargs).fit(x.iloc[window_begin: window_begin + window_size]).explained_variance_[0] pca_list.append(current_pca) pca_df = pd.DataFrame({f"windowed_{user_func.__name__}": pca_list}) return pca_df
In [125]:
window=28 pca_df = window_PCA(x=df.set_index("date"), window_size=window).set_index(df.loc[:df.shape[0]-window-1, "date"])
In [126]:
pca_df.head()
Out[126]:
<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>
| windowed_PCA | |
|---|---|
| date | |
| 2021-06-09 | 39.595692 |
| 2021-06-10 | 37.623526 |
| 2021-06-11 | 38.198733 |
| 2021-06-12 | 38.344276 |
| 2021-06-13 | 39.809786 |
In [127]:
_ = plt.figure(figsize=(15, 5)) _ = sns.lineplot(data=pca_df.reset_index(), x="date", y="windowed_PCA", legend=True) _ = plt.title("Windowed PCA") _ = plt.ylabel("Explained Variance of\nfirst principal component")
In [128]:
tmp = apply_scaling(pca_df, method="Standard")
In [129]:
significance_level=0.1 sig_plot_df = tmp.mask( tmp > norm.ppf(1 - significance_level), 1 ).mask(tmp <= norm.ppf(1 - significance_level), 0)
In [130]:
norm.ppf(1 - significance_level) norm.ppf(1 - significance_level)
Out[130]:
1.2815515655446004
In [131]:
# _ = sns.heatmap(sig_plot_df.T)
In [132]:
_ = plt.figure(figsize=(15, 5)) _ = plt.plot(tmp) _ = plt.title("Windowed PCA (significant peaks in yellow based on z-test)") _ = plt.ylabel("Scaled explained variance of\nfirst principal component") for sig_date in sig_plot_df[sig_plot_df==True].dropna().index: _ = plt.axvline(sig_date, ls="--", c=JmsColors.YELLOW)
In [133]:
window_size=7 window_range = np.arange(0, len(tmp)-window_size, window_size) coefs_list = [] for window_begin in window_range: win_df = tmp.iloc[window_begin: window_begin + window_size].iloc[:, 0] mod = LinearRegression() _ = mod.fit(X=win_df.reset_index().index.to_numpy().reshape(-1, 1), y=win_df) coefs_list.append(mod.coef_[0])
In [134]:
_ = plt.plot(pd.Series(coefs_list))
In [135]:
fi_df = fluctuation_intensity(df.set_index("date").pipe(apply_scaling), win=28, xmin=0, xmax=1, col_first=1, col_last=df.set_index("date").shape[1])
In [136]:
_ = complexity_resonance_diagram(fi_df, plot_title="Fluctuation Intensity Plot", figsize = (20, 15))
In [137]:
ccp_df, sig_ccps_df = cumulative_complexity_peaks(df=fi_df, significant_level_item=0.05, significant_level_time=0.05)
In [138]:
_ = cumulative_complexity_peaks_plot(cumulative_complexity_peaks_df=ccp_df, significant_peaks_df=sig_ccps_df, plot_title="Cumulative Fluctuation Peaks Plot", figsize = (20, 15))
In [139]:
_ = plt.figure(figsize=(15, 5)) _ = plt.plot(fi_df.sum(axis=1).replace(0, np.nan)) _ = plt.title("Sum of Fluctuation Intensity (significant fluctuation peaks in yellow)") for sig_date in sig_ccps_df[sig_ccps_df==True].dropna().index: _ = plt.axvline(sig_date, ls="--", c=JmsColors.YELLOW)
In [140]:
print(f"time elapsed: {time.time()-start} seconds")
time elapsed: 115.71455883979797 seconds