Fitting the population dynamics statistical model
clc;
clearvars -except E D; % clear base workspace except for data struct/tables
if exist('E','var') == 0
E = getfield(load(defaults.files('exported_jpca_matfile')),'E');
end
Define model specifications
We use Wilkinson Notation to define two generalized linear mixed effects models; one simple model that ignores the behavioral covariate and one that accounts for Reach Duration.
Simple Model
We may include some or all of the following terms in our model (select from the drop-down box):
Response
- R2_Skew : This is the average coefficient of determination for planes fit via the skew-symmetric constrained regression. If this value is high, it means that we can accurately linearize the dynamical system that is designed to capture rotatory projections.
Fixed Effects
- GroupID : Categorical fixed main effect with two levels: "Ischemia" or "Intact"
- PostOpDay : Continuous fixed main effect on the range [3, 30].
- GroupID*PostOpDay : Main interaction effect, which is really what we're interested in.
- Intercept : Note that MATLAB automatically adds an intercept term unless it's explicitly removed from the expression.
Random Effects
- R2_Best : This is a covariate that describes how much of the least-squares optimal (unconstrained) regression linearized system we can capture; it is a control that accounts for cases where we simply can't really fit the derivative of the observed trajectories very well.
- Explained_Skew : This describes the fraction of the total data (including data compression from the PCA dimensionality reduction step) explained by the skew-symmetric planar projection. This is the other part that must be accounted for--it is pretty meaningless if we can reconstruct the data well, if the data only lives in a very small subspace of the total observed basis!
Simple Model List (Wilkinson Model Specifications)
- "R2_Skew~GroupID*PostOpDay+(1+R2_Best*Explained_Skew|AnimalID)"
- "R2_Skew~GroupID*PostOpDay+R2_Best+Explained_Skew+(1|AnimalID)"
- "R2_Skew~GroupID*PostOpDay+R2_Best*Explained_Skew+(1|AnimalID)"
- "R2_Skew~GroupID*PostOpDay+(1|AnimalID)+(-1+R2_Best|GroupID)+(-1+Explained_Skew|GroupID)"
simple_model = "R2_Skew~GroupID*PostOpDay+(1+R2_Best*Explained_Skew|AnimalID)";
Complex Model
The complex model includes all the terms of the simple_model but additionally has one more fixed effect in the main interaction:
- Duration : This is the reach duration, which is known to be an important covariate of neural dynamics. It is clear that reach duration on average becomes shorter in rats as they recover from focal ischemia.
Complex Model List (Wilkinson Model Specifications)
- "R2_Skew~GroupID*PostOpDay*Area+(1+R2_Best*Explained_Skew*AverageDuration|GroupID:AnimalID)"
- "R2_Skew~GroupID*PostOpDay*AverageDuration+(1|AnimalID)+(-1+R2_Best*Explained_Skew|PlaneIndex)"
- "R2_Skew~GroupID*PostOpDay+AverageDuration+R2_Best+Explained_Skew+(1|AnimalID)"
- "R2_Skew~GroupID*PostOpDay+AverageDuration+R2_Best*Explained_Skew+(1|AnimalID)"
- "R2_Skew~GroupID*PostOpDay+AverageDuration+(1|AnimalID)+(-1+R2_Best|GroupID)+(-1+Explained_Skew|GroupID)"
- "R2_Skew~AverageDuration*GroupID+(1+PostOpDay|AnimalID)+(-1+R2_Best|GroupID)+(-1+Explained_Skew|GroupID)"
- "R2_Skew~GroupID*PostOpDay*R2_Best*Explained_Skew+(PostOpDay|GroupID:AnimalID)"
complex_model = "R2_Skew~GroupID*PostOpDay+(1+R2_Best*Explained_Skew*AverageDuration|GroupID:AnimalID)";
Specify a list of parameters to use with the fitglme method.
fit_params = {...
'Verbose',0,... % [0, 1, 2] -- how much is reported by optimizer
'FitMethod','ApproximateLaplace', ... {'Laplace', 'REMPL', or 'MPL'}
'DummyVarCoding','effects', ... {'effects','reference', or 'full'}
'UseSequential',false ...
};
Once the desired models have been selected, click Run Fit
fprintf(1,'Current simple model: <strong>"%s"</strong>\n',simple_model);
fprintf(1,'Current complex model: <strong>"%s"</strong>\n',complex_model);
Fit models
Use the fitglme function to recover two generalized linear mixed-effects models, which are then compared in order to determine if it is necessary to include Duration as a model term, since this has a substantial effect on interpreting the model. Here, we do not suppress the output to the Command Window so that we can see what statistics are returned. Note that in order to use a likelihood ratio test for comparison of the two models, we must use the Laplace (we use the approximate Laplace, which approximates the integrand of the likelihood over all latent variables in the model; the approximation is generally good unless the observations are sparse for a given level) fitting method, as MPL or REMPL methods do not allow models to be compared since in MPL the maximized log-likelihood of the pseudodata from the final pseudo likelihood iteration is used in the likelihood ratio test.
Fit simple models
glme_simple = fitglme(E,simple_model,fit_params{:})
To interpret this result:
- (Intercept) : There is a statistically significant effect, and this is easily the largest portion of the model in itself. The easiest way to think about this is it's the model average output; all other fixed effects should sum to zero. This value plus the random intercept value should sum to the average observed response (the average ability to reconstruct a given skew-symmetric plane).
Note that the (Fixed-Effect) Intercept is different from the average of the Response variable:
disp(mean(E.R2_Skew));
This difference is accounted for by the Random-Effect Intercept terms, which are grouped by AnimalID.
disp(glme_simple.Coefficients.Estimate(1)-mean(E.R2_Skew));
disp(glme_simple.Rsquared)
Fit complex models
glme_complex = fitglme(E,complex_model,fit_params{:})
disp(glme_complex.Rsquared)
[~,~,stats] = randomEffects(glme_complex) %#ok<ASGLU>
Compare simple and complex models
if glme_simple.LogLikelihood > glme_complex.LogLikelihood
glme = glme_simple;
mdl = simple_model;
else
glme = glme_complex;
mdl = complex_model;
end
fprintf(1,'Selected Model: <strong>%s</strong>\n',mdl);
If the pvalue is large, we should retain glme_simple; otherwise, the additional model complexity adds enough meaningful information that we should keep glme_complex.
Check the residuals
We can plot the residuals using a method of the glme object that we've recovered:
analyze.stat.panelized_residuals(glme);
The fitted residuals plot looks reasonable. It's balanced on both the X- and Y-axes. Our residuals are not going to pass a test of normality probably, but they look good enough for government work.
Plot residuals grouped by AnimalID
We should verify that there isn't an outlier Animal (as opposed to single-trial) that is biasing the regression.
[groupings,TID] = findgroups(E(:,{'AnimalID'}));
analyze.stat.plot_glme_residuals(glme,E,groupings,TID,'yVar','R2_Skew','scatterAlpha',0.45,'scatterSize',15);
It looks like RC-02 hand RC-26 are really mucking up our dataset. Both have very few observations, and since our model depends on estimating regressions from, among other terms, data grouped by each Animal, we should probably just exclude those two from our model. We should re-fit the model without these two categories.
Re-fit the model, excluding RC-02 and RC-26 observations
excludeRat = "RC-02";
Once we have re-fit, we should double-check our residuals once again and then we can begin interpreting our model.
exclusions = ismember(string(E.AnimalID),excludeRat);
E.Properties.UserData.ModelExclusions = exclusions;
E.Properties.UserData.Excluded.glme_outlier_residuals = sum(exclusions);
E.Properties.UserData.FinalGLMEparams = fit_params{:};
glme2 = fitglme(E,mdl,'Exclude',exclusions,fit_params{:})
disp(glme2.Rsquared)
[~,~,stats] = randomEffects(glme2) %#ok<ASGLU>
[~,~,stats_fixed] = fixedEffects(glme2)
analyze.stat.panelized_residuals(glme2);
[groupings,TID] = findgroups(E(~exclusions,{'AnimalID'}));
analyze.stat.plot_glme_residuals(glme2,E(~exclusions,:),groupings,TID,'yVar','R2_Skew','scatterAlpha',0.33);
The fits have improved after removing the two outlier rats (each of which only had very few observations. Clearly, the number of observations for RC-14 isn't great, but after exclusions only 2 animals remain in the Intact group to begin with so we will leave it in.
Verify that PostOpDay also doesn't have some funny residuals
[groupings,TID] = findgroups(E(~exclusions,{'PostOpDay'}));
analyze.stat.plot_glme_residuals(glme2,E(~exclusions,:),groupings,TID,'yVar','R2_Skew','scatterAlpha',0.33);
The residuals look pretty close to the reference line for the remaining post-operative days; therefore, we are done verifying the model. We could also try to re-run it by excluding some of the latest days, restricting the analysis to the first 20 days, which we have more data points for. Excluding the first week of data, when the rat in general is very inactive, may also be a good idea.
Define a second set of exclusions, and re-run the model
minPostOpDay = 8;
maxPostOpDay = 22;
combineExclusions = false;
Click Run once the range is set.
if combineExclusions
exclusions2 = exclusions | (E.PostOpDay < minPostOpDay) | (E.PostOpDay > maxPostOpDay);
else
exclusions2 = (E.PostOpDay < minPostOpDay) | (E.PostOpDay > maxPostOpDay); %#ok<*UNRCH>
end
E.Properties.UserData.ModelExclusions2 = exclusions2;
E.Properties.UserData.Excluded.glme_outlier_residuals2 = sum(exclusions2);
glme3 = fitglme(E,mdl,'Exclude',exclusions2,fit_params{:})
disp(glme3.Rsquared)
[~,~,stats] = randomEffects(glme3)
analyze.stat.panelized_residuals(glme3);
[groupings,TID] = findgroups(E(~exclusions2,{'AnimalID'}));
analyze.stat.plot_glme_residuals(glme3,E(~exclusions2,:),groupings,TID,'yVar','R2_Skew','scatterAlpha',0.33);
[groupings,TID] = findgroups(E(~exclusions2,{'PostOpDay'}));
analyze.stat.plot_glme_residuals(glme3,E(~exclusions2,:),groupings,TID,'yVar','R2_Skew','scatterAlpha',0.33);
Interpreting the jPCA planes
This model tells us that there is something about the jPCA planes that separates the two groups, particularly during the early post-operative phase; however, we can use the projection matrices to track differences back to the original spike rates so that we can get a better "physical" interpretation of what they mean.
1) Find example Rat from Ischemia group with differences during Early and Late Post-Op recordings
RC-04 or RC-30 show low ability to reconstruct rotatory dimensions during early post-op recordings, but stronger rotations during late post-op recordings.
example_animal = "RC-04";
e = E(E.AnimalID==example_animal,:);
analyze.jPCA.scatter_R2skew(e);
2) Recover dynamic channel weightings as PCs evolve over time
The jPCA method allows us to recover M_skew, such that dX = M_skew * X; here, X is not the observed spike rates, but rather, the scores of principal components of the spike rates in terms of all their possible simultaneous combinations of "rate states." Indexing makes it easier for me to think about this: take M_skew(1,:) and M_skew(2,:) -- these should correspond to the primary jPC plane, assuming the vectors are sorted in ordered pairs. Each column of M_skew(1,:) corresponds to a weighting for each of the (probably 12 total) PCs (rows). In turn, columns of each row give us a time-series: how the scores vary across all the observed samples from each trial. Therefore, to get a single-channel weighting through time for a single jPC vector, we recover the original rates using the projected principal components: W = X*M_skew*(coeff') + repmat(mu,size(X,1),1) where coeff is the matrix of principal component vectors, and mu is the (row) vector of individual channel rates.
% if exist('D','var')==0 % Load in `multi_jPCA` exported data table
% D = getfield(load(defaults.files('multi_jpca_default_matfile')),'D');
% end
% d = D(D.AnimalID==example_animal,:);
% [groupings,tid] = findgroups(d(:,{'AnimalID','Alignment','PostOpDay'}));
% P = splitapply(@(P,S)analyze.jPCA.recover_channel_weights(P{:},S{:}),...
% d.Projection,d.Summary,groupings);
3) Approach as normal with "rosette" or "movie" strategy for visualization
% batch.exportMovies(P);