% experiment_results.m
% Loads measurement data from measurements.csv, splits by condition,
% fits a linear regression of treatment-A measurements against control
% measurements, and plots the result with the fitted line.
%
% Companion file for the MathJet tutorial:
%   "Bringing your MATLAB scripts into MathJet"
%   https://mathjet.com/docs/tutorials/matlab-scripts/

% --- 1. Load the data ----------------------------------------------------
data = readtable('measurements.csv');

control   = data(strcmp(data.condition, 'control'),     :);
treatment = data(strcmp(data.condition, 'treatment_a'), :);

% Pair the first N rows from each group (data shuffle order is independent).
n = min(height(control), height(treatment));
x = control.measurement_value(1:n);
y = treatment.measurement_value(1:n);

% --- 2. Fit a linear regression ------------------------------------------
p     = polyfit(x, y, 1);   % p(1) = slope, p(2) = intercept
y_fit = polyval(p, x);

ss_res = sum((y - y_fit).^2);
ss_tot = sum((y - mean(y)).^2);
r2     = 1 - ss_res / ss_tot;

% --- 3. Plot scatter + fit -----------------------------------------------
figure;
scatter(x, y, 36, 'filled');
hold on;
plot(x, y_fit, 'r-', 'LineWidth', 2);
hold off;

xlabel('Control measurement');
ylabel('Treatment A measurement');
title(sprintf('Linear fit: y = %.3fx + %.3f   (R^2 = %.3f)', p(1), p(2), r2));
grid on;
legend({'Observed', 'Linear fit'}, 'Location', 'best');

% --- 4. Print summary ----------------------------------------------------
fprintf('\nLinear regression: treatment_a ~ control\n');
fprintf('  Slope:     %.4f\n', p(1));
fprintf('  Intercept: %.4f\n', p(2));
fprintf('  R-squared: %.4f\n', r2);
fprintf('  N pairs:   %d\n\n', n);