matti_jms_collabs/covid19/uni_maryland_fb_covid19_windowed_PCA.ipynb

Windowed PCA and fluctuation intensity in COVID19 data¶

In [1]:

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

start = time.time()

In [3]:

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

In [4]:

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

df = pd.read_csv(filepath)

In [6]:

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

In [7]:

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

(215, 67)

In [8]:

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

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

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

In [11]:

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

In [12]:

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

In [13]:

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

In [14]:

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

In [15]:

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

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, svd_solver="full", **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 [17]:

pd.Series(decomposition.PCA(n_components=None, svd_solver="full").fit(df).explained_variance_).to_frame(name="explained_variance")

Out[17]:

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

	explained_variance
0	31.185597
1	10.071345
2	8.045120
3	3.212244
4	1.982784
...	...
61	0.002219
62	0.002029
63	0.001589
64	0.001458
65	0.001144

66 rows × 1 columns

In [18]:

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

In [19]:

pca_df.head()

Out[19]:

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

_ = 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 [21]:

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

In [22]:

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

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

Out[23]:

1.2815515655446004

In [24]:

# _ = sns.heatmap(sig_plot_df.T)

In [25]:

_ = 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 [26]:

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

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

In [28]:

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

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

In [30]:

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

In [31]:

_ = 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 [32]:

_ = 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 [33]:

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

time elapsed: 80.39891409873962 seconds

In [34]:

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

<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.507446	-1.524815	-0.109627	0.038661	1.190169	-0.588766	-1.790490	0.174861	0.484895	1.191590	...	0.058584	1.305392	0.600960	-0.021390	1.100561	0.033990	-0.978640	-1.719107	-0.654882	-1.423091
1	-0.346043	0.059419	0.654253	0.918169	-0.040968	1.379308	-0.084455	1.689244	-0.629499	1.772785	...	-0.174882	0.683060	-0.393272	0.149387	0.066825	2.535854	-0.054839	-2.147302	0.260677	-0.382394
2	0.541965	0.828549	1.061066	-0.147975	-0.047946	0.837137	-0.670662	-2.241606	0.453212	-1.841148	...	-0.363934	-0.387437	-1.033888	-1.180247	-0.055773	-1.288102	-0.290769	0.152713	-0.225541	0.222705
3	1.702650	0.459099	-2.694137	0.440802	0.311015	0.532554	-1.024340	0.840414	0.808548	0.030647	...	-1.043615	0.426505	0.816201	-1.196336	0.321254	0.021498	-1.376914	1.577555	0.018009	0.600904
4	-0.798074	-0.750160	1.148647	1.064333	0.751707	0.989130	0.289990	1.119342	-0.701148	0.456162	...	-0.731079	0.803282	-0.031749	-0.728007	0.463038	1.237527	-1.437393	0.159985	-2.195214	-0.673937

5 rows × 66 columns

In [35]:

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

In [36]:

pca_df.head()

Out[36]:

<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.541275
1	6.541521
2	6.093961
3	6.299040
4	6.262574

In [37]:

_ = 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")

Jms for loop

In [38]:

window_size=28

window_range = np.arange(0, len(df)-window_size)

window_range = np.arange(0, 3)


for window_begin in window_range:
    print(f"window_begin: {window_begin}, window_end: {window_begin + window_size}")
    # display(df.iloc[window_begin: window_begin + window_size, :].head())
    print(df.iloc[window_begin: window_begin + window_size, :].index.min())
    print(df.iloc[window_begin: window_begin + window_size, :].index.max())

window_begin: 0, window_end: 28
2021-06-09
2021-07-11
window_begin: 1, window_end: 29
2021-06-10
2021-07-12
window_begin: 2, window_end: 30
2021-06-11
2021-07-13

Matti for loop

In [39]:

# for(i in (win:nrow(mattidata_withdate))){
#   res.pca <- prcomp(
#     mattidata_withdate %>% 
#       # For e.g. window of 28, keep observations 1-28, 2-29, 3-30 etc.
#       dplyr::filter(row_number() >= i-win+1 &
#                       row_number() <= i) %>% 
#       dplyr::select(-all_of(omitted_from_pca)),
#     scale = TRUE)

In [40]:

win = 28
for i in range(win, len(df)):
    print(f"window_begin: {i-win+1}, window_end: {i}")

