Exploratory multivariate technique for contingency tables. Decomposes the chi-squared distance between rows and columns of a frequency table into principal axes:
$$
\chi^2 = N \sum_{ij} \frac{(p_{ij} - p_{i+}p_{+j})^2}{p_{i+}p_{+j}}
$$
Biplot shows rows (types) and columns (attributes) in low-dimensional space, revealing typological patterns.
Step 3: Lithic Analysis — Reduction Sequence and Debitage
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import spearmanr, kruskal
# ------------------------------------------------------------------ ## Lithic assemblage analysis: reduction sequence from debitage# ------------------------------------------------------------------ #
np.random.seed(42)
n_lithics = 400# Flake attributes (simulated)
lithics_data = {
'artifact_id': [f"L{i:04d}"for i inrange(n_lithics)],
'site': np.random.choice(['Area 1', 'Area 2', 'Area 3'], n_lithics, p=[0.5, 0.3, 0.2]),
'material': np.random.choice(['Flint', 'Obsidian', 'Chert', 'Quartzite'], n_lithics,
p=[0.50, 0.20, 0.25, 0.05]),
'type': np.random.choice(['Primary flake', 'Secondary flake', 'Interior flake',
'Blade', 'Core', 'Biface', 'Tool'], n_lithics,
p=[0.15, 0.25, 0.30, 0.10, 0.05, 0.05, 0.10]),
'cortex_pct': np.random.choice([0, 10, 25, 50, 75, 100], n_lithics),
'length_mm': np.random.lognormal(3.2, 0.5, n_lithics),
'width_mm': np.random.lognormal(2.9, 0.5, n_lithics),
'thickness_mm': np.random.lognormal(1.8, 0.5, n_lithics),
'platform': np.random.choice(['cortical', 'plain', 'faceted', 'punctiform', 'absent'], n_lithics,
p=[0.15, 0.40, 0.25, 0.10, 0.10]),
'flake_scars': np.random.randint(0, 15, n_lithics),
'retouched': np.random.choice([True, False], n_lithics, p=[0.15, 0.85]),
}
lithics_df = pd.DataFrame(lithics_data)
# Compute reduction index (proxy: log weight ÷ dorsal scar count)
lithics_df['weight_g'] = (lithics_df['length_mm'] * lithics_df['width_mm'] * lithics_df['thickness_mm']) * 0.002
lithics_df['reduction_index'] = (lithics_df['flake_scars'] + 1) / (lithics_df['cortex_pct']/100 + 0.1)
print(f"Lithic assemblage: {len(lithics_df)} artifacts")
print("\nType distribution:")
print(lithics_df['type'].value_counts())
# ---- Reduction stage by cortex percentage ---------------------- #
cortex_by_type = lithics_df.groupby('type')['cortex_pct'].describe()
print("\nCortex % by artifact type:")
print(cortex_by_type[['mean','std','50%']].round(1))
# ---- Platform × flake type association ------------------------- #
ct_plat = pd.crosstab(lithics_df['type'], lithics_df['platform'])
from scipy.stats import chi2_contingency
chi2, p, _, _ = chi2_contingency(ct_plat)
print(f"\nPlatform × type: χ²={chi2:.2f}, p={p:.4f}")
# ---- Metric analysis by material -------------------------------- ## Kruskal-Wallis test for length differences by material
materials = lithics_df['material'].unique()
groups = [lithics_df[lithics_df['material']==m]['length_mm'].values for m in materials]
h_stat, p_kw = kruskal(*groups)
print(f"\nLength by material (Kruskal-Wallis): H={h_stat:.2f}, p={p_kw:.4f}")
# ---- Visualization --------------------------------------------- #
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
# Type distribution
type_counts = lithics_df['type'].value_counts()
colors_lit = ['#8B4513','#D2691E','#A0522D','#CD853F','#DEB887','#F4A460','#FFDEAD']
axes[0][0].barh(type_counts.index, type_counts.values, color=colors_lit[:len(type_counts)],
edgecolor='black', linewidth=0.7)
axes[0][0].set_xlabel("Count"); axes[0][0].set_title("Lithic Type Distribution")
axes[0][0].grid(axis='x', alpha=0.3)
# Reduction sequence: cortex by type
lithics_df.boxplot(column='cortex_pct', by='type', ax=axes[0][1])
axes[0][1].set_title("Cortex % by Artifact Type (Reduction Stage)")
axes[0][1].set_xlabel("Artifact type"); axes[0][1].set_ylabel("Cortex %")
plt.sca(axes[0][1]); plt.xticks(rotation=30, fontsize=7)
# Length distribution by materialfor mat in materials:
subset = lithics_df[lithics_df['material']==mat]['length_mm']
axes[1][0].hist(subset, bins=20, alpha=0.5, density=True, label=mat)
axes[1][0].set_xlabel("Length (mm)"); axes[1][0].set_ylabel("Density")
axes[1][0].set_title("Flake Length Distribution by Raw Material")
axes[1][0].legend(fontsize=8); axes[1][0].grid(True, alpha=0.3)
# Length × width scatter (flake shape)
axes[1][1].scatter(lithics_df['length_mm'], lithics_df['width_mm'],
c=['#3498db'if r else'#e74c3c'for r in lithics_df['retouched']],
s=15, alpha=0.5)
axes[1][1].set_xlabel("Length (mm)"); axes[1][1].set_ylabel("Width (mm)")
axes[1][1].set_title("Flake Dimensions (blue=unretouched, red=retouched)")
axes[1][1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("lithic_analysis.png", dpi=150)
plt.show()
Advanced Usage
Hierarchical Cluster Analysis of Assemblages
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import pdist
# Build site-level attribute matrix (ware proportions)
ct_site_ware = pd.crosstab(df['site'], df['ware'], normalize='index')
# Bray-Curtis dissimilarity (for compositional data)try:
from scipy.spatial.distance import braycurtis
bc_matrix = np.zeros((len(ct_site_ware), len(ct_site_ware)))
for i inrange(len(ct_site_ware)):
for j inrange(len(ct_site_ware)):
bc_matrix[i,j] = braycurtis(ct_site_ware.iloc[i], ct_site_ware.iloc[j])
except Exception:
bc_matrix = pdist(ct_site_ware.values, metric='euclidean')
from scipy.spatial.distance import squareform
try:
dist_condensed = squareform(bc_matrix)
except Exception:
dist_condensed = pdist(ct_site_ware.values, metric='euclidean')
Z = linkage(dist_condensed, method='ward')
fig, ax = plt.subplots(figsize=(8, 5))
dendrogram(Z, labels=ct_site_ware.index.tolist(), ax=ax)
ax.set_title("Assemblage Cluster Analysis (Ward linkage, Bray-Curtis dissimilarity)")
ax.set_xlabel("Site"); ax.set_ylabel("Dissimilarity")
ax.grid(axis='y', alpha=0.3)
plt.tight_layout(); plt.savefig("assemblage_cluster.png", dpi=150); plt.show()
Troubleshooting
CA axes flip sign unexpectedly
SVD sign is arbitrary. Multiply both row and column coordinates by -1 for a given axis if needed for interpretability — the relative positions remain the same.
Chi-square test invalid with small expected frequencies
Fix: Merge rare categories before testing:
# Merge rare wares (<10 counts) into "Other"
rare = ware_counts[ware_counts < 10].index
df['ware_merged'] = df['ware'].replace({w: 'Other'for w in rare})
Outlier sherds distort CA
Fix: Remove outlier rows/columns with very low frequencies before CA:
ct_filtered = ct_ca[ct_ca.sum(axis=1) >= 5] # Min 5 sherds per context
ct_filtered = ct_filtered.loc[:, ct_filtered.sum(axis=0) >= 5]
