matti_jms_collabs/covid19/uni_maryland_fb_covid19_windowed_PCA.ipynb

Windowed PCA and fluctuation intensity in COVID19 data¶

In [99]:

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 [100]:

start = time.time()

In [101]:

if "jms_style_sheet" in plt.style.available:
    plt.style.use("jms_style_sheet")

In [102]:

country_choice="Netherlands"
filepath=f"../shield-complexity/data/{country_choice}_worldsurvey_nonmissing_c_of_v_since_2021-06-08_to_2022-01-21.csv"

In [103]:

df = pd.read_csv(filepath)

In [104]:

df["date"] = pd.to_datetime(df["date"]).dt.date

In [105]:

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_concerned_sideeffects	pct_community_cli	pct_activity_work_outside_home	...	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_work_outside_home_1d	pct_mask_shop_1d	pct_mask_spent_time_1d
0	2021-06-09	14.535989	13.195395	0.770978	3.346972	6.422849	1.337626	1.292528	4.387491	1.258690	...	5.271337	4.509109	3.382587	3.528127	4.871057	1.119098	3.155650	1.074371	0.432519	3.148366
1	2021-06-10	14.773657	14.773657	0.715611	3.498316	7.144337	1.267841	1.550583	4.649098	1.376345	...	4.687061	3.669607	3.150838	3.858300	5.208682	1.071412	3.435815	1.034580	0.358411	3.208651
2	2021-06-11	18.572413	18.203818	0.649700	3.456743	9.709611	1.245488	1.449014	4.170443	1.320875	...	5.485248	3.985184	3.796074	3.540908	4.152605	1.011222	3.139859	0.967428	0.346627	2.720786
3	2021-06-13	17.418052	17.418052	0.703517	3.761218	10.277039	1.348971	1.526305	5.509964	1.503175	...	5.254944	3.342272	3.168605	3.342676	4.585373	1.245500	3.266020	1.129369	0.294228	3.224900
4	2021-06-14	20.332173	17.400103	0.716563	3.497908	7.931059	1.214978	1.482273	5.389113	1.490736	...	5.430036	3.805503	3.120780	3.441871	4.359385	1.047207	3.108143	0.911346	0.471995	3.524901

5 rows × 67 columns

Out[105]:

(215, 67)

In [106]:

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 215 entries, 0 to 214
Data columns (total 67 columns):
 #   Column                             Non-Null Count  Dtype  
