{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":112509,"databundleVersionId":14254895,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# =====================================================================\n# CSIRO PASTURE BIOMASS - COMPREHENSIVE EDA WITH DEEP INSIGHTS\n# =====================================================================\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom pathlib import Path\nimport warnings\nfrom scipy import stats\nfrom scipy.stats import pearsonr, spearmanr, normaltest, shapiro, anderson\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.ensemble import IsolationForest\nimport itertools\n\nwarnings.filterwarnings('ignore')\n\n# Enhanced plotting style\nsns.set_style(\"whitegrid\")\nplt.rcParams['figure.figsize'] = (12, 6)\nplt.rcParams['font.size'] = 10\n%matplotlib inline\n\n# Configuration\nCSV_PATH = \"/kaggle/input/csiro-biomass/train.csv\"\nEXPECTED_TARGETS = [\"Dry_Green_g\", \"Dry_Dead_g\", \"Dry_Clover_g\", \"GDM_g\", \"Dry_Total_g\"]\nTARGET_WEIGHTS = {\"Dry_Green_g\": 0.1, \"Dry_Dead_g\": 0.1, \"Dry_Clover_g\": 0.1, \"GDM_g\": 0.2, \"Dry_Total_g\": 0.5}\n\npd.set_option(\"display.max_columns\", 150)\npd.set_option(\"display.width\", 200)\npd.set_option(\"display.float_format\", '{:.4f}'.format)\n\nprint(\"=\"*80)\nprint(\"CSIRO PASTURE BIOMASS PREDICTION - EXTENSIVE EDA\")\nprint(\"=\"*80)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load data\ndf_long = pd.read_csv(CSV_PATH)\nprint(f\"\\n{'='*80}\")\nprint(\"1. DATA STRUCTURE OVERVIEW\")\nprint(f\"{'='*80}\")\nprint(f\"Shape: {df_long.shape}\")\nprint(f\"Memory usage: {df_long.memory_usage(deep=True).sum() / 1024**2:.2f} MB\")\nprint(f\"\\nFirst few rows:\")\ndisplay(df_long.head(10))\n\nprint(f\"\\n{'='*80}\")\nprint(\"2. COLUMN TYPES & MISSING VALUES\")\nprint(f\"{'='*80}\")\ninfo_df = pd.DataFrame({\n    'dtype': df_long.dtypes,\n    'non_null': df_long.count(),\n    'null_count': df_long.isna().sum(),\n    'null_pct': (df_long.isna().sum() / len(df_long) * 100).round(2),\n    'unique_values': df_long.nunique(),\n    'sample_value': [df_long[col].dropna().iloc[0] if df_long[col].dropna().size > 0 else None \n                     for col in df_long.columns]\n})\ndisplay(info_df)\n\n# Parse dates\nif \"Sampling_Date\" in df_long.columns:\n    df_long[\"Sampling_Date\"] = pd.to_datetime(df_long[\"Sampling_Date\"], errors='coerce')\n    df_long['Year'] = df_long['Sampling_Date'].dt.year\n    df_long['Month'] = df_long['Sampling_Date'].dt.month\n    df_long['Quarter'] = df_long['Sampling_Date'].dt.quarter\n    df_long['Season'] = df_long['Month'].map({12: 'Summer', 1: 'Summer', 2: 'Summer',\n                                               3: 'Autumn', 4: 'Autumn', 5: 'Autumn',\n                                               6: 'Winter', 7: 'Winter', 8: 'Winter',\n                                               9: 'Spring', 10: 'Spring', 11: 'Spring'})\n    print(\"\\n✓ Date features extracted: Year, Month, Quarter, Season\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"3. IMAGE-LEVEL DATA QUALITY CHECKS\")\nprint(f\"{'='*80}\")\n\n# Check rows per image\nn_unique_images = df_long['image_path'].nunique()\nrows_per_image = df_long.groupby('image_path').size()\n\nprint(f\"\\nUnique images: {n_unique_images}\")\nprint(f\"Expected total rows (5 per image): {n_unique_images * 5}\")\nprint(f\"Actual total rows: {len(df_long)}\")\nprint(f\"\\nRows per image distribution:\")\nprint(rows_per_image.value_counts().sort_index())\n\n# Find problematic images\nproblematic_images = rows_per_image[rows_per_image != 5]\nif len(problematic_images) > 0:\n    print(f\"\\n⚠️  WARNING: {len(problematic_images)} images don't have exactly 5 rows!\")\n    display(df_long[df_long['image_path'].isin(problematic_images.index)].sort_values('image_path'))\nelse:\n    print(\"\\n✓ All images have exactly 5 rows\")\n\n# Check target completeness\nprint(f\"\\n{'='*80}\")\nprint(\"4. TARGET COMPLETENESS CHECK\")\nprint(f\"{'='*80}\")\n\ndef check_targets_complete(group):\n    return set(group['target_name'].values) == set(EXPECTED_TARGETS)\n\ncomplete_mask = df_long.groupby('image_path').apply(check_targets_complete)\nincomplete_images = complete_mask[~complete_mask]\n\nif len(incomplete_images) > 0:\n    print(f\"⚠️  {len(incomplete_images)} images missing some targets!\")\n    for img in incomplete_images.index[:5]:\n        img_data = df_long[df_long['image_path'] == img]\n        missing = set(EXPECTED_TARGETS) - set(img_data['target_name'].values)\n        print(f\"  {img}: missing {missing}\")\nelse:\n    print(\"✓ All images have all 5 target types\")\n\n# Metadata consistency check\nprint(f\"\\n{'='*80}\")\nprint(\"5. METADATA CONSISTENCY ACROSS SAME IMAGE\")\nprint(f\"{'='*80}\")\n\nmeta_cols = ['Sampling_Date', 'State', 'Pre_GSHH_NDVI', 'Height_Ave_cm', 'Species']\nmeta_cols = [c for c in meta_cols if c in df_long.columns]\n\ninconsistencies = {}\nfor col in meta_cols:\n    unique_per_image = df_long.groupby('image_path')[col].nunique()\n    inconsistent = unique_per_image[unique_per_image > 1]\n    inconsistencies[col] = len(inconsistent)\n    \n    if len(inconsistent) > 0:\n        print(f\"\\n⚠️  {col}: {len(inconsistent)} images have inconsistent values\")\n        sample_img = inconsistent.index[0]\n        display(df_long[df_long['image_path'] == sample_img][['image_path', 'target_name', col]])\n    else:\n        print(f\"✓ {col}: All images have consistent values\")\n\n# Duplicate check\nprint(f\"\\n{'='*80}\")\nprint(\"6. DUPLICATE DETECTION\")\nprint(f\"{'='*80}\")\n\nduplicates = df_long.duplicated(subset=['image_path', 'target_name'], keep=False)\nn_duplicates = duplicates.sum()\n\nif n_duplicates > 0:\n    print(f\"⚠️  Found {n_duplicates} duplicate (image, target) pairs!\")\n    display(df_long[duplicates].sort_values(['image_path', 'target_name']))\nelse:\n    print(\"✓ No duplicate (image, target) pairs found\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"7. CONVERTING TO WIDE FORMAT & FEATURE ENGINEERING\")\nprint(f\"{'='*80}\")\n\n# Pivot to wide format\ndf = df_long.pivot_table(\n    index='image_path',\n    columns='target_name',\n    values='target',\n    aggfunc='mean'\n).reindex(columns=EXPECTED_TARGETS)\n\n# Add metadata\nmeta = df_long.drop_duplicates('image_path').set_index('image_path')\nmeta_cols_to_add = [c for c in ['Sampling_Date', 'State', 'Species', 'Pre_GSHH_NDVI', \n                                 'Height_Ave_cm', 'Year', 'Month', 'Quarter', 'Season'] \n                    if c in meta.columns]\ndf = df.join(meta[meta_cols_to_add]).reset_index()\n\nprint(f\"Wide format shape: {df.shape}\")\nprint(f\"Columns: {df.columns.tolist()}\")\n\n# FEATURE ENGINEERING: Derived features\nprint(\"\\n🔧 Creating derived features...\")\n\n# Ratios and percentages\ndf['Green_Pct'] = (df['Dry_Green_g'] / df['Dry_Total_g'] * 100).replace([np.inf, -np.inf], np.nan)\ndf['Dead_Pct'] = (df['Dry_Dead_g'] / df['Dry_Total_g'] * 100).replace([np.inf, -np.inf], np.nan)\ndf['Clover_Pct'] = (df['Dry_Clover_g'] / df['Dry_Total_g'] * 100).replace([np.inf, -np.inf], np.nan)\ndf['GDM_Pct'] = (df['GDM_g'] / df['Dry_Total_g'] * 100).replace([np.inf, -np.inf], np.nan)\n\n# Ratios between components\ndf['Green_Dead_Ratio'] = (df['Dry_Green_g'] / (df['Dry_Dead_g'] + 1)).replace([np.inf, -np.inf], np.nan)\ndf['Green_Clover_Ratio'] = (df['Dry_Green_g'] / (df['Dry_Clover_g'] + 1)).replace([np.inf, -np.inf], np.nan)\n\n# Sum check (should equal Dry_Total_g)\ndf['Component_Sum'] = df['Dry_Green_g'] + df['Dry_Dead_g'] + df['Dry_Clover_g']\ndf['Sum_Error'] = df['Component_Sum'] - df['Dry_Total_g']\ndf['Sum_Error_Pct'] = (df['Sum_Error'] / df['Dry_Total_g'] * 100).replace([np.inf, -np.inf], np.nan)\n\n# NDVI per height (if available)\nif 'Pre_GSHH_NDVI' in df.columns and 'Height_Ave_cm' in df.columns:\n    df['NDVI_per_Height'] = (df['Pre_GSHH_NDVI'] / (df['Height_Ave_cm'] + 1)).replace([np.inf, -np.inf], np.nan)\n\n# Biomass per height\nif 'Height_Ave_cm' in df.columns:\n    df['Biomass_per_Height'] = (df['Dry_Total_g'] / (df['Height_Ave_cm'] + 1)).replace([np.inf, -np.inf], np.nan)\n\nprint(f\"✓ Created {df.shape[1] - len(EXPECTED_TARGETS) - len(meta_cols_to_add) - 1} derived features\")\nprint(f\"\\nNew shape: {df.shape}\")\ndisplay(df.head())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"8. PHYSICAL CONSTRAINT VALIDATION\")\nprint(f\"{'='*80}\")\n\n# Check if components sum to total\nprint(\"\\n📐 Component Sum Validation:\")\nprint(f\"  Mean sum error: {df['Sum_Error'].mean():.4f} g\")\nprint(f\"  Median sum error: {df['Sum_Error'].median():.4f} g\")\nprint(f\"  Std sum error: {df['Sum_Error'].std():.4f} g\")\nprint(f\"  Max absolute error: {df['Sum_Error'].abs().max():.4f} g\")\n\n# Percentage error distribution\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\naxes[0].hist(df['Sum_Error'].dropna(), bins=50, edgecolor='black')\naxes[0].axvline(0, color='red', linestyle='--', label='Perfect Match')\naxes[0].set_xlabel('Component Sum - Dry_Total_g (grams)')\naxes[0].set_ylabel('Frequency')\naxes[0].set_title('Distribution of Sum Error')\naxes[0].legend()\n\naxes[1].hist(df['Sum_Error_Pct'].dropna(), bins=50, edgecolor='black')\naxes[1].axvline(0, color='red', linestyle='--', label='Perfect Match')\naxes[1].set_xlabel('Sum Error (%)')\naxes[1].set_ylabel('Frequency')\naxes[1].set_title('Distribution of Sum Error (Percentage)')\naxes[1].legend()\n\nplt.tight_layout()\nplt.show()\n\n# Identify problematic samples\nlarge_errors = df[df['Sum_Error'].abs() > 10].copy()\nif len(large_errors) > 0:\n    print(f\"\\n⚠️  {len(large_errors)} samples have |sum_error| > 10g:\")\n    display(large_errors[['image_path', 'Dry_Green_g', 'Dry_Dead_g', 'Dry_Clover_g', \n                          'Component_Sum', 'Dry_Total_g', 'Sum_Error']].head(10))\n\n# Check for negative values\nprint(f\"\\n🔍 Negative Value Check:\")\nfor col in EXPECTED_TARGETS:\n    n_negative = (df[col] < 0).sum()\n    if n_negative > 0:\n        print(f\"  ⚠️  {col}: {n_negative} negative values\")\n    else:\n        print(f\"  ✓ {col}: No negative values\")\n\n# Check GDM vs Green\ndf['GDM_Green_Diff'] = df['GDM_g'] - df['Dry_Green_g']\nprint(f\"\\n📊 GDM vs Dry_Green relationship:\")\nprint(f\"  GDM > Green in {(df['GDM_Green_Diff'] > 0).sum()} samples\")\nprint(f\"  GDM = Green in {(df['GDM_Green_Diff'].abs() < 0.1).sum()} samples\")\nprint(f\"  GDM < Green in {(df['GDM_Green_Diff'] < 0).sum()} samples\")\nprint(f\"  Mean difference: {df['GDM_Green_Diff'].mean():.4f} g\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"9. COMPREHENSIVE STATISTICAL ANALYSIS OF TARGETS\")\nprint(f\"{'='*80}\")\n\n# Extended descriptive statistics\nextended_stats = []\nfor target in EXPECTED_TARGETS:\n    s = df[target].dropna()\n    if len(s) > 0:\n        stats_dict = {\n            'Target': target,\n            'Weight': TARGET_WEIGHTS[target],\n            'Count': len(s),\n            'Missing': df[target].isna().sum(),\n            'Mean': s.mean(),\n            'Median': s.median(),\n            'Std': s.std(),\n            'Min': s.min(),\n            'Max': s.max(),\n            'Q1': s.quantile(0.25),\n            'Q3': s.quantile(0.75),\n            'IQR': s.quantile(0.75) - s.quantile(0.25),\n            'Range': s.max() - s.min(),\n            'CV': (s.std() / s.mean() * 100) if s.mean() != 0 else np.nan,  # Coefficient of variation\n            'Skewness': s.skew(),\n            'Kurtosis': s.kurtosis(),\n        }\n        \n        # Outliers (IQR method)\n        Q1, Q3 = s.quantile(0.25), s.quantile(0.75)\n        IQR = Q3 - Q1\n        lower_bound = Q1 - 1.5 * IQR\n        upper_bound = Q3 + 1.5 * IQR\n        outliers = ((s < lower_bound) | (s > upper_bound)).sum()\n        stats_dict['Outliers_IQR'] = outliers\n        stats_dict['Outliers_Pct'] = (outliers / len(s) * 100)\n        \n        # Zeros\n        stats_dict['Zeros'] = (s == 0).sum()\n        stats_dict['Zeros_Pct'] = (s == 0).sum() / len(s) * 100\n        \n        extended_stats.append(stats_dict)\n\nextended_stats_df = pd.DataFrame(extended_stats)\ndisplay(extended_stats_df.round(4))\n\n# Normality tests\nprint(f\"\\n{'='*80}\")\nprint(\"10. NORMALITY TESTS FOR TARGETS\")\nprint(f\"{'='*80}\")\n\nnormality_results = []\nfor target in EXPECTED_TARGETS:\n    s = df[target].dropna()\n    if len(s) > 3:\n        # Shapiro-Wilk test (good for n < 5000)\n        if len(s) < 5000:\n            shapiro_stat, shapiro_p = shapiro(s)\n        else:\n            shapiro_stat, shapiro_p = np.nan, np.nan\n        \n        # D'Agostino-Pearson test\n        k2_stat, k2_p = normaltest(s)\n        \n        # Anderson-Darling test\n        anderson_result = anderson(s)\n        \n        normality_results.append({\n            'Target': target,\n            'Shapiro_Stat': shapiro_stat,\n            'Shapiro_p': shapiro_p,\n            'Shapiro_Normal': 'Yes' if shapiro_p > 0.05 else 'No',\n            'K2_Stat': k2_stat,\n            'K2_p': k2_p,\n            'K2_Normal': 'Yes' if k2_p > 0.05 else 'No',\n            'Anderson_Stat': anderson_result.statistic,\n            'Anderson_Critical_5%': anderson_result.critical_values[2],\n            'Anderson_Normal': 'Yes' if anderson_result.statistic < anderson_result.critical_values[2] else 'No'\n        })\n\nnormality_df = pd.DataFrame(normality_results)\ndisplay(normality_df.round(4))\n\nprint(\"\\n💡 Interpretation:\")\nprint(\"  • p-value > 0.05: Data is likely normally distributed\")\nprint(\"  • p-value < 0.05: Data is NOT normally distributed (reject null hypothesis)\")\nprint(\"  • Most biomass data is right-skewed → Consider log transformation\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"11. ADVANCED DISTRIBUTION VISUALIZATIONS\")\nprint(f\"{'='*80}\")\n\n# Multi-panel distribution analysis\nfig, axes = plt.subplots(5, 4, figsize=(20, 20))\n\nfor idx, target in enumerate(EXPECTED_TARGETS):\n    s = df[target].dropna()\n    \n    # Histogram\n    axes[idx, 0].hist(s, bins=50, edgecolor='black', alpha=0.7)\n    axes[idx, 0].axvline(s.mean(), color='red', linestyle='--', label=f'Mean: {s.mean():.1f}')\n    axes[idx, 0].axvline(s.median(), color='green', linestyle='--', label=f'Median: {s.median():.1f}')\n    axes[idx, 0].set_title(f'{target} - Histogram')\n    axes[idx, 0].set_xlabel('Value (g)')\n    axes[idx, 0].legend()\n    \n    # Box plot\n    bp = axes[idx, 1].boxplot(s, vert=True, patch_artist=True, showfliers=True)\n    bp['boxes'][0].set_facecolor('lightblue')\n    axes[idx, 1].set_title(f'{target} - Boxplot')\n    axes[idx, 1].set_ylabel('Value (g)')\n    \n    # Violin plot\n    parts = axes[idx, 2].violinplot([s], vert=True, showmeans=True, showmedians=True)\n    axes[idx, 2].set_title(f'{target} - Violin Plot')\n    axes[idx, 2].set_ylabel('Value (g)')\n    \n    # Q-Q plot for normality\n    stats.probplot(s, dist=\"norm\", plot=axes[idx, 3])\n    axes[idx, 3].set_title(f'{target} - Q-Q Plot')\n    axes[idx, 3].grid(True)\n\nplt.tight_layout()\nplt.show()\n\n# Log-transformed distributions (for skewed data)\nprint(\"\\n📊 Log-Transformed Distributions (for skewed targets):\")\nfig, axes = plt.subplots(2, 3, figsize=(18, 10))\naxes = axes.flatten()\n\nfor idx, target in enumerate(EXPECTED_TARGETS):\n    s = df[target].dropna()\n    if s.min() > 0 and abs(s.skew()) > 1:\n        s_log = np.log1p(s)\n        \n        axes[idx].hist(s_log, bins=50, edgecolor='black', alpha=0.7, color='green')\n        axes[idx].axvline(s_log.mean(), color='red', linestyle='--', label=f'Mean: {s_log.mean():.2f}')\n        axes[idx].axvline(s_log.median(), color='blue', linestyle='--', label=f'Median: {s_log.median():.2f}')\n        axes[idx].set_title(f'log1p({target}) - Skew: {s_log.skew():.2f}')\n        axes[idx].set_xlabel('log1p(Value)')\n        axes[idx].legend()\n    else:\n        axes[idx].text(0.5, 0.5, 'Not heavily skewed\\nor contains zeros', \n                       ha='center', va='center', transform=axes[idx].transAxes)\n        axes[idx].set_title(f'{target}')\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"12. CORRELATION ANALYSIS WITH STATISTICAL SIGNIFICANCE\")\nprint(f\"{'='*80}\")\n\n# Pearson correlation\ncorr_pearson = df[EXPECTED_TARGETS].corr(method='pearson')\n\n# Spearman correlation (for non-linear relationships)\ncorr_spearman = df[EXPECTED_TARGETS].corr(method='spearman')\n\n# Calculate p-values for correlations\ndef correlation_with_pvalue(df, method='pearson'):\n    \"\"\"Calculate correlation matrix with p-values\"\"\"\n    cols = df.columns\n    n = len(cols)\n    corr_matrix = np.zeros((n, n))\n    p_matrix = np.zeros((n, n))\n    \n    for i in range(n):\n        for j in range(n):\n            if i == j:\n                corr_matrix[i, j] = 1.0\n                p_matrix[i, j] = 0.0\n            else:\n                s1 = df.iloc[:, i].dropna()\n                s2 = df.iloc[:, j].dropna()\n                mask = df.iloc[:, i].notna() & df.iloc[:, j].notna()\n                if mask.sum() > 3:\n                    if method == 'pearson':\n                        corr, pval = pearsonr(df.iloc[:, i][mask], df.iloc[:, j][mask])\n                    else:\n                        corr, pval = spearmanr(df.iloc[:, i][mask], df.iloc[:, j][mask])\n                    corr_matrix[i, j] = corr\n                    p_matrix[i, j] = pval\n                else:\n                    corr_matrix[i, j] = np.nan\n                    p_matrix[i, j] = np.nan\n    \n    return pd.DataFrame(corr_matrix, index=cols, columns=cols), \\\n           pd.DataFrame(p_matrix, index=cols, columns=cols)\n\ncorr_pearson, p_pearson = correlation_with_pvalue(df[EXPECTED_TARGETS], method='pearson')\ncorr_spearman, p_spearman = correlation_with_pvalue(df[EXPECTED_TARGETS], method='spearman')\n\n# Visualization\nfig, axes = plt.subplots(1, 2, figsize=(16, 7))\n\n# Pearson correlation heatmap\nsns.heatmap(corr_pearson, annot=True, fmt='.3f', cmap='RdBu_r', center=0, \n            vmin=-1, vmax=1, square=True, ax=axes[0], cbar_kws={'label': 'Correlation'})\naxes[0].set_title('Pearson Correlation Matrix\\n(Linear Relationships)', fontsize=14, fontweight='bold')\n\n# Spearman correlation heatmap\nsns.heatmap(corr_spearman, annot=True, fmt='.3f', cmap='RdBu_r', center=0, \n            vmin=-1, vmax=1, square=True, ax=axes[1], cbar_kws={'label': 'Correlation'})\naxes[1].set_title('Spearman Correlation Matrix\\n(Monotonic Relationships)', fontsize=14, fontweight='bold')\n\nplt.tight_layout()\nplt.show()\n\n# Print significant correlations\nprint(\"\\n🔍 Significant Correlations (p < 0.05, |r| > 0.3):\")\nprint(\"\\nPearson:\")\nfor i in range(len(EXPECTED_TARGETS)):\n    for j in range(i+1, len(EXPECTED_TARGETS)):\n        r = corr_pearson.iloc[i, j]\n        p = p_pearson.iloc[i, j]\n        if abs(r) > 0.3 and p < 0.05:\n            print(f\"  {EXPECTED_TARGETS[i]} vs {EXPECTED_TARGETS[j]}: r={r:.3f}, p={p:.4f}\")\n\nprint(\"\\nSpearman:\")\nfor i in range(len(EXPECTED_TARGETS)):\n    for j in range(i+1, len(EXPECTED_TARGETS)):\n        r = corr_spearman.iloc[i, j]\n        p = p_spearman.iloc[i, j]\n        if abs(r) > 0.3 and p < 0.05:\n            print(f\"  {EXPECTED_TARGETS[i]} vs {EXPECTED_TARGETS[j]}: r={r:.3f}, p={p:.4f}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"12. CORRELATION ANALYSIS WITH STATISTICAL SIGNIFICANCE\")\nprint(f\"{'='*80}\")\n\n# Pearson correlation\ncorr_pearson = df[EXPECTED_TARGETS].corr(method='pearson')\n\n# Spearman correlation (for non-linear relationships)\ncorr_spearman = df[EXPECTED_TARGETS].corr(method='spearman')\n\n# Calculate p-values for correlations\ndef correlation_with_pvalue(df, method='pearson'):\n    \"\"\"Calculate correlation matrix with p-values\"\"\"\n    cols = df.columns\n    n = len(cols)\n    corr_matrix = np.zeros((n, n))\n    p_matrix = np.zeros((n, n))\n    \n    for i in range(n):\n        for j in range(n):\n            if i == j:\n                corr_matrix[i, j] = 1.0\n                p_matrix[i, j] = 0.0\n            else:\n                s1 = df.iloc[:, i].dropna()\n                s2 = df.iloc[:, j].dropna()\n                mask = df.iloc[:, i].notna() & df.iloc[:, j].notna()\n                if mask.sum() > 3:\n                    if method == 'pearson':\n                        corr, pval = pearsonr(df.iloc[:, i][mask], df.iloc[:, j][mask])\n                    else:\n                        corr, pval = spearmanr(df.iloc[:, i][mask], df.iloc[:, j][mask])\n                    corr_matrix[i, j] = corr\n                    p_matrix[i, j] = pval\n                else:\n                    corr_matrix[i, j] = np.nan\n                    p_matrix[i, j] = np.nan\n    \n    return pd.DataFrame(corr_matrix, index=cols, columns=cols), \\\n           pd.DataFrame(p_matrix, index=cols, columns=cols)\n\ncorr_pearson, p_pearson = correlation_with_pvalue(df[EXPECTED_TARGETS], method='pearson')\ncorr_spearman, p_spearman = correlation_with_pvalue(df[EXPECTED_TARGETS], method='spearman')\n\n# Visualization\nfig, axes = plt.subplots(1, 2, figsize=(16, 7))\n\n# Pearson correlation heatmap\nsns.heatmap(corr_pearson, annot=True, fmt='.3f', cmap='RdBu_r', center=0, \n            vmin=-1, vmax=1, square=True, ax=axes[0], cbar_kws={'label': 'Correlation'})\naxes[0].set_title('Pearson Correlation Matrix\\n(Linear Relationships)', fontsize=14, fontweight='bold')\n\n# Spearman correlation heatmap\nsns.heatmap(corr_spearman, annot=True, fmt='.3f', cmap='RdBu_r', center=0, \n            vmin=-1, vmax=1, square=True, ax=axes[1], cbar_kws={'label': 'Correlation'})\naxes[1].set_title('Spearman Correlation Matrix\\n(Monotonic Relationships)', fontsize=14, fontweight='bold')\n\nplt.tight_layout()\nplt.show()\n\n# Print significant correlations\nprint(\"\\n🔍 Significant Correlations (p < 0.05, |r| > 0.3):\")\nprint(\"\\nPearson:\")\nfor i in range(len(EXPECTED_TARGETS)):\n    for j in range(i+1, len(EXPECTED_TARGETS)):\n        r = corr_pearson.iloc[i, j]\n        p = p_pearson.iloc[i, j]\n        if abs(r) > 0.3 and p < 0.05:\n            print(f\"  {EXPECTED_TARGETS[i]} vs {EXPECTED_TARGETS[j]}: r={r:.3f}, p={p:.4f}\")\n\nprint(\"\\nSpearman:\")\nfor i in range(len(EXPECTED_TARGETS)):\n    for j in range(i+1, len(EXPECTED_TARGETS)):\n        r = corr_spearman.iloc[i, j]\n        p = p_spearman.iloc[i, j]\n        if abs(r) > 0.3 and p < 0.05:\n            print(f\"  {EXPECTED_TARGETS[i]} vs {EXPECTED_TARGETS[j]}: r={r:.3f}, p={p:.4f}\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"13. PAIRWISE RELATIONSHIPS ANALYSIS\")\nprint(f\"{'='*80}\")\n\n# Enhanced scatter matrix with regression lines\nfrom pandas.plotting import scatter_matrix\n\nfig = plt.figure(figsize=(18, 18))\naxes = scatter_matrix(df[EXPECTED_TARGETS].dropna(), figsize=(18, 18), \n                     alpha=0.3, diagonal='kde', density_kwds={'color': 'blue'})\n\n# Add regression lines\nfor i in range(len(EXPECTED_TARGETS)):\n    for j in range(len(EXPECTED_TARGETS)):\n        if i != j:\n            ax = axes[i, j]\n            x_data = df[EXPECTED_TARGETS[j]].dropna()\n            y_data = df[EXPECTED_TARGETS[i]].dropna()\n            mask = df[EXPECTED_TARGETS[j]].notna() & df[EXPECTED_TARGETS[i]].notna()\n            \n            if mask.sum() > 1:\n                x = df.loc[mask, EXPECTED_TARGETS[j]]\n                y = df.loc[mask, EXPECTED_TARGETS[i]]\n                z = np.polyfit(x, y, 1)\n                p = np.poly1d(z)\n                ax.plot(x.sort_values(), p(x.sort_values()), \"r--\", alpha=0.8, linewidth=2)\n\nplt.suptitle('Scatter Matrix of 5 Targets with Regression Lines', \n             y=1.0, fontsize=16, fontweight='bold')\nplt.tight_layout()\nplt.show()\n\n# Joint plots for key relationships\nprint(\"\\n📊 Detailed Joint Plots for Key Relationships:\")\n\nkey_pairs = [\n    ('Dry_Total_g', 'Dry_Green_g'),\n    ('Dry_Total_g', 'GDM_g'),\n    ('Dry_Green_g', 'GDM_g'),\n]\n\nfor x_col, y_col in key_pairs:\n    g = sns.jointplot(data=df, x=x_col, y=y_col, kind='reg', height=8,\n                      marginal_kws=dict(bins=30, fill=True))\n    \n    # Add correlation info\n    mask = df[x_col].notna() & df[y_col].notna()\n    if mask.sum() > 0:\n        r_pearson, p_pearson = pearsonr(df.loc[mask, x_col], df.loc[mask, y_col])\n        r_spearman, p_spearman = spearmanr(df.loc[mask, x_col], df.loc[mask, y_col])\n        \n        g.fig.suptitle(f'{x_col} vs {y_col}\\n' + \n                       f'Pearson r={r_pearson:.3f} (p={p_pearson:.4f}), ' +\n                       f'Spearman ρ={r_spearman:.3f} (p={p_spearman:.4f})',\n                       y=1.02, fontsize=12, fontweight='bold')\n    plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(f\"\\n{'='*80}\")\nprint(\"14. FEATURE RELATIONSHIPS - NDVI & HEIGHT ANALYSIS\")\nprint(f\"{'='*80}\")\n\nif 'Pre_GSHH_NDVI' in df.columns and 'Height_Ave_cm' in df.columns:\n    \n    # NDVI distribution and stats\n    print(\"\\n📊 NDVI Statistics:\")\n    ndvi_stats = df['Pre_GSHH_NDVI'].describe(percentiles=[.01, .05, .25, .5, .75, .95, .99])\n    display(ndvi_stats.to_frame().T)\n    \n    # Height distribution and stats\n    print(\"\\n📊 Height Statistics:\")\n    height_stats = df['Height_Ave_cm'].describe(percentiles=[.01, .05, .25, .5, .75, .95, .99])\n    display(height_stats.to_frame().T)\n    \n    # Correlations with all targets\n    print(\"\\n🔗 Correlations with Targets:\")\n    feature_corr = pd.DataFrame({\n        'NDVI_Pearson': [df[['Pre_GSHH_NDVI', t]].corr().iloc[0, 1] for t in EXPECTED_TARGETS],\n        'NDVI_Spearman': [df[['Pre_GSHH_NDVI', t]].corr(method='spearman').iloc[0, 1] for t in EXPECTED_TARGETS],\n        'Height_Pearson': [df[['Height_Ave_cm', t]].corr().iloc[0, 1] for t in EXPECTED_TARGETS],\n        'Height_Spearman': [df[['Height_Ave_cm', t]].corr(method='spearman').iloc[0, 1] for t in EXPECTED_TARGETS],\n    }, index=EXPECTED_TARGETS)\n    display(feature_corr.round(3))\n    \n    # Scatter plots: NDVI vs all targets\n    fig, axes = plt.subplots(2, 3, figsize=(18, 10))\n    axes = axes.flatten()\n    \n    for idx, target in enumerate(EXPECTED_TARGETS):\n        mask = df['Pre_GSHH_NDVI'].notna() & df[target].notna()\n        if mask.sum() > 0:\n            axes[idx].scatter(df.loc[mask, 'Pre_GSHH_NDVI'], df.loc[mask, target], \n                            alpha=0.4, s=20)\n            \n            # Add regression line\n            x = df.loc[mask, 'Pre_GSHH_NDVI']\n            y = df.loc[mask, target]\n            z = np.polyfit(x, y, 1)\n            p = np.poly1d(z)\n            axes[idx].plot(x.sort_values(), p(x.sort_values()), \"r--\", linewidth=2)\n            \n            # Add correlation\n            r, p_val = pearsonr(x, y)\n            axes[idx].set_title(f'NDVI vs {target}\\nr={r:.3f}, p={p_val:.4f}')\n            axes[idx].set_xlabel('Pre_GSHH_NDVI')\n            axes[idx].set_ylabel(f'{target} (g)')\n            axes[idx].grid(True, alpha=0.3)\n    \n    axes[5].axis('off')\n    plt.tight_layout()\n    plt.show()\n    \n    # Scatter plots: Height vs all targets\n    fig, axes = plt.subplots(2, 3, figsize=(18, 10))\n    axes = axes.flatten()\n    \n    for idx, target in enumerate(EXPECTED_TARGETS):\n        mask = df['Height_Ave_cm'].notna() & df[target].notna()\n        if mask.sum() > 0:\n            axes[idx].scatter(df.loc[mask, 'Height_Ave_cm'], df.loc[mask, target], \n                            alpha=0.4, s=20, color='green')\n            \n            # Add regression line\n            x = df.loc[mask, 'Height_Ave_cm']\n            y = df.loc[mask, target]\n            z = np.polyfit(x, y, 1)\n            p = np.poly1d(z)\n            axes[idx].plot(x.sort_values(), p(x.sort_values()), \"r--\", linewidth=2)\n            \n            # Add correlation\n            r, p_val = pearsonr(x, y)\n            axes[idx].set_title(f'Height vs {target}\\nr={r:.3f}, p={p_val:.4f}')\n            axes[idx].set_xlabel('Height_Ave_cm')\n            axes[idx].set_ylabel(f'{target} (g)')\n            axes[idx].grid(True, alpha=0.3)\n    \n    axes[5].axis('off')\n    plt.tight_layout()\n    plt.show()\n    \n    # NDVI vs Height relationship\n    print(\"\\n📊 NDVI vs Height Relationship:\")\n    fig, ax = plt.subplots(figsize=(10, 6))\n    mask = df['Pre_GSHH_NDVI'].notna() & df['Height_Ave_cm'].notna()\n    \n    if mask.sum() > 0:\n        scatter = ax.scatter(df.loc[mask, 'Height_Ave_cm'], \n                           df.loc[mask, 'Pre_GSHH_NDVI'],\n                           c=df.loc[mask, 'Dry_Total_g'], \n                           cmap='viridis', alpha=0.6, s=50)\n        \n        x = df.loc[mask, 'Height_Ave_cm']\n        y = df.loc[mask, 'Pre_GSHH_NDVI']\n        z = np.polyfit(x, y, 1)\n        p = np.poly1d(z)\n        ax.plot(x.sort_values(), p(x.sort_values()), \"r--\", linewidth=2)\n        \n        r, p_val = pearsonr(x, y)\n        ax.set_title(f'NDVI vs Height (colored by Dry_Total_g)\\nr={r:.3f}, p={p_val:.4f}', \n                    fontsize=14, fontweight='bold')\n        ax.set_xlabel('Height_Ave_cm', fontsize=12)\n        ax.set_ylabel('Pre_GSHH_NDVI', fontsize=12)\n        plt.colorbar(scatter, label='Dry_Total_g', ax=ax)\n        ax.grid(True, alpha=0.3)\n        plt.show()\n\nelse:\n    print(\"⚠️  NDVI and/or Height columns not available\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}