Statistical Models — MATLAB lab work

Multivariate statistics coursework (UPF, Statistical Models, 3rd year, autumn 2023), written in MATLAB. Three labs, each a self-contained hypothesis-testing exercise:

Almost everything is implemented from first principles rather than by calling MATLAB's built-in manova1 / pca / bootstrp: Mahalanobis distances, SSB/SSW decomposition, Wilks' Lambda and its F transformation, the bootstrap loops, the covariance eigen-decomposition, the ISOMAP geodesic embedding and the Kernel PCA Gram matrix are all written out explicitly. That was the point of the course — the built-ins appear only for distribution quantiles (norminv, chi2inv, finv, tinv, prctile) and for plotting.

Taken from the students' own commentary in the scripts:

• Normality is checked before, not assumed. Every test is preceded by univariate Q-Q correlation checks against the critical-point table, then chi-square plots of Mahalanobis distances for every pair of variables and for all variables jointly. A chi-square plot is read through the slope of a fitted line: a slope far from 1 is evidence against normality.
• When normality fails, fix the data rather than the test. In Lab 2 exercise 1, outlier removal was tried and rejected (the best corrected Q-Q correlation reachable was 0.9279 against a critical value of 0.9591), and a log transformation was used instead — which did bring Group A into a range where normality could not be rejected.
• Equal covariance is verified before MANOVA, by checking that the ratio of the squared diagonal elements of the two group covariance matrices lies inside [1/4, 4].
• Bootstrap is used as a cross-check, not a replacement. In Labs 4 and 6, every parametric conclusion is repeated with a bootstrap percentile interval, and the two are compared explicitly.
• Multiple comparisons are Bonferroni-corrected (Lab 4 uses alpha_corrected = alpha/3).

src/
  lab2_ex1_gait_manova.m            MANOVA on gait-analysis data, 4 patient groups x 3 variables
  lab2_ex2_repeated_measures.m      Repeated-measures ANOVA on 5 lab grades for 57 students
  lab2_ex3_fish_manova.m            MANOVA on a 3-method x 4-attribute fish tasting panel
  lab4_bootstrap.m                  Lab 4, converted from the .mlx live script
  lab6_dimensionality_reduction.m   Lab 6, converted from the .mlx live script
livescripts/
  lab4_bootstrap.mlx                Original MATLAB live scripts (the submitted artifacts)
  lab6_dimensionality_reduction.mlx
data/                               The five course-provided .xlsx datasets (all anonymous)
docs/
  RESULTS.md                        Every numeric result recovered from the stored outputs
  SM2023_Lab02.pdf                  Course-provided lab statements
  SM2023_Lab04.pdf
  SM2023_Lab06.pdf

None of these contain names or identifiers — they are anonymous course datasets.

Requires MATLAB with the Statistics and Machine Learning Toolbox (norminv, chi2inv, finv, tinv, prctile, skewness, randsample, anova1). No add-ons beyond that.

cd data          % the scripts call xlsread with bare filenames, so run from data/
run('../src/lab2_ex1_gait_manova.m')
run('../src/lab2_ex2_repeated_measures.m')
run('../src/lab2_ex3_fish_manova.m')
run('../src/lab4_bootstrap.m')
run('../src/lab6_dimensionality_reduction.m')

Lab 4 takes 2-3 minutes end to end (the authors' own note in the live script) because of the number of bootstrap replicates and plots. Lab 6 runs 1000 column-wise permutations of a 285x137 matrix and is also slow.

To open the originals as live scripts instead, open livescripts/*.mlx in MATLAB — those keep the rendered equations, the formatted prose and the stored figures, none of which survive the conversion to .m.

All numbers below were read out of the stored outputs embedded in the original .mlx files (matlab/output.xml) or out of comments the authors wrote after running the code. Nothing here was re-run or recomputed.

Conclusions as written by the authors: none of the three bootstrap p-values approach the corrected alpha, so multivariate normality of Group A cannot be rejected; all three observed statistics fall inside their percentile intervals, which agrees. For Problem 2, T² = 32.2547 exceeds the critical value and the p-value is far below alpha, so the means of Groups A and D differ. For Problem 3 the Student-t interval and the bootstrap percentile interval for skewness disagree — the authors note the t interval assumes normality, which does not hold here, so it is not accurate.

Conclusions as written by the authors: the boy/girl mean difference is detectable in the 41-dimensional PCA space and in the Kernel PCA space, but not in the 3-dimensional ISOMAP embedding — 3 dimensions simply do not carry enough of the variance. Kernel PCA retains essentially the same variance as linear PCA (0.57083 vs 0.57133) and reaches the same conclusion, which the authors read as evidence that the relevant structure in this dataset is linear and KPCA buys nothing here.

The permutation-test component count (41) and the bootstrap intervals will shift slightly between runs — the scripts do not set a random seed. The authors flag this themselves ("41 may vary slightly between executions due to the permutations").

Lab 2 was submitted as plain .m files, which do not store their outputs, so the only recorded numbers are the ones the authors wrote into their comments:

• Best Q-Q correlation reachable for Group D variable 2 after outlier removal: 0.9279, against a critical value of 0.9591 → normality still rejected.
• Group A slopes were "very far from 1" before a log transform and "almost 0.8" after it (the authors' words; the exact per-plot slopes were not recorded).
• Lab 2 ex. 2: at 99% confidence the variable pairs whose difference intervals exclude 0 are 1-2, 1-3, 1-5 and 3-4.
• Lab 2 ex. 3: the multivariate post-hoc test rejects equality for pairs A-C and B-C, but not A-B.

• Lab 4 and Lab 6 .m files are machine-converted from the .mlx live scripts by extracting matlab/document.xml and emitting code paragraphs verbatim with prose as comments. The code is byte-for-byte the original; what is lost is: rendered LaTeX equations (Lab 4 Problem 2b's Hotelling's T² formula was an embedded image and is simply absent from the .m), inline figures, and text formatting. The original .mlx files are shipped in livescripts/ for that reason.
• An operator-precedence bug was fixed. The original code wrote f_transform = ((1-(wilks_lambda)^1/2) / 6) / ((wilks_lambda^1/2)/df2). In MATLAB ^ binds tighter than /, so wilks_lambda^1/2 evaluates as (wilks_lambda^1)/2 — half of Wilks' Lambda — rather than the square root the Wilks-to-F transformation requires. It now reads wilks_lambda^(1/2) in both lab2_ex1_gait_manova.m and lab2_ex3_fish_manova.m. Consequence: the F-statistics those two lines produce are now different from the values stored in docs/RESULTS.md, which were recorded with the original code. The fix is correct; the stored outputs are the historical record. They are not expected to agree.
• The PDFs in docs/ are course-provided lab statements, not student work.

• José Mª Pérez Clar
• Guillem Sans

Universitat Pompeu Fabra, Statistical Models, 3rd year, Trimester 1, 2023.