---  ------                             --------------  -----  
 0   date                               215 non-null    object 
 1   pct_covid                          215 non-null    float64
 2   pct_flu                            215 non-null    float64
 3   percent_mc                         215 non-null    float64
 4   percent_hf                         215 non-null    float64
 5   percent_anos                       215 non-null    float64
 6   pct_twodoses                       215 non-null    float64
 7   pct_concerned_sideeffects          215 non-null    float64
 8   pct_community_cli                  215 non-null    float64
 9   pct_activity_work_outside_home     215 non-null    float64
 10  pct_activity_shop                  215 non-null    float64
 11  pct_activity_restaurant_bar        215 non-null    float64
 12  pct_activity_spent_time            215 non-null    float64
 13  pct_activity_large_event           215 non-null    float64
 14  pct_activity_public_transit        215 non-null    float64
 15  pct_food_security                  215 non-null    float64
 16  pct_anxious_7d                     215 non-null    float64
 17  pct_depressed_7d                   215 non-null    float64
 18  pct_symp_fever                     215 non-null    float64
 19  pct_symp_cough                     215 non-null    float64
 20  pct_symp_diff_breathing            215 non-null    float64
 21  pct_symp_fatigue                   215 non-null    float64
 22  pct_symp_stuffy_nose               215 non-null    float64
 23  pct_symp_aches                     215 non-null    float64
 24  pct_symp_sore_throat               215 non-null    float64
 25  pct_symp_chest_pain                215 non-null    float64
 26  pct_symp_nausea                    215 non-null    float64
 27  pct_symp_headache                  215 non-null    float64
 28  pct_symp_chills                    215 non-null    float64
 29  pct_testing_rate                   215 non-null    float64
 30  pct_avoid_contact                  215 non-null    float64
 31  pct_worried_catch_covid            215 non-null    float64
 32  pct_belief_distancing_effective    215 non-null    float64
 33  pct_belief_masking_effective       215 non-null    float64
 34  pct_others_distanced_public        215 non-null    float64
 35  pct_others_masked_public           215 non-null    float64
 36  pct_belief_children_immune         215 non-null    float64
 37  pct_belief_no_spread_hot_humid     215 non-null    float64
 38  pct_received_news_local_health     215 non-null    float64
 39  pct_received_news_experts          215 non-null    float64
 40  pct_received_news_who              215 non-null    float64
 41  pct_received_news_govt_health      215 non-null    float64
 42  pct_received_news_politicians      215 non-null    float64
 43  pct_received_news_journalists      215 non-null    float64
 44  pct_received_news_friends          215 non-null    float64
 45  pct_received_news_religious        215 non-null    float64
 46  pct_received_news_none             215 non-null    float64
 47  pct_trust_covid_info_local_health  215 non-null    float64
 48  pct_trust_covid_info_experts       215 non-null    float64
 49  pct_trust_covid_info_who           215 non-null    float64
 50  pct_trust_covid_info_govt_health   215 non-null    float64
 51  pct_trust_covid_info_politicians   215 non-null    float64
 52  pct_trust_covid_info_journalists   215 non-null    float64
 53  pct_trust_covid_info_friends       215 non-null    float64
 54  pct_trust_covid_info_religious     215 non-null    float64
 55  pct_want_info_covid_treatment      215 non-null    float64
 56  pct_want_info_covid_variants       215 non-null    float64
 57  pct_want_info_children_education   215 non-null    float64
 58  pct_want_info_economic_impact      215 non-null    float64
 59  pct_want_info_mental_health        215 non-null    float64
 60  pct_want_info_relationships        215 non-null    float64
 61  pct_want_info_employment           215 non-null    float64
 62  pct_want_info_none                 215 non-null    float64
 63  pct_delayed_care_cost              215 non-null    float64
 64  pct_mask_work_outside_home_1d      215 non-null    float64
 65  pct_mask_shop_1d                   215 non-null    float64
 66  pct_mask_spent_time_1d             215 non-null    float64
dtypes: float64(66), object(1)
memory usage: 112.7+ KB

In [107]:

shuffle_data=False
if shuffle_data:
    df = df.set_index("date").sample(frac=1, random_state=69420).set_index(df["date"])
else:
    df = df.set_index("date")

In [108]:

_ = plt.figure(figsize=(20, 10))
_ = sns.heatmap(df
# .set_index("date")
.T
)

In [109]:

_ = plt.figure(figsize=(20, 10))
_ = sns.heatmap(df
# .set_index("date")
.pipe(apply_scaling, method="Standard")
.T
)

In [110]:

df_melt = (df
# .set_index("date")
.pipe(apply_scaling, method="Standard")
.reset_index()
.melt(id_vars="date"))

In [111]:

_ = plt.figure(figsize=(30, 7))
_ = sns.lineplot(data=df_melt, x="date", y="value", hue="variable", legend=False)

In [112]:

g = sns.FacetGrid(df_melt, col="variable", col_wrap=3, aspect=2)
_ = g.map(sns.lineplot, "date", "value")

In [113]:

# 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 [114]:

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 [115]:

window=28
pca_df = window_PCA(x=df, window_size=window).set_index(df.reset_index().loc[:df.shape[0]-window-1, "date"])

In [116]:

pca_df.head()

Out[116]:

<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	35.991569
2021-06-10	37.060099
2021-06-11	38.635545
2021-06-13	38.527492
2021-06-14	38.777388

In [117]:

_ = plt.figure(figsize=(15, 5))
_ = sns.lineplot(data=pca_df.reset_index(), x="date", y="windowed_PCA", legend=True)
_ = plt.title(f"Windowed PCA: {country_choice}")
_ = plt.ylabel("Explained Variance of\nfirst principal component")

In [118]:

tmp = apply_scaling(pca_df, 
method="Standard")

In [119]:

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 [120]:

