Gallery

// Which schools do better than their students' prior scores predict?
// Fit one regression on all students, then average each school's residuals.
$schools = DataFrame::fromCsv('schools.csv');
$fit = phplyr($schools)->lm('score ~ prior_score + study_hours');
[$b0, $b_prior, $b_hours] = $fit['coefficients'];

echo phplyr($schools)
    ->mutate(expected: fn ($r) => $b0 + $b_prior * $r->prior_score + $b_hours * $r->study_hours)
    ->mutate(residual: fn ($r) => $r->score - $r->expected)
    ->groupBy('school')
    ->summarize(
        students: n(),
        avg_score: mean('score'),
        avg_residual: mean('residual'),
        share_above: fn ($g) => round(\count(array_filter($g->column('residual'), fn ($e) => $e > 0)) / $g->nrow(), 2),
    )
    ->mutate(rank: min_rank(desc('avg_residual')), avg_score: fn ($r) => round($r->avg_score, 1), avg_residual: fn ($r) => round($r->avg_residual, 2))
    ->arrange('rank');

// A random-intercept model: students nested in schools. The fixed effects
// are shared; each school gets its own intercept shift.
$schools = DataFrame::fromCsv('schools.csv');
$fit = phplyr($schools)->lmer('score ~ study_hours + prior_score + (1 | school)');

foreach ($fit['columns'] as $i => $term) {
    printf("%-12s %8.4f  (se %.4f)\n", $term, $fit['coefficients'][$i], $fit['std_errors'][$i]);
}
$schoolSd = $fit['theta'][0] * $fit['sigma'];
printf("school sd %.3f, residual sd %.3f, ICC %.3f\n\n", $schoolSd, $fit['sigma'], $schoolSd ** 2 / ($schoolSd ** 2 + $fit['sigma'] ** 2));

// Compare each school's raw mean with its estimated intercept shift.
$random = $fit['random'][0];
$effects = array_combine($random['group_labels'], array_column($random['effects'], 0));
echo phplyr($schools)
    ->groupBy('school')
    ->summarize(students: n(), raw_mean: mean('score'))
    ->mutate(intercept_shift: fn ($r) => round($effects[$r->school], 2), raw_mean: fn ($r) => round($r->raw_mean, 1))
    ->arrange(desc('intercept_shift'))
    ->head(5);

// Price per square meter by number of bedrooms, with each flat's
// position relative to its bedroom group.
$housing = phplyr(DataFrame::fromCsv('housing.csv'))
    ->mutate(per_m2: fn ($r) => round($r->monthly_rent / $r->area_m2, 2))
    ->get();

echo phplyr($housing)
    ->groupBy('bedrooms')
    ->summarize(across(['area_m2', 'monthly_rent', 'per_m2'], ['mean' => mean(...), 'max' => max(...)]), flats: n())
    ->mutate(across(where('double'), fn ($v) => round($v, 1)));

echo "\n", phplyr($housing)
    ->mutate(vs_group: fn ($r, $g) => round($r->per_m2 - mean('per_m2')($g), 2), _by: 'bedrooms')
    ->arrange('bedrooms', desc('vs_group'));

// Pass rates by study-hour quartile, then a logistic regression.
$learning = DataFrame::fromCsv('learning.csv');

echo phplyr($learning)
    ->mutate(hours_quartile: ntile('study_hours', 4))
    ->groupBy('hours_quartile')
    ->summarize(
        students: n(),
        min_hours: min('study_hours'),
        max_hours: max('study_hours'),
        pass_rate: mean('passed'),
    );

$fit = phplyr($learning)->glm('passed ~ study_hours + practice_score', family: 'binomial');
printf("\nOdds ratio per extra study hour: %.3f\n", exp($fit['coefficients'][1]));

// Calls by staffing level and visitor volume, as a cross-tabulation.
echo phplyr(DataFrame::fromCsv('clinic.csv'))
    ->mutate(
        staffing: case_when([fn ($r) => $r->staff < 4, 'low'], [fn ($r) => $r->staff < 6, 'mid'], default: 'high'),
        volume: if_else(fn ($r) => $r->visitors >= 250, 'busy', 'quiet'),
    )
    ->groupBy('staffing', 'volume')
    ->summarize(mean_calls: mean('calls'), _groups: 'drop')
    ->mutate(mean_calls: fn ($r) => round($r->mean_calls, 2))
    ->pivotWider(namesFrom: 'volume', valuesFrom: 'mean_calls', namesSort: true)
    ->arrange(fn ($r) => array_search($r->staffing, ['low', 'mid', 'high']));