%% MATLAB basics and regression
% Table of Contents
% 1. Why MATLAB?
% 2. Getting started
% 3. Data types and operations
%     Numeric types
%     Characters and Strings
%     Tables
%     Cell Arrays
%     Structures
% 4. Regression
% 5. Conclusion
%%
%% 1. Why MATLAB?
% MATLAB is a proprietary programming language and computing enviorment developed by MathWorks. It is specifically desinged for engineers and scientists who work with data matrices of large sizes. Some of its advantages over other options (e.g., R, Python) include:
% having a graphical user interface (GUI)
% (relatively) easy-to-understand data types
% comprehensive tried-and-tested built-in functions
% various toolboxes for brain research and beyond
% Here are some famous/useful fMRI packages that you might use:
% SPM
% one of the most popular packages for preprocessing and analyzing neuroimaging data, besides FSL and AFNI
% CONN
% MATLAB- and SPM-based interactive toolbox for analyzing brain connectivity
% EEGLAB
% interactive toolbox for processing EEG, MEG, and other electrophysiological data
% FieldTrip
% non-interactive toolbox for processing MEG, EEG, and iEEG data
% PsychToolbox
% toolbox for controling visual and auditory stimulus presentation
%
% And here are some previous DIBS Methods Meeting posts that used MATLAB to some extent:
% MATLAB basics
% Multivariate pattern analysis
% EEG Preprocessing using EEGLAB
% Decomposing Fourier transforms
%% 2. Getting started
% Let's get started! If you haven't installed MATLAB or want to familiarize yourself with the GUI, please refer to my earlier post on MATLAB basics OR try out the online version here: https://matlab.mathworks.com/.
%
%% 3. Data types and operations
% This was covered in a bit more depths in my earlier post, but I thought I would do a quick review of some important things here because they are going to come up later as we talk about regression.
%% Numeric types

% double	Double-precision arrays (64-bit) <-- MATLAB default
% single	Single-precision arrays (32-bit)
% int8	    8-bit signed integer arrays
% uint8	    8-bit unsigned integer arrays
% ...
% for example:
class(1) % use the function class to check data types
% use space or comma to separate columns
% use semicolon to separate rows
% use square brackets [ ] to put things together
a = [1    2,    NaN;
    .5    Inf   nan] % this is a 2 x 3 double array
% note the two special numbers: Inf and NaN (not a number)
class(a) % the data type is still 'double'

%% Characters and Strings

% use single quotes for characters
char1 = 'a';
char2 = 'abc';
% use double quotes for strings
str1 = "a";
str2 = "abc";
% check data types
class(char1);
class(str1);
% difference
length(char2)
length(str2)
% if we want to concatnate different things
[char1 char2] % correctly concatenated character array
str1 + str2 % correctly concatenated strings
char1 + char2 % data type converted to numeric for addition
[char1 str2] % character coerced to string

%% Tables

T = table;
T.col1 = [1; 2];
T.col2 = ["hello"; "world"];
T
% dimensions must match
T.col3 = repmat(3, height(T), 1);
try T.col4 = 4; catch err; disp(err); end
% data type must be consistent within each row
T.col2(2) = 1;
T % see how 1 is coerced into "1"
class(T.col2)

%% Cell Arrays

% use curly brackets to put together things of different sizes and types
C = {42, "abcd"; table(nan), [1 2 3]; Inf, {}}
C(2) % paratheses indexing retrieves the cell 
C{2} % curly braces indexing retrieves the cell *content*

% Arrays are effcient, but make sure you access the in the correct order.

% double arrays are indexed with parentheses
a = [11 12 13; 21 23 24]; % 2 x 3 double array
a
a(1, 2) % row 1, column 2
a(1:end, 2) % all rows, column 2
a(2, :) % row 2, all columns
% here's something that may seem odd
a(4) % 4th entry
a(:) % all entries

% In Python and R, array data is stored row-wise. In contrast, MATLAB arrays are stored column-wise, even though they are easily defined row-wise.
% ┑(￣Д ￣)┍
%% Structures

% group data using fields
% each field can be of any data type, including structures
S = struct;
S.field1_struct = struct;
S.field2_cell = C;
S.field3_table = T;
S.field4_char = char1;
S

%
%% 4. Regression
% Let's talk about regression! Through the official Statistics and Machine Learning Toolbox, we have access to several built-in MATLAB functions for regression.
% First, let's load some data. We are going to use an open dataset on Kaggle on life expectancy. The original data came from the World Health Organization (WHO), who has been keeping track of the life expectancy and many other health factors of all countries. The final dataset consists of 20 predictor variables and 2938 rows, containing information for 193 countries between 2000 and 2015.

data_table = readtable("2024-03-01-matlab-regression-Life-Expectancy.csv", ...
    "VariableNamingRule", "preserve"); % preserve white space in column names for readability
head(data_table);

% Let's validate the information in the description above.

fprintf( ...
    "Number of rows = %d \n" + ...
    "Number of years = %d \n" + ...
    "Number of countries = %d \n", ...
    height(data_table), ...
    length(unique(data_table.Year)), ...
    length(unique(data_table.Country)));

% For simplicity, we are going to focus on the most recent complete sample (year 2014) and on the following variables:
% Life expectancy (in years)
% Status: "Developed", "Developing"
% Total expenditure: General government expenditure on health as a percentage of total government expenditure (%)

data_2014 = data_table(data_table.Year==2014, ["Country" "Status" "Total expenditure" "Life expectancy"]);
head(data_2014);

% Before actually fitting a linear regression model, let's plot the data.

close all
figure
gscatter( ...
    data_2014.("Total expenditure"), ... x-axis
    data_2014.("Life expectancy"), ... y-axis
    data_2014.Status ... color
    );
title("Life Expectancy against Healthcare Expenditure in 2014")

% We make 3 observations:
% There are much fewer developed countries (orange) than developing countries (blue).
% Developed countries tend to have higher life expectancy than developing countries.
% Life expectance MAYBE is positively correlated with heathcare expenditure for developing countries, but less so for developed countries.
% Now let's fit some linear regression models!
% We're going to use the `fitlm` function in MATLAB. This function provides very detailed outputs.

data_2014.Status = categorical(data_2014.Status, ["Developing" "Developed"]); %%% note the order
data_2014.LE = data_2014.("Life expectancy");
data_2014.HE = data_2014.("Total expenditure");
m1 = fitlm(data_2014, "LE ~ HE * Status");
% anova(m1, "component", 3) %%% Type III anova
m1 %%% regression coefficients and stats

% Note that it's critically important to know how to correctly interpret these results, e.g., what is the "Intercept" and whether "HE" is a simple effect or a main effect. See more in Kevin's post on Interpreting Regression Coefficients. Briefly, the order of the categorical variable AND whether the continuous variable is mean-centered matters. Let's see:

data_2014.Status_rev = categorical(data_2014.Status, ["Developed" "Developing"]);
data_2014.HE_mc = data_2014.HE - mean(data_2014.HE, "omitmissing");
m2 = fitlm(data_2014, "LE ~ HE_mc * Status");
m3 = fitlm(data_2014, "LE ~ HE_mc * Status_rev");
m1.Coefficients
m2.Coefficients
m3.Coefficients

% Let's quickly plot the effects of HE (with 95% confidence intervals) on LE, separately for developed and developing countries. This is essentially doing `plot(ggemmeans(m, ~ HE + Status))` in R.

plotSlice(m2); %%% alternative syntax: `m2.plotSlice;` 

%
% However, If we need to fit a large number of models and don't really need detailed statistics from each model, we can use the `regress` function to trade comprehensiveness for speed.

% define outcome variable
y = data_2014.LE;
% define predictor variables
X = nan(length(y), 4); %%% initialize 3 columns
X(:, 1) = 1; %%% constant term is NOT automatically put into the model!
X(:, 2) = double(data_2014.Status=="Developing"); %%% manually code the categorical variable
X(:, 3) = data_2014.("Total expenditure") - mean(data_2014.("Total expenditure"), "omitmissing");
X(:, 4) = X(:, 2) .* X(:, 3); %%% interaction term
[b, bint, r, rint, stats] = regress(y, X);
b
m3.Coefficients
all(m3.Coefficients.Estimate == b)

% Ta-da! As we can see the regression coefficients we obtained using `regress` are exactly the same as what we had from `fitlm` earlier, though we no longer have the nice-looking table filled with stats. However, the (very crude) test below shows that `regress` is indeed a lot faster than `fitlm`.
%

tic
for i=1:500
    y = data_2014.LE;
    % define predictor variables
    X = nan(length(y), 4); %%% initialize 3 columns
    X(:, 1) = 1; %%% constant term is NOT automatically put into the model!
    X(:, 2) = double(data_2014.Status=="Developing"); %%% manually code the categorical variable
    X(:, 3) = data_2014.("Total expenditure") - mean(data_2014.("Total expenditure"), "omitmissing");
    X(:, 4) = X(:, 2) .* X(:, 3); %%% interaction term
    [b, bint, r, rint, stats] = regress(y, X);
end
toc
tic
for i=1:500
    m = fitlm(data_2014, "LE ~ HE_mc * Status_rev");
end
toc

%% 5. Conclusion
% We've made it to the end of MATLAB basics and regression (yay!), though this is just the tip of the iceberge of regression models out in the wild world, such as regularized regression, nonlinear regression, and mixed-effects regression. Feel free to check out how they are can be implemented in MATLAB here.
%
% Enjoying MATLAB-ing!