norm.ppf(1 - significance_level)
norm.ppf(1 - significance_level)

Out[120]:

1.2815515655446004

In [121]:

# _ = sns.heatmap(sig_plot_df.T)

In [122]:

_ = plt.figure(figsize=(15, 5))
_ = plt.plot(tmp)
_ = plt.title(f"Windowed PCA: {country_choice} (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 [123]:

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 [124]:

_ = plt.plot(pd.Series(coefs_list))

In [125]:

fi_df = fluctuation_intensity(df.pipe(apply_scaling, method="Standard"), 
win=28, 
xmin=0, 
xmax=1, 
col_first=1, 
col_last=df.shape[1])

In [126]:

_ = complexity_resonance_diagram(fi_df, plot_title=f"Fluctuation Intensity Plot: {country_choice}", figsize = (20, 15))

In [127]:

ccp_df, sig_ccps_df = cumulative_complexity_peaks(df=fi_df, significant_level_item=0.05, significant_level_time=0.05)

In [128]:

_ = cumulative_complexity_peaks_plot(cumulative_complexity_peaks_df=ccp_df, 
significant_peaks_df=sig_ccps_df, plot_title=f"Cumulative Fluctuation Peaks Plot: {country_choice}", figsize = (20, 15))

In [129]:

_ = plt.figure(figsize=(15, 5))
_ = plt.plot(fi_df.sum(axis=1).replace(0, np.nan))
_ = plt.title(f"Sum of Fluctuation Intensity: {country_choice} (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 [130]:

print(f"time elapsed: {time.time()-start} seconds")

time elapsed: 62.25389528274536 seconds

In [137]:

rand_df = pd.concat([pd.Series(np.random.normal(size=df.shape[0]), name=f"feat_{x}") for x in range(0, df.shape[1])], axis=1)
rand_df.head()

Out[137]:

<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>

	feat_0	feat_1	feat_2	feat_3	feat_4	feat_5	feat_6	feat_7	feat_8	feat_9	...	feat_56	feat_57	feat_58	feat_59	feat_60	feat_61	feat_62	feat_63	feat_64	feat_65
0	0.607970	1.973006	-0.584729	-0.306673	-0.732000	0.968532	1.185900	-1.497914	-1.186765	-0.898712	...	0.237570	0.423949	0.048828	0.858256	0.347837	-0.032485	-0.569684	0.363079	-0.621208	-0.467781
1	-0.555724	0.567568	0.318248	1.266914	-3.469426	-1.642858	-0.398362	-0.142489	-1.533189	-0.196776	...	-1.362928	1.869612	-0.300996	0.071474	0.527987	-0.737835	0.030544	3.015852	0.993087	-0.067536
2	-0.619340	-0.976191	0.576521	-0.619316	-0.072111	-0.837100	-0.414019	-0.028701	-0.605296	1.424502	...	-0.581591	0.562998	-0.179906	-0.140870	0.682107	-1.232831	0.225519	-0.788377	1.165745	-0.659170
3	-0.077220	-1.735391	0.234120	-0.645802	0.410166	0.727649	-0.317359	-0.263082	0.321844	-0.990930	...	1.617682	-1.190932	0.098023	0.425726	-0.851748	-1.134010	0.646297	-2.004089	0.460868	0.385467
4	0.616478	-0.623573	1.501336	0.635861	-0.231083	0.669580	-1.104844	-1.754440	1.210693	0.642050	...	0.385674	-0.648049	-0.168010	1.173464	0.289020	0.538613	0.940948	-1.086279	0.141587	0.714852

5 rows × 66 columns

In [143]:

pca_df = window_PCA(x=rand_df, window_size=window)

In [144]:

pca_df.head()

Out[144]:

<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
0	6.207433
1	6.082589
2	5.976907
3	5.947471
4	6.008262

In [145]:

_ = plt.figure(figsize=(15, 5))
_ = sns.lineplot(data=pca_df.reset_index(), x="index", y="windowed_PCA", legend=True)
_ = plt.title("Windowed PCA: random data")
_ = plt.ylabel("Explained Variance of\nfirst principal component")

In [ ]: