chuyển đổi cấu trúc

2026-07-08 15:11:34 +07:00
parent b041ccc11a
commit b22d327c5f
45 changed files with 201 additions and 4 deletions
@@ -0,0 +1,14 @@
+,PC1,PC2,Uniqueness
+Temperature,0.7859143713690758,-0.15342458640921602,0.35879949716071136
+Noise,0.7392275315643728,-0.14504349077151615,0.43250504236205733
+Stormwind,0.7983317716636736,-0.15246305454979966,0.33942139934968485
+Respiratory,0.8141977005065351,-0.16078705379178215,0.31122962782282915
+Exercises,0.7861049487051457,-0.22688720348435373,0.3305612065163297
+Culture,0.7978897043499087,-0.05932506148429913,0.35985255677229944
+Beauty,0.8700711219058058,-0.1556431549883994,0.2187514511308295
+Education,0.8379709091142482,-0.2606539630840815,0.22986426700680263
+Society,0.8178019516337711,-0.097610506717968,0.3216721568822566
+Spirit,0.8499003710179431,-0.1980420960982776,0.23844868751656323
+Dirty,-0.12911798501226357,0.9075618096917915,0.15966010753533333
+Unsafe,-0.12502417318133396,0.8926423245694102,0.1875586365076436
+Danger,-0.19294042642503595,0.8920213064079019,0.1670719807652653
@@ -0,0 +1,4 @@
+Cluster,1,2,3,4,5,6,1_Pct,2_Pct,3_Pct,4_Pct,5_Pct,6_Pct
+1,23,4,7,9,0,4,48.93617021276596,8.51063829787234,14.893617021276595,19.148936170212767,0.0,8.51063829787234
+2,10,16,3,12,3,17,16.39344262295082,26.229508196721312,4.918032786885246,19.672131147540984,4.918032786885246,27.86885245901639
+3,43,25,37,72,0,22,21.608040201005025,12.562814070351758,18.592964824120603,36.18090452261307,0.0,11.055276381909549
@@ -0,0 +1,4 @@
+Cluster,1,2,3,4,5,6,1_Pct,2_Pct,3_Pct,4_Pct,5_Pct,6_Pct
+1,9,10,10,9,9,0,19.148936170212767,21.27659574468085,21.27659574468085,19.148936170212767,19.148936170212767,0.0
+2,16,13,11,11,10,0,26.229508196721312,21.311475409836063,18.0327868852459,18.0327868852459,16.39344262295082,0.0
+3,71,53,33,29,12,1,35.678391959798994,26.633165829145728,16.582914572864322,14.572864321608039,6.030150753768844,0.5025125628140703
@@ -0,0 +1,4 @@
+Cluster,1,2,3,4,5,1_Pct,2_Pct,3_Pct,4_Pct,5_Pct
+1,8,12,9,9,9,17.02127659574468,25.53191489361702,19.148936170212767,19.148936170212767,19.148936170212767
+2,21,18,6,9,7,34.42622950819672,29.508196721311474,9.836065573770492,14.754098360655737,11.475409836065573
+3,18,24,30,65,62,9.045226130653267,12.060301507537687,15.07537688442211,32.663316582914575,31.155778894472363
@@ -0,0 +1,4 @@
+Cluster,1,2,3,1_Pct,2_Pct,3_Pct
+1,18,27,2,38.297872340425535,57.446808510638306,4.25531914893617
+2,27,34,0,44.26229508196721,55.73770491803278,0.0
+3,91,108,0,45.7286432160804,54.2713567839196,0.0
@@ -0,0 +1,4 @@
+Cluster,1,2,3,4,5,6,1_Pct,2_Pct,3_Pct,4_Pct,5_Pct,6_Pct
+1,0,3,8,8,0,28,0.0,6.382978723404255,17.02127659574468,17.02127659574468,0.0,59.57446808510638
+2,0,7,12,25,4,13,0.0,11.475409836065573,19.672131147540984,40.98360655737705,6.557377049180328,21.311475409836063
+3,9,4,30,65,0,91,4.522613065326634,2.0100502512562812,15.07537688442211,32.663316582914575,0.0,45.7286432160804
@@ -0,0 +1,4 @@
+Cluster,1,2,3,4,5,1_Pct,2_Pct,3_Pct,4_Pct,5_Pct
+1,0,3,8,23,13,0.0,6.382978723404255,17.02127659574468,48.93617021276596,27.659574468085108
+2,6,9,15,17,14,9.836065573770492,14.754098360655737,24.59016393442623,27.86885245901639,22.950819672131146
+3,0,10,61,68,60,0.0,5.025125628140704,30.65326633165829,34.17085427135678,30.15075376884422
@@ -0,0 +1,4 @@
+Cluster,1,2,3,4,5,1_Pct,2_Pct,3_Pct,4_Pct,5_Pct
+1,8,10,21,6,2,17.02127659574468,21.27659574468085,44.680851063829785,12.76595744680851,4.25531914893617
+2,10,14,31,4,2,16.39344262295082,22.950819672131146,50.81967213114754,6.557377049180328,3.278688524590164
+3,32,37,123,6,1,16.08040201005025,18.592964824120603,61.80904522613066,3.015075376884422,0.5025125628140703
@@ -0,0 +1,239 @@
+import sys
+sys.stdout.reconfigure(encoding='utf-8')
+import pandas as pd
+import numpy as np
+import scipy.stats as stats
+from scipy.cluster.hierarchy import linkage, fcluster
+import matplotlib.pyplot as plt
+import seaborn as sns
+from factor_analyzer.factor_analyzer import calculate_kmo, calculate_bartlett_sphericity
+import os
+
+# --- Helper function for Varimax Rotation ---
+def varimax(loadings, max_iter=500, tolerance=1e-6):
+    X = loadings.copy()
+    n_rows, n_cols = X.shape
+    if n_cols < 2:
+        return X, np.eye(n_cols)
+    
+    R = np.eye(n_cols)
+    d = 0
+    for i in range(max_iter):
+        d_old = d
+        Lambda = np.dot(X, R)
+        grad = np.dot(X.T, Lambda**3 - (1.0 / n_rows) * np.dot(Lambda, np.diag(np.sum(Lambda**2, axis=0))))
+        u, s, vh = np.linalg.svd(grad)
+        R = np.dot(u, vh)
+        d = np.sum(s)
+        if d_old != 0 and (d - d_old) / d_old < tolerance:
+            break
+            
+    rotated_loadings = np.dot(X, R)
+    return rotated_loadings, R
+
+# --- Setup Paths ---
+base_dir = r"c:\Users\NASPC\Documents\Du án tại SG tháng 8"
+output_dir = os.path.join(base_dir, "PCA_HCA_Results")
+if not os.path.exists(output_dir):
+    os.makedirs(output_dir)
+
+file_path = os.path.join(base_dir, 'Data_VN_filter_v5.xlsx')
+df = pd.read_excel(file_path)
+
+# --- 1. Variables Definition ---
+pca_vars = [
+    'Temperature', 'Noise', 'Stormwind', 'Respiratory', # RES
+    'Exercises', 'Culture', 'Beauty', 'Education', 'Society', 'Spirit', # CES
+    'Dirty', 'Unsafe', 'Danger' # DES
+]
+demographic_vars = ['Gender', 'Career', 'Literacy']
+habit_vars = ['Distance', 'Frequency', 'Time', 'Transportation']
+
+# Ensure variables exist and drop NaNs for PCA
+df_pca = df[pca_vars].dropna()
+n_samples = len(df_pca)
+print(f"--- Bước 2: Phân tích thành phần chính (PCA) ---")
+print(f"Cỡ mẫu (N): {n_samples}")
+
+# --- 2. KMO & Bartlett's Test ---
+kmo_all, kmo_model = calculate_kmo(df_pca)
+bartlett_stat, bartlett_p = calculate_bartlett_sphericity(df_pca)
+
+print("\nBảng 1: Kaiser-Meyer-Olkin Test")
+print(f"Overall MSA: {kmo_model:.5f}")
+for var, kmo_val in zip(pca_vars, kmo_all):
+    print(f"  {var}: {kmo_val:.5f}")
+
+print("\nBảng 2: Bartlett's Test of Sphericity")
+df_bartlett = len(pca_vars) * (len(pca_vars) - 1) / 2
+print(f"X^2: {bartlett_stat:.5f}, df: {int(df_bartlett)}, p: {bartlett_p}")
+
+# --- 3. PCA with Varimax Rotation ---
+# Standardize with ddof=1 to match JASP
+X_mean = df_pca.mean()
+X_std = df_pca.std(ddof=1)
+X_scaled = (df_pca - X_mean) / X_std
+
+# Correlation matrix
+R_corr = np.corrcoef(df_pca.T)
+eigenvalues, eigenvectors = np.linalg.eigh(R_corr)
+
+# Sort descending
+idx = np.argsort(eigenvalues)[::-1]
+eigenvalues = eigenvalues[idx]
+eigenvectors = eigenvectors[:, idx]
+
+print("\nBảng 6: Component Characteristics (Unrotated)")
+for i, ev in enumerate(eigenvalues[:5]): # Print top 5 for illustration
+    print(f"Component {i+1}: Eigenvalue = {ev:.5f}, Proportion = {ev/sum(eigenvalues):.5f}")
+
+# Extract 2 components
+unrotated_loadings = eigenvectors[:, :2] * np.sqrt(eigenvalues[:2])
+rotated_loadings, R = varimax(unrotated_loadings)
+
+# Align signs with JASP: PC1 positive for Beauty, PC2 positive for Dirty
+idx_beauty = pca_vars.index('Beauty')
+idx_dirty = pca_vars.index('Dirty')
+
+if rotated_loadings[idx_beauty, 0] < 0:
+    rotated_loadings[:, 0] = -rotated_loadings[:, 0]
+    R[:, 0] = -R[:, 0]
+if rotated_loadings[idx_dirty, 1] < 0:
+    rotated_loadings[:, 1] = -rotated_loadings[:, 1]
+    R[:, 1] = -R[:, 1]
+
+uniqueness = 1 - np.sum(rotated_loadings**2, axis=1)
+
+print("\nBảng 5: Component Loadings (Varimax Rotated)")
+loadings_df = pd.DataFrame(rotated_loadings, index=pca_vars, columns=['PC1', 'PC2'])
+loadings_df['Uniqueness'] = uniqueness
+print(loadings_df.round(5).to_string())
+
+# Save Loadings to CSV
+loadings_df.to_csv(os.path.join(output_dir, 'PCA_Component_Loadings.csv'))
+
+# --- 4. Plot PCA Loading Biplot ---
+plt.figure(figsize=(10, 8))
+plt.scatter(rotated_loadings[:, 0], rotated_loadings[:, 1], color='blue', alpha=0.5)
+for i, txt in enumerate(pca_vars):
+    plt.annotate(txt, (rotated_loadings[i, 0], rotated_loadings[i, 1]), xytext=(5,5), textcoords='offset points')
+plt.axhline(0, color='black',linewidth=1, ls='--')
+plt.axvline(0, color='black',linewidth=1, ls='--')
+plt.xlabel('Component 1 (ESS)')
+plt.ylabel('Component 2 (DES)')
+plt.title('PCA Loading Plot (Varimax Rotated)')
+plt.grid(True, linestyle=':', alpha=0.6)
+plt.savefig(os.path.join(output_dir, 'PCA_Loading_Plot.png'), dpi=300)
+plt.close()
+
+# --- 5. Component Scores for HCA ---
+# Calculate Standardized Component Scores
+scores_std = (np.dot(X_scaled, eigenvectors[:, :2]) / np.sqrt(eigenvalues[:2])).dot(R)
+
+print("\n--- Bước 3: Phân tích cụm phân cấp (HCA) & Chi-square ---")
+# --- 6. HCA (Ward's Method) ---
+Z = linkage(scores_std, method='ward')
+clusters = fcluster(Z, 3, criterion='maxclust')
+
+# Map clusters to original dataframe
+# Note: Since we dropped NaNs for PCA, we need to carefully assign back
+# Assuming original df has no NaNs in these columns based on previous logs (N=307 valid).
+# Just to be safe, we assign via index.
+df.loc[df_pca.index, 'Cluster'] = clusters
+df.loc[df_pca.index, 'PC1_Score'] = scores_std[:, 0]
+df.loc[df_pca.index, 'PC2_Score'] = scores_std[:, 1]
+
+# Align cluster labels with JASP based on sizes (199, 61, 47) and Means
+cluster_sizes = df['Cluster'].value_counts()
+print(f"\nGiai đoạn 3.1 - Gom cụm (HCA)")
+print("Cluster Sizes (Before Label Alignment):")
+print(cluster_sizes)
+
+# Map our generic cluster IDs (1, 2, 3) to JASP's cluster IDs
+# From our previous test:
+# Our Cluster with size 199 -> JASP Cluster 3 (Hài hòa & Thụ hưởng)
+# Our Cluster with size 61 -> JASP Cluster 2 (Thờ ơ)
+# Our Cluster with size 47 -> JASP Cluster 1 (Thực dụng & Lo ngại)
+mapping = {}
+for clst, size in cluster_sizes.items():
+    if size == 199:
+        mapping[clst] = 3
+    elif size == 61:
+        mapping[clst] = 2
+    elif size == 47:
+        mapping[clst] = 1
+    else:
+        mapping[clst] = clst # fallback
+
+df['Cluster'] = df['Cluster'].map(mapping)
+print("\nCluster Sizes (Aligned with JASP):")
+print(df['Cluster'].value_counts())
+
+cluster_means = df.groupby('Cluster')[['PC1_Score', 'PC2_Score']].mean()
+print("\nBảng 3: Cluster Means")
+print(cluster_means.round(5))
+
+# Plot Cluster Means
+cluster_means.plot(kind='bar', figsize=(10, 6))
+plt.title('Cluster Means for PC1 and PC2')
+plt.ylabel('Mean Standardized Score')
+plt.xlabel('Cluster')
+plt.axhline(0, color='black', linewidth=0.8, ls='--')
+plt.xticks(rotation=0)
+plt.savefig(os.path.join(output_dir, 'HCA_Cluster_Means_Plot.png'), dpi=300)
+plt.close()
+
+# Save DataFrame with Clusters
+output_dataset_path = os.path.join(base_dir, 'Data_VN_filter_v5_with_clusters.xlsx')
+df.to_excel(output_dataset_path, index=False)
+print(f"\nĐã lưu dataset mới kèm nhãn Cụm và Điểm nhân tố tại: {output_dataset_path}")
+
+# --- 7. Profiling (Chi-Square) ---
+def compute_cramer_v(chi2, n, shape):
+    return np.sqrt(chi2 / (n * (min(shape) - 1)))
+
+def profile_clusters(df, variables, title):
+    print(f"\n{title}")
+    results = []
+    
+    for var in variables:
+        if var not in df.columns:
+            continue
+        # Contingency table (cross-tabulation)
+        ct = pd.crosstab(df['Cluster'], df[var])
+        
+        # Chi-square test
+        chi2, p, dof, expected = stats.chi2_contingency(ct)
+        n = ct.sum().sum()
+        
+        # Cramer's V and Phi
+        v = compute_cramer_v(chi2, n, ct.shape)
+        
+        print(f"\n--- {var} ---")
+        print("Contingency Table (Counts):")
+        print(ct)
+        print(f"Chi-Square: {chi2:.5f}, df: {dof}, p: {p:.5f}")
+        print(f"Cramer's V (Contingency coefficient approx): {v:.5f}")
+        
+        results.append({
+            'Variable': var,
+            'Chi_Square': chi2,
+            'df': dof,
+            'p_value': p,
+            'Cramers_V': v
+        })
+        
+        # Detailed Table (like JASP) with row percentages
+        ct_pct = ct.div(ct.sum(axis=1), axis=0) * 100
+        # save detailed to csv just in case
+        ct_detailed = pd.concat([ct, ct_pct.add_suffix('_Pct')], axis=1)
+        ct_detailed.to_csv(os.path.join(output_dir, f'Profiling_{var}_Contingency.csv'))
+
+    res_df = pd.DataFrame(results)
+    res_df.to_csv(os.path.join(output_dir, f'Profiling_{title.replace(" ", "_")}_Summary.csv'), index=False)
+
+profile_clusters(df, demographic_vars, "Lần so sánh 1: Nhân khẩu học (Demographics)")
+profile_clusters(df, habit_vars, "Lần so sánh 2: Thói quen tương tác (Habits)")
+
+print("\n--- HOÀN THÀNH ---")
+print(f"Các bảng kết quả và biểu đồ đã được lưu tại: {output_dir}")