window_begin: 1, window_end: 28
window_begin: 2, window_end: 29
window_begin: 3, window_end: 30
window_begin: 4, window_end: 31
window_begin: 5, window_end: 32
window_begin: 6, window_end: 33
window_begin: 7, window_end: 34
window_begin: 8, window_end: 35
window_begin: 9, window_end: 36
window_begin: 10, window_end: 37
window_begin: 11, window_end: 38
window_begin: 12, window_end: 39
window_begin: 13, window_end: 40
window_begin: 14, window_end: 41
window_begin: 15, window_end: 42
window_begin: 16, window_end: 43
window_begin: 17, window_end: 44
window_begin: 18, window_end: 45
window_begin: 19, window_end: 46
window_begin: 20, window_end: 47
window_begin: 21, window_end: 48
window_begin: 22, window_end: 49
window_begin: 23, window_end: 50
window_begin: 24, window_end: 51
window_begin: 25, window_end: 52
window_begin: 26, window_end: 53
window_begin: 27, window_end: 54
window_begin: 28, window_end: 55
window_begin: 29, window_end: 56
window_begin: 30, window_end: 57
window_begin: 31, window_end: 58
window_begin: 32, window_end: 59
window_begin: 33, window_end: 60
window_begin: 34, window_end: 61
window_begin: 35, window_end: 62
window_begin: 36, window_end: 63
window_begin: 37, window_end: 64
window_begin: 38, window_end: 65
window_begin: 39, window_end: 66
window_begin: 40, window_end: 67
window_begin: 41, window_end: 68
window_begin: 42, window_end: 69
window_begin: 43, window_end: 70
window_begin: 44, window_end: 71
window_begin: 45, window_end: 72
window_begin: 46, window_end: 73
window_begin: 47, window_end: 74
window_begin: 48, window_end: 75
window_begin: 49, window_end: 76
window_begin: 50, window_end: 77
window_begin: 51, window_end: 78
window_begin: 52, window_end: 79
window_begin: 53, window_end: 80
window_begin: 54, window_end: 81
window_begin: 55, window_end: 82
window_begin: 56, window_end: 83
window_begin: 57, window_end: 84
window_begin: 58, window_end: 85
window_begin: 59, window_end: 86
window_begin: 60, window_end: 87
window_begin: 61, window_end: 88
window_begin: 62, window_end: 89
window_begin: 63, window_end: 90
window_begin: 64, window_end: 91
window_begin: 65, window_end: 92
window_begin: 66, window_end: 93
window_begin: 67, window_end: 94
window_begin: 68, window_end: 95
window_begin: 69, window_end: 96
window_begin: 70, window_end: 97
window_begin: 71, window_end: 98
window_begin: 72, window_end: 99
window_begin: 73, window_end: 100
window_begin: 74, window_end: 101
window_begin: 75, window_end: 102
window_begin: 76, window_end: 103
window_begin: 77, window_end: 104
window_begin: 78, window_end: 105
window_begin: 79, window_end: 106
window_begin: 80, window_end: 107
window_begin: 81, window_end: 108
window_begin: 82, window_end: 109
window_begin: 83, window_end: 110
window_begin: 84, window_end: 111
window_begin: 85, window_end: 112
window_begin: 86, window_end: 113
window_begin: 87, window_end: 114
window_begin: 88, window_end: 115
window_begin: 89, window_end: 116
window_begin: 90, window_end: 117
window_begin: 91, window_end: 118
window_begin: 92, window_end: 119
window_begin: 93, window_end: 120
window_begin: 94, window_end: 121
window_begin: 95, window_end: 122
window_begin: 96, window_end: 123
window_begin: 97, window_end: 124
window_begin: 98, window_end: 125
window_begin: 99, window_end: 126
window_begin: 100, window_end: 127
window_begin: 101, window_end: 128
window_begin: 102, window_end: 129
window_begin: 103, window_end: 130
window_begin: 104, window_end: 131
window_begin: 105, window_end: 132
window_begin: 106, window_end: 133
window_begin: 107, window_end: 134
window_begin: 108, window_end: 135
window_begin: 109, window_end: 136
window_begin: 110, window_end: 137
window_begin: 111, window_end: 138
window_begin: 112, window_end: 139
window_begin: 113, window_end: 140
window_begin: 114, window_end: 141
window_begin: 115, window_end: 142
window_begin: 116, window_end: 143
window_begin: 117, window_end: 144
window_begin: 118, window_end: 145
window_begin: 119, window_end: 146
window_begin: 120, window_end: 147
window_begin: 121, window_end: 148
window_begin: 122, window_end: 149
window_begin: 123, window_end: 150
window_begin: 124, window_end: 151
window_begin: 125, window_end: 152
window_begin: 126, window_end: 153
window_begin: 127, window_end: 154
window_begin: 128, window_end: 155
window_begin: 129, window_end: 156
window_begin: 130, window_end: 157
window_begin: 131, window_end: 158
window_begin: 132, window_end: 159
window_begin: 133, window_end: 160
window_begin: 134, window_end: 161
window_begin: 135, window_end: 162
window_begin: 136, window_end: 163
window_begin: 137, window_end: 164
window_begin: 138, window_end: 165
window_begin: 139, window_end: 166
window_begin: 140, window_end: 167
window_begin: 141, window_end: 168
window_begin: 142, window_end: 169
window_begin: 143, window_end: 170
window_begin: 144, window_end: 171
window_begin: 145, window_end: 172
window_begin: 146, window_end: 173
window_begin: 147, window_end: 174
window_begin: 148, window_end: 175
window_begin: 149, window_end: 176
window_begin: 150, window_end: 177
window_begin: 151, window_end: 178
window_begin: 152, window_end: 179
window_begin: 153, window_end: 180
window_begin: 154, window_end: 181
window_begin: 155, window_end: 182
window_begin: 156, window_end: 183
window_begin: 157, window_end: 184
window_begin: 158, window_end: 185
window_begin: 159, window_end: 186
window_begin: 160, window_end: 187
window_begin: 161, window_end: 188
window_begin: 162, window_end: 189
window_begin: 163, window_end: 190
window_begin: 164, window_end: 191
window_begin: 165, window_end: 192
window_begin: 166, window_end: 193
window_begin: 167, window_end: 194
window_begin: 168, window_end: 195
window_begin: 169, window_end: 196
window_begin: 170, window_end: 197
window_begin: 171, window_end: 198
window_begin: 172, window_end: 199
window_begin: 173, window_end: 200
window_begin: 174, window_end: 201
window_begin: 175, window_end: 202
window_begin: 176, window_end: 203
window_begin: 177, window_end: 204
window_begin: 178, window_end: 205
window_begin: 179, window_end: 206
window_begin: 180, window_end: 207
window_begin: 181, window_end: 208
window_begin: 182, window_end: 209
window_begin: 183, window_end: 210
window_begin: 184, window_end: 211
window_begin: 185, window_end: 212
window_begin: 186, window_end: 213
window_begin: 187, window_end: 214

In [ ]: