Behavior Endpoint: Trial Durations

Live script to plot different filters of total trial segment duration by day and outcome.

clearvars -except T S
if exist('T','var')==0
   T = getfield(load(defaults.files('rate_table_default_matfile'),'T'),'T');
end

Plot Total Duration

"Duration" is the total number of seconds between Reach and Complete.

T.Duration = T.Complete - T.Reach;

T.Properties.VariableUnits{find(strcmp(T.Properties.VariableNames,'Duration'),1,'first')} = 'sec';

analyze.behavior.durations(T);

Dashed lines are individual animal regressions (for all subsequent plots as well). Black lines are group (subplot) median regressions. Regression lines are fit using a logistic link function; see analyze.stat for more details. In general, pellet retrievals became faster in Ischemic rats (left) over the course of the first 3-4 weeks after the injury.

analyze.behavior.durations(T,'Outcome','Successful');

To see whether this is due to an improved Success Rate, we broke down the trial durations by Outcome. Note that we only included trials in which both a Reach and a Grasp were present. We also required that the pellet was present (sometimes rats get frustrated and flail out of the box even when no reward pellet is there to retrieve). The trend holds true for Successful results: Ischemic rats complete Successful trials more quickly as recovery proceeds.

analyze.behavior.durations(T,'Outcome','Unsuccessful');

Interestingly, the reduction in trial duration by post-op day does not seem to be strongly influenced by Outcome. This suggests that, in general, Ischemic rats perform the motor behavior more quickly as recovery progresses.

Is there a segment of the retrieval that biases the speed increases?

If Ischemic rats retrieve pellets faster as they recover from injury, is there a particular segment of the behavior this can be attributed to? For example, is the rat extending its paw and finding the pellet more quickly (i.e. "targeting" improvements) or is the rat sensing that it has come into contact with the pellet more quickly and then able to return its paw to the box accordingly with less latency (i.e. "sensing" improvements).

Plot Extension Duration

"Extension" is the phase between Reach and Grasp.

T.Extension_Duration = T.Grasp - T.Reach;

T.Properties.VariableUnits{find(strcmp(T.Properties.VariableNames,'Extension_Duration'),1,'first')} = 'sec';

analyze.behavior.show(T,'Extension_Duration');

analyze.behavior.show(T,'Extension_Duration','Outcome','Successful');

analyze.behavior.show(T,'Extension_Duration','Outcome','Unsuccessful');

In each case, the Extension duration decreases as a function of post-op day, as was observed with full-trial duration.

Plot Retraction Duration

"Retraction" is the phase between Grasp and Complete.

T.Retraction_Duration = T.Complete - T.Grasp;

T.Properties.VariableUnits{find(strcmp(T.Properties.VariableNames,'Retraction_Duration'),1,'first')} = 'sec';

analyze.behavior.show(T,'Retraction_Duration');

analyze.behavior.show(T,'Retraction_Duration','Outcome','Successful');

analyze.behavior.show(T,'Retraction_Duration','Outcome','Unsuccessful');

In each case, the Retraction duration decreases as a function of post-op day, as was observed with full-trial duration. Beyond what is evident from the scatter trends, to get a more nuanced idea of what is happening we should use a statistical model to quantify these results.

Trial/Segment Duration Generalized Linear Models

First, filter to ensure our table has unique trials only (this is done in the analyze.behavior functions already).

[~,uniqueTrials] = unique(T.Trial_ID);
U = T(uniqueTrials,:);

Initialize a data struct to keep all the generalized linear models grouped together

glm = struct('duration',struct,'outcome',struct);
glm.duration.data = U; % Associates data that was used with the struct as well

Next, use fitglme to estimate a Generalized Linear Mixed-Effects model for each of the responses (Total, Extension, and Retraction durations). We use Wilkinson Notation, where a three-way interaction defined as A*B*C includes the main fixed effects of each as well as the pairwise-interaction effects. Random effects are denoted as (1|R), where R is the grouping variable.

glm.duration.Duration = fitglme(U,"Duration~Group*Outcome*PostOpDay+(1|AnimalID)",'FitMethod',"Laplace");
disp(glm.duration.Duration);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations            4290
    Fixed effects coefficients           8
    Random effects coefficients         10
    Covariance parameters                2
    Distribution                    Normal
    Link                            Identity
    FitMethod                       Laplace

Formula:
    Duration ~ 1 + Group*PostOpDay + Group*Outcome + PostOpDay*Outcome + Group:PostOpDay:Outcome + (1 | AnimalID)

Model fit statistics:
    AIC        BIC       LogLikelihood    Deviance
    -726.34    -662.7    373.17           -746.34 

Fixed effects coefficients (95% CIs):
    Name                                                   Estimate     SE            tStat       DF      pValue        Lower         Upper     
    {'(Intercept)'                                }          0.86176      0.046875      18.384    4282    1.0328e-72       0.76986       0.95366
    {'Group_Intact'                               }        0.0012898      0.076416    0.016879    4282       0.98653      -0.14852        0.1511
    {'PostOpDay'                                  }        -0.012765    0.00087407     -14.604    4282    3.5502e-47     -0.014478     -0.011051
    {'Outcome_Unsuccessful'                       }        -0.054048      0.021857     -2.4728    4282      0.013444     -0.096899     -0.011197
    {'Group_Intact:PostOpDay'                     }         0.020976      0.002032      10.323    4282    1.0804e-24      0.016992      0.024959
    {'Group_Intact:Outcome_Unsuccessful'          }        -0.081696      0.035727     -2.2867    4282      0.022263      -0.15174     -0.011653
    {'PostOpDay:Outcome_Unsuccessful'             }        0.0031369     0.0012282      2.5539    4282      0.010685    0.00072887     0.0055449
    {'Group_Intact:PostOpDay:Outcome_Unsuccessful'}        -0.012732     0.0023525      -5.412    4282    6.5722e-08     -0.017344    -0.0081197

Random effects covariance parameters:
Group: AnimalID (10 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.1095  

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        0.22064 

Result: Total Duration

There is a significant main fixed effect of intercept (p << 0.001), indicating that an average trial duration can be expected to be approximately 860-ms.
There is no significant main fixed effect of Group (p = 0.99).
There is a significant main fixed effect of Post-Op Day (p << 0.001), which indicates that there is an overall tendency to decrease reach duration by approximately 12-ms per day.
There is a significant main fixed effect of Outcome (p = 0.013), which indicates that Unsuccessful trials are on average about 54-ms faster than Successful ones.
There is a significant fixed interaction of Group*Post-Op Day (p << 0.001), which indicates that Intact rats tend to have a slight increase in trial duration (about 21-ms per day) by comparison to Ischemic rats.
There is a significant fixed interaction of Group*Outcome (p = 0.022), which indicates that Intact rats tend to perform Unsuccessful retrievals about 80-ms faster on average.
There is a significant fixed interaction effect of PostOpDay*Outcome (p = 0.011), which indicates that Unsuccessful retrievals take about 3-ms longer per day on average.
There is a significant three-way fixed interaction of Group*PostOpDay*Outcome (p << 0.001), which indicates that on average, Intact rats decrease duration of unsuccessful pellet retrievals as time progresses.

disp(glm.duration.Duration.Rsquared);
    Ordinary: 0.2628
    Adjusted: 0.2616
glm.duration.Extension_Duration = fitglme(U,"Extension_Duration~Group*Outcome*PostOpDay+(1|AnimalID)",'FitMethod',"Laplace");
disp(glm.duration.Extension_Duration);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations            4352
    Fixed effects coefficients           8
    Random effects coefficients         10
    Covariance parameters                2
    Distribution                    Normal
    Link                            Identity
    FitMethod                       Laplace

Formula:
    Extension_Duration ~ 1 + Group*PostOpDay + Group*Outcome + PostOpDay*Outcome + Group:PostOpDay:Outcome + (1 | AnimalID)

Model fit statistics:
    AIC        BIC        LogLikelihood    Deviance
    -4427.6    -4363.8    2223.8           -4447.6 

Fixed effects coefficients (95% CIs):
    Name                                                   Estimate       SE            tStat       DF      pValue        Lower         Upper     
    {'(Intercept)'                                }            0.38235       0.03179      12.027    4344    8.4281e-33       0.32003       0.44468
    {'Group_Intact'                               }           0.069223      0.051713      1.3386    4344       0.18077     -0.032161       0.17061
    {'PostOpDay'                                  }         -0.0043886      0.000572     -7.6724    4344    2.0709e-14    -0.0055101    -0.0032672
    {'Outcome_Unsuccessful'                       }          0.0013424      0.014303    0.093849    4344       0.92523       -0.0267      0.029384
    {'Group_Intact:PostOpDay'                     }          0.0054365     0.0013297      4.0885    4344    4.4206e-05     0.0028296     0.0080434
    {'Group_Intact:Outcome_Unsuccessful'          }           -0.04489      0.023375     -1.9204    4344      0.054868     -0.090716    0.00093666
    {'PostOpDay:Outcome_Unsuccessful'             }        -0.00016386    0.00080377    -0.20386    4344       0.83847    -0.0017397     0.0014119
    {'Group_Intact:PostOpDay:Outcome_Unsuccessful'}          0.0003544     0.0015386     0.23034    4344       0.81784     -0.002662     0.0033708

Random effects covariance parameters:
Group: AnimalID (10 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.074517

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        0.14439 

Result: Extension Duration

There is a significant main fixed effect of intercept (p << 0.001), indicating that an average Extension can be expected to last approximately 380-ms.
There is no significant main fixed effect of Group (p = 0.181).
There is a significant main fixed effect of Post-Op Day (p << 0.001), which indicates that there is an overall tendency to decrease reach duration by approximately 4-ms per day.
Unlike the model for Total Duration, there is no main effect of Outcome (p = 0.925).
There is a significant fixed interaction of Group*Post-Op Day (p<<0.001), which indicates that Intact rats tend to have a slight increase in extension duration (about 5-ms per day) by comparison to Ischemic rats.
There is a trending significant fixed interaction of Group*Outcome (p = 0.055), which indicates that Intact rats tend to perform Unsuccessful retrievals about 45-ms faster on average.
Unlike the model for Total Duration, there is no significant fixed interaction effect of PostOpDay*Outcome (p = 0.838).
Unlike the model for Total Duration, there is no significant three-way fixed interaction of Group*PostOpDay*Outcome (p = 0.818).

disp(glm.duration.Extension_Duration.Rsquared);
    Ordinary: 0.2769
    Adjusted: 0.2758
glm.duration.Retraction_Duration = fitglme(U,"Retraction_Duration~Group*Outcome*PostOpDay+(1|AnimalID)",'FitMethod',"Laplace");
disp(glm.duration.Retraction_Duration);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations            4290
    Fixed effects coefficients           8
    Random effects coefficients         10
    Covariance parameters                2
    Distribution                    Normal
    Link                            Identity
    FitMethod                       Laplace

Formula:
    Retraction_Duration ~ 1 + Group*PostOpDay + Group*Outcome + PostOpDay*Outcome + Group:PostOpDay:Outcome + (1 | AnimalID)

Model fit statistics:
    AIC        BIC        LogLikelihood    Deviance
    -2593.8    -2530.1    1306.9           -2613.8 

Fixed effects coefficients (95% CIs):
    Name                                                   Estimate      SE            tStat      DF      pValue         Lower         Upper     
    {'(Intercept)'                                }           0.47903       0.02174     22.034    4282    5.1694e-102       0.43641       0.52166
    {'Group_Intact'                               }         -0.051369      0.037502    -1.3698    4282        0.17083      -0.12489      0.022155
    {'PostOpDay'                                  }        -0.0083463    0.00070371     -11.86    4282     5.9829e-32    -0.0097259    -0.0069667
    {'Outcome_Unsuccessful'                       }          -0.05555      0.017596     -3.157    4282      0.0016052     -0.090048     -0.021053
    {'Group_Intact:PostOpDay'                     }          0.015703     0.0016365     9.5949    4282     1.3827e-21      0.012494      0.018911
    {'Group_Intact:Outcome_Unsuccessful'          }         -0.047157      0.028772     -1.639    4282        0.10129      -0.10357     0.0092514
    {'PostOpDay:Outcome_Unsuccessful'             }         0.0033098    0.00098932     3.3456    4282     0.00082824     0.0013703     0.0052494
    {'Group_Intact:PostOpDay:Outcome_Unsuccessful'}         -0.011929     0.0018951    -6.2949    4282     3.3851e-10     -0.015645    -0.0082139

Random effects covariance parameters:
Group: AnimalID (10 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.04542 

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        0.17776 

Result: Retraction Duration

There is a significant main fixed effect of intercept (p << 0.001), indicating that an average Retraction can be expected to last approximately 480-ms.
Therefore, on average, we expect that Retraction takes longer than Extension.
There is no significant main fixed effect of Group (p = 0.171).
There is a significant main fixed effect of Post-Op Day (p << 0.001), which indicates that there is an overall tendency to decrease Retraction duration by approximately 8-ms per day.
Therefore, on average, the reduction in Total Duration due to Retraction is approximately twice as much as the reduction during Extension.
Unlike the model for Extension, but similar to the model for Total Duration, there is a main effect of Outcome (p = 0.002), indicating that the retraction phase during Unsuccessful grasps is on average about 50-ms faster.
There is a significant fixed interaction of Group*Post-Op Day (p<<0.001), which indicates that Intact rats tend to have a slight increase in Retraction duration (about 15-ms per day) by comparison to Ischemic rats.
Unlike the models for Extension or Total Duration, there is no significant interaction of Group*Outcome (p = 0.101).
Unlike the model for Extension, but similar to the model for Total Duration, there is a significant fixed interaction effect of PostOpDay*Outcome (p < 0.001), which indicates that there is a slight increase in Retraction duration (about 3-ms per day) for Unsuccessful retrievals.
Unlike the model for Extension, but similar to the model for Total Duration, there is no significant three-way fixed interaction of Group*PostOpDay*Outcome (p << 0.001), which indicates that the Retraction becomes faster (by about 12-ms per day) for Intact rats during Unsuccessful retrievals.

disp(glm.duration.Retraction_Duration.Rsquared);
    Ordinary: 0.1342
    Adjusted: 0.1328

Statistics for Behavior Score vs Time

We should also make plots of the actual behavioral performance based on the included trial dataset, so we can see how the animals progressed over time.

Load the behavior scoring file separately

Note: there are two separate databases used for the neurophysiological rate trials (main data table is T), and the overall behavioral performance (main data table is S). We had to use a reduced subset of the total behavior, primarily because the neurophysiological rate analyses require the trials under consideration to be as stereotyped as possible: since the behavior is relatively unconstrained to begin with (by comparison with a head-fixed animal), we at least had to give the best chance of tracking a conserved behavior by requiring qualitative features of the behavior to remain consistent. That is, for neurophysiological rate analyses, a trial was required to have a Reach, Grasp, and Complete event. This means that only trials in which the rat approached the opening, reached for the pellet, closed its digits either on or near the pellet, and fully retracted its paw into the box were considered in such analyses. Attempts in which the rat "flailed" out of the box and then immediately re-attempted to capture the pellet on a "second attempt" without fully retracting into the box were discarded, as this would clearly cause an unwanted bias in the neurophysiological dynamics under consideration. Such attempts might be included with a larger group size and more trials, but given the scope of the present study, we thought it was more informative to exclude such cases.

However, in light of the manner in which rehabilitation in the reaching rat is typically characterized, we kept a separate table quantifying the total successful retrievals and attempts that were scored in total for any videos used during the metadata tagging process (the process of assigning the category of Reach, Grasp, Complete, or Support to a given video frame, along with denoting whether the pellet retrieval attempt was successful or unsuccessful). This table (S) did not exclude trials that lacked the Reach or Complete component. In this way, we include trials in which there were multiple "flails" in which the rat immediately realized it did not grasp the pellet, reached out, and then retrieved the pellet to the box. Retaining such trials makes a large difference in the overall percent of apparently successful trials: this may be due to the fact that (particularly in Intact rats), as the rat becomes more proficient (and therefore more motivated to attempt to retrieve the pellets), they sometimes develop the "bad habit" of reaching too quickly (almost as if they are impatient to retrieve the pellet). This may explain why there is an apparent increase in the duration of Successful pellet retrievals in Intact rats: they realize that they have to slow down and be a little more careful during the retrieval, or simply because of how the data is included, the Successful trials are more likely to be the ones in which they were a little more "careful" during reaching out for the pellet and that's reflected in the overall reach duration.

As a final note, these two separate data tables are ultimately derived through a different set of exclusions, but from the same scored source data. To see how the table S is obtained, refer to group.getBlockSummary() while for T refer to group.getRateTable().

if exist('S','var')==0
   S = getfield(load(defaults.files('neurophys_behavior_data_file'),'S'),'S');
end

Plot figures of behavior performance over time

[~,Sf] = analyze.behavior.outcomes(S,'BatchIndex',1); % Sf: `S` but formatted (as used for figures)

Left: There is a tendency for an increased overall number of trials in Ischemic rats as recovery proceeds. Right: Trends are not as strong for Intact rats, but there may be a slight overall decline in total number of trials as time progresses further from the electrode implant.

analyze.behavior.outcomes(S,'BatchIndex',2);

Left: Along with an increase in overall motivation, Ischemic rats have an overall tendency to have an increase in the total number of Successful pellet retrievals as the time from injury increases. Right: There is a somewhat substantial drop-off in total number of Successful pellet retrievals as time increases for the Intact group, which is not what is expected. Anecdotally, the rats become "less-patient," often waiting in front of the door and immediately flailing at the pellet as soon as it is placed with much more rapid forelimb movements. This result suggests two hypotheses for the decline in success of Intact rats: in the first alternative, the rats knock more pellets off the platform due to reaching in a less-accurate manner; in the second alternative, the Intact rats simply perform less reaches altogether. Note: these hypotheses are not mutually-exclusive; if the rat knocks the pellet off the platform, he does not have the opportunity for a second- or third-reach attempt, so the total number of available attempts would decline as well. The relative reach (within a pellet being placed on the platform) was not tagged during video metadata scoring.

analyze.behavior.outcomes(S,'BatchIndex',3);

Left: In Ischemic rats, there appear to be two "types" of trends with respect to the overall number of Unsuccessful pellet retrievals: three of the animals experience a decrease in the total number of Unsuccessful pellet retrievals as time from the injury increases, while three of the animals experience an increase in the total number of Unsuccessful pellet retrievals as time from the injury increases. Of these latter three, two of the animals demonstrate a substantially increased number of Unsuccessful pellet retrievals as time from the injury progresses; qualitatively, these animals seemed to get healthier and were quite motivated to attempt recovery of the pellet, but they often knocked the pellet off the platform in a fashion that was similar to the Intact animals. Right: An overall decline in total number of Unsuccessful attempts for Intact rats suggests that the latter hypotheses suggested previously (that there are less reaches altogether, whether due to lower motivation or knocking pellets off the platform).

analyze.behavior.outcomes(S,'BatchIndex',4);

Left: Except for one rat in the Ischemic group, injured rats tended to improve the relative ratio of successful pellet retrievals to overall pellet retrieval attempts as a function of time from the injury. Right: Rats in the Intact group tended to maintain a comparable overall ratio of successful pellet retrievals to the total number of retrieval attempts. The overall low ratio (many other studies report values from 50-80%) may be because all movements toward the pellet were counted independently as retrieval attempts. The important normalization to consider in this case is that the Intact rats maintained a comparable level of accuracy over the duration of the recovery period, while the Ischemic rats did not.

Fit Generalized Linear Model for each response plotted above

Initialize the glm.outcome struct with the original data (glm.outcome.data) used as well as the logical vector of excluded (glm.outcome.exclusions; any values that are 1). We will exclude any recordings in which any of the following conditions are met:

No Successful retrievals
No Unsuccessful retrievals
Less than 10 retrieval attempts

glm.outcome.data = Sf;
glm.outcome.exclusions = (glm.outcome.data.Successful_Trials==0) | (glm.outcome.data.Unsuccessful_Trials==0) | (glm.outcome.data.N < 10);

Fit model for trend in total number of retrievals attempted

This may be an important factor to indicate how motivated, overall, the animal was. We fit the model using the fixed effects of Intercept, Group, PostOpDay, and Group:PostOpDay, and a random intercept and coefficient term grouped by AnimalID. In Wilkinson Notation, this is denoted "Response~Group*PostOpDay + (1+PostOpDay|AnimalID)". Since our Response is a count, we fit a Poisson distribution, and use the corresponding canonical link function (log). We exclude any recordings that had either no Successful trials or no Unsuccessful trials. We fit using the Laplace method so that we can compare model likelihoods for similar models.

glm.outcome.Total_Trials = fitglme(glm.outcome.data,...
   "Total_Trials~Group*PostOpDay+(1+PostOpDay|AnimalID)",...
   'FitMethod',"Laplace",...
   'Distribution',"Poisson",...
   'Link',"log",...
   'Exclude',glm.outcome.exclusions);
disp(glm.outcome.Total_Trials);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations             170
    Fixed effects coefficients           4
    Random effects coefficients         20
    Covariance parameters                3
    Distribution                    Poisson
    Link                            Log   
    FitMethod                       Laplace

Formula:
    Total_Trials ~ 1 + Group*PostOpDay + (1 + PostOpDay | AnimalID)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    3391.1    3413.1    -1688.6          3377.1  

Fixed effects coefficients (95% CIs):
    Name                              Estimate     SE          tStat      DF     pValue        Lower        Upper     
    {'(Intercept)'           }           4.3461     0.20144     21.575    166    4.9753e-50       3.9484        4.7438
    {'Group_Intact'          }           0.5153     0.31746     1.6232    166       0.10644     -0.11148        1.1421
    {'PostOpDay'             }         0.013316    0.011851     1.1236    166       0.26279    -0.010082      0.036714
    {'Group_Intact:PostOpDay'}        -0.040857    0.018709    -2.1838    166      0.030379    -0.077795    -0.0039187

Random effects covariance parameters:
Group: AnimalID (10 Levels)
    Name1                  Name2                  Type            Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std' }         0.48733
    {'PostOpDay'  }        {'(Intercept)'}        {'corr'}        -0.80349
    {'PostOpDay'  }        {'PostOpDay'  }        {'std' }        0.028678

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        1       
disp(glm.outcome.Total_Trials.Rsquared);
    Ordinary: 0.4895
    Adjusted: 0.4803

Fit model for trend in total number of successful retrievals

Regardless of how many retrievals were attempted, it's possible that the only thing that matters is how many times the behavior was successfully executed. We fit the model using the fixed effects of Intercept, Group, PostOpDay, and Group:PostOpDay, and a random intercept and coefficient term grouped by AnimalID. In Wilkinson Notation, this is denoted "Response~Group*PostOpDay + (1+PostOpDay|AnimalID)". Since our Response is a count, we fit a Poisson distribution, and use the corresponding canonical link function (log). We exclude any recordings that had either no Successful trials or no Unsuccessful trials. We fit using the Laplace method so that we can compare model likelihoods for similar models.

glm.outcome.Successful_Trials = fitglme(glm.outcome.data,...
   "Successful_Trials~Group*PostOpDay+(1+PostOpDay|AnimalID)",...
   'FitMethod',"Laplace",...
   'Distribution',"Poisson",...
   'Link',"log",...
   'Exclude',glm.outcome.exclusions);
disp(glm.outcome.Successful_Trials);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations             170
    Fixed effects coefficients           4
    Random effects coefficients         20
    Covariance parameters                3
    Distribution                    Poisson
    Link                            Log   
    FitMethod                       Laplace

Formula:
    Successful_Trials ~ 1 + Group*PostOpDay + (1 + PostOpDay | AnimalID)

Model fit statistics:
    AIC       BIC     LogLikelihood    Deviance
    1335.1    1357    -660.54          1321.1  

Fixed effects coefficients (95% CIs):
    Name                              Estimate     SE          tStat      DF     pValue        Lower         Upper    
    {'(Intercept)'           }            2.272     0.26025       8.73    166    2.6151e-15        1.7581       2.7858
    {'Group_Intact'          }           1.2561     0.40308     3.1162    166     0.0021594       0.46024       2.0519
    {'PostOpDay'             }         0.022422    0.012996     1.7253    166      0.086326    -0.0032361     0.048081
    {'Group_Intact:PostOpDay'}        -0.050591     0.02013    -2.5132    166      0.012918     -0.090336    -0.010847

Random effects covariance parameters:
Group: AnimalID (10 Levels)
    Name1                  Name2                  Type            Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std' }         0.60562
    {'PostOpDay'  }        {'(Intercept)'}        {'corr'}        -0.84439
    {'PostOpDay'  }        {'PostOpDay'  }        {'std' }         0.02978

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        1       
disp(glm.outcome.Successful_Trials.Rsquared);
    Ordinary: 0.5663
    Adjusted: 0.5585

Fit model for trend in total number of unsuccessful retrievals

We can also get information from the raw number of unsuccessful retrievals by day. The information provided should be similar to the information provided by the Total_Trials model when interpreted in combination with the Successful_Trials or Score results, except it tells us quantitatively if there are trends in the overall number of unsuccessful reaches, as opposed to the total overall motivation. So it is just a different way of looking at the data that we might not ultimately use but if it's more convenient (e.g. based on how the discussion of the paper is framed) to describe or understand the data using this model then we will just have it ready. We fit the model using the fixed effects of Intercept, Group, PostOpDay, and Group:PostOpDay, and a random intercept term grouped by AnimalID. In Wilkinson Notation, this is denoted "Response~Group*PostOpDay + (1|AnimalID)". Since our Response is a count, we fit a Poisson distribution, and use the corresponding canonical link function (log). We exclude any recordings that had either no Successful trials or no Unsuccessful trials. We fit using the Laplace method so that we can compare model likelihoods for similar models.

glm.outcome.Unsuccessful_Trials = fitglme(glm.outcome.data,...
   "Unsuccessful_Trials~Group*PostOpDay+(1+PostOpDay|AnimalID)",...
   'FitMethod',"Laplace",...
   'Distribution',"Poisson",...
   'Link',"log",...
   'Exclude',glm.outcome.exclusions);
disp(glm.outcome.Unsuccessful_Trials);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations             170
    Fixed effects coefficients           4
    Random effects coefficients         20
    Covariance parameters                3
    Distribution                    Poisson
    Link                            Log   
    FitMethod                       Laplace

Formula:
    Unsuccessful_Trials ~ 1 + Group*PostOpDay + (1 + PostOpDay | AnimalID)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    3283.6    3305.6    -1634.8          3269.6  

Fixed effects coefficients (95% CIs):
    Name                              Estimate     SE          tStat      DF     pValue        Lower        Upper    
    {'(Intercept)'           }           4.2312     0.19459     21.744    166    1.8994e-50        3.847       4.6154
    {'Group_Intact'          }            0.311     0.30671      1.014    166       0.31206     -0.29455      0.91655
    {'PostOpDay'             }        0.0098519    0.012242    0.80474    166       0.42212    -0.014319     0.034023
    {'Group_Intact:PostOpDay'}        -0.036939    0.019344    -1.9096    166      0.057914    -0.075132    0.0012534

Random effects covariance parameters:
Group: AnimalID (10 Levels)
    Name1                  Name2                  Type            Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std' }         0.46928
    {'PostOpDay'  }        {'(Intercept)'}        {'corr'}        -0.76167
    {'PostOpDay'  }        {'PostOpDay'  }        {'std' }        0.029585

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        1       
disp(glm.outcome.Unsuccessful_Trials.Rsquared);
    Ordinary: 0.4903
    Adjusted: 0.4811

Fit model for trend in ratio of successful retrievals to total retrieval attempts

Regardless of how many retrievals were attempted, it's possible that the only thing that matters is how many times the behavior was successfully executed. To give us more flexibility in fitting the model, we will fit a primarily random effects model, then conduct coefficient tests on the estimated random effects. We fit the model using the fixed effect of Intercept (only), and a random intercept and coefficient term grouped by AnimalID (to take care of within-group heterogeneity) nested in Group.

In Wilkinson Notation, this is denoted "Successful_Trials~1+(1+PostOpDay|GroupID:AnimalID)".

Instead of fitting a Poisson distribution, since our Success_Rate is a Bernoulli variable (i.e. takes either a zero or one value for every trial), we instead fit using a related distribution, setting the Distribution parameter to Binomial. We also use the corresponding canonical link function for the Binomial distribution (logit). We fit using the Laplace method so that we can compare model likelihoods for similar models.

Our strategy is to recover the Random Effects matrix, then specify hypothesis tests for coefficient estimates that should sum to zero: we can do this for the Intercept as well as for the Coefficient term for both Group levels (Intact, Ischemia).

glm.outcome.Success_Rate = fitglme(glm.outcome.data,...

"Successful_Trials~1+(1+PostOpDay|Group:AnimalID)",...

'FitMethod',"Laplace",...

'Distribution',"Binomial",... % Note that this is different than the first three

'BinomialSize',glm.outcome.data.Total_Trials, ... % For Binomial distribution, we also need to specify the total number of trials

'Link',"logit",... % We should use the corresponding canonical link function, since optimizer is best-suited for that case

'DummyVarCoding',"full",...

'Exclude',glm.outcome.exclusions);

disp(glm.outcome.Success_Rate);

Generalized linear mixed-effects model fit by ML Model information: Number of observations 170 Fixed effects coefficients 1 Random effects coefficients 20 Covariance parameters 3 Distribution Binomial Link Logit FitMethod Laplace Formula: Successful_Trials ~ 1 + (1 + PostOpDay | Group:AnimalID) Model fit statistics: AIC BIC LogLikelihood Deviance 1221.4 1234 -606.71 1213.4 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue Lower Upper {'(Intercept)'} -1.4839 0.15091 -9.8328 169 2.5123e-18 -1.7818 -1.186 Random effects covariance parameters: Group: Group:AnimalID (10 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.56709 {'PostOpDay' } {'(Intercept)'} {'corr'} -0.55632 {'PostOpDay' } {'PostOpDay' } {'std' } 0.023881 Group: Error Name Estimate {'sqrt(Dispersion)'} 1

analyze.stat.panelized_residuals(glm.outcome.Success_Rate);

R-squared Ordinary: 0.6085 Adjusted: 0.6085

This model captures the data very well. We can verify that this is not a spurious relation by plotting the fitted residuals. We can see that the residuals histogram (top panel) isn't perfectly normal, but it is approximately "okay" and the fitted residual scatter plot (bottom panel) indicates that there is no systematic bias toward over-predicting or under-predicting at one end or the other of the line. The fitted residuals also indicate that the data is unlikely likely to be heteroscedastic, except for possibly for fitted values above 0.35, which are only covered sparsely in these data.

Recover and display the Random Effects table

[~,~,glm.outcome.randomStats] = randomEffects(glm.outcome.Success_Rate);
disp(glm.outcome.randomStats);
    Random effect coefficients: DFMethod = 'residual', Alpha = 0.05

    Group                     Level                     Name                   Estimate      SEPred       tStat        DF     pValue        Lower        Upper     
    {'Group:AnimalID'}        {'Ischemia RC-02'}        {'(Intercept)'}          -0.36963      0.23396      -1.5799    169       0.11601      -0.8315      0.092237
    {'Group:AnimalID'}        {'Ischemia RC-02'}        {'PostOpDay'  }           0.05125    0.0096701       5.2998    169    3.5799e-07      0.03216       0.07034
    {'Group:AnimalID'}        {'Ischemia RC-04'}        {'(Intercept)'}          -0.75779      0.22435      -3.3777    169    0.00090742      -1.2007      -0.31489
    {'Group:AnimalID'}        {'Ischemia RC-04'}        {'PostOpDay'  }          0.022008    0.0098305       2.2388    169      0.026475     0.002602      0.041415
    {'Group:AnimalID'}        {'Ischemia RC-05'}        {'(Intercept)'}        -0.0057139      0.20242    -0.028228    169       0.97751      -0.4053       0.39388
    {'Group:AnimalID'}        {'Ischemia RC-05'}        {'PostOpDay'  }        -0.0086102    0.0078271      -1.1001    169       0.27287    -0.024062     0.0068413
    {'Group:AnimalID'}        {'Ischemia RC-08'}        {'(Intercept)'}          -0.49975      0.23034      -2.1697    169      0.031427     -0.95446      -0.04505
    {'Group:AnimalID'}        {'Ischemia RC-08'}        {'PostOpDay'  }         -0.009265    0.0086336      -1.0731    169       0.28474    -0.026309     0.0077785
    {'Group:AnimalID'}        {'Ischemia RC-26'}        {'(Intercept)'}          -0.30312      0.21889      -1.3848    169       0.16794     -0.73523       0.12899
    {'Group:AnimalID'}        {'Ischemia RC-26'}        {'PostOpDay'  }        -0.0055906    0.0097743     -0.57197    169        0.5681    -0.024886      0.013705
    {'Group:AnimalID'}        {'Ischemia RC-30'}        {'(Intercept)'}          -0.70934      0.33912      -2.0917    169      0.037961      -1.3788     -0.039878
    {'Group:AnimalID'}        {'Ischemia RC-30'}        {'PostOpDay'  }          0.016804     0.018702      0.89852    169       0.37019    -0.020116      0.053724
    {'Group:AnimalID'}        {'Intact RC-14'  }        {'(Intercept)'}           0.18456      0.20383      0.90546    169       0.36651     -0.21782       0.58695
    {'Group:AnimalID'}        {'Intact RC-14'  }        {'PostOpDay'  }         -0.011486    0.0098862      -1.1619    169       0.24693    -0.031003     0.0080299
    {'Group:AnimalID'}        {'Intact RC-18'  }        {'(Intercept)'}             0.802      0.20882       3.8407    169    0.00017328      0.38977        1.2142
    {'Group:AnimalID'}        {'Intact RC-18'  }        {'PostOpDay'  }        -0.0072763     0.011055      -0.6582    169       0.51131      -0.0291      0.014547
    {'Group:AnimalID'}        {'Intact RC-21'  }        {'(Intercept)'}           0.80194      0.18847        4.255    169    3.4557e-05      0.42988         1.174
    {'Group:AnimalID'}        {'Intact RC-21'  }        {'PostOpDay'  }         -0.015751    0.0069133      -2.2784    169      0.023958    -0.029398    -0.0021034
    {'Group:AnimalID'}        {'Intact RC-43'  }        {'(Intercept)'}            0.1336      0.17351      0.76997    169       0.44239     -0.20893       0.47612
    {'Group:AnimalID'}        {'Intact RC-43'  }        {'PostOpDay'  }           0.02459    0.0051287       4.7946    169    3.5564e-06     0.014465      0.034715

Recover and display the Random Effects Covariance matrix for the Experimental Group grouped effect

[~,~,glm.outcome.covarianceStats] = covarianceParameters(glm.outcome.Success_Rate);
disp(glm.outcome.covarianceStats{1});
    Covariance Type: FullCholesky

    Group          Name1                  Name2                  Type            Estimate    Lower       Upper   
    [1x14 char]    {'(Intercept)'}        {'(Intercept)'}        {'std' }         0.56709     0.33945     0.94737
    [1x14 char]    {'PostOpDay'  }        {'(Intercept)'}        {'corr'}        -0.55632    -0.88065     0.12304
    [1x14 char]    {'PostOpDay'  }        {'PostOpDay'  }        {'std' }        0.023881    0.013858    0.041154

The covariance matrix tells us that on average the intercept and slope terms for our random effects are approximately uncorrelated (the 95%-confidence bound for the interaction between Intercept and PostOpDay contains zero, indicating that this effect is not significant at that level). However, both the Intercept and PostOpDay terms do not contain zero in their confidence interval bounds, indicating that they are both significant random effects.

Conduct coefficient test to compare Intercept and PostOpDay by Group

The last step is to identify the term indices that correspond to matched Intercept or PostOpDay parameters, by Group, then use those terms to conduct a hypothesis test against the null hypothesis that the sum of coefficients for variables in those groupings

First, design the contrasts matrices. Since there are 4 animals in the Intact group and 6 animals in the Ischemia group, we will swap the assignment of those values, setting the Ischemia design matrix elements to -4 and the Intact design matrix elements to 6 (so that they should sum to zero).

K_Int = zeros(1,length(glm.outcome.randomStats));
K_Day = zeros(1,length(glm.outcome.randomStats)); 
index = struct('Intercept',struct,'PostOpDay',struct);
index.Intercept.Intact= string(glm.outcome.randomStats.Name) == "(Intercept)" & startsWith(string(glm.outcome.randomStats.Level),"Intact");
index.PostOpDay.Intact = (string(glm.outcome.randomStats.Name) == "PostOpDay") & startsWith(string(glm.outcome.randomStats.Level),"Intact");
index.Intercept.Ischemia = string(glm.outcome.randomStats.Name) == "(Intercept)" & startsWith(string(glm.outcome.randomStats.Level),"Ischemia");
index.PostOpDay.Ischemia = (string(glm.outcome.randomStats.Name) == "PostOpDay") & startsWith(string(glm.outcome.randomStats.Level),"Ischemia");
K_Int(index.Intercept.Intact) = 6;
K_Day(index.PostOpDay.Intact) = 6;
K_Int(index.Intercept.Ischemia) = -4;
K_Day(index.PostOpDay.Ischemia) = -4;

Next, conduct the hypothesis tests.

[glm.result.Intercept.pVal,glm.result.Intercept.F,glm.result.Intercept.DF1,glm.result.Intercept.DF2] = coefTest(glm.outcome.Success_Rate,0,0,'REContrast',K_Int,'DFMethod','residual');
[glm.result.PostOpDay.pVal,glm.result.PostOpDay.F,glm.result.PostOpDay.DF1,glm.result.PostOpDay.DF2] = coefTest(glm.outcome.Success_Rate,0,0,'REContrast',K_Day,'DFMethod','residual');

And associate the grouped estimates.

glm.result.Intercept.Intact = sum(glm.outcome.randomStats.Estimate(index.Intercept.Intact));
glm.result.Intercept.Ischemia = sum(glm.outcome.randomStats.Estimate(index.Intercept.Ischemia));
glm.result.PostOpDay.Intact = sum(glm.outcome.randomStats.Estimate(index.PostOpDay.Intact));
glm.result.PostOpDay.Ischemia = sum(glm.outcome.randomStats.Estimate(index.PostOpDay.Ischemia));

Last, display the results.

disp(glm.result.Intercept);
        pVal: 6.0606e-16
           F: 80.1594
         DF1: 1
         DF2: 169
      Intact: 1.9221
    Ischemia: -2.6453

The p-value for Intercept is statistically significant, indicating that we reject the null hypothesis that the random coefficients for Intercept sum to zero. This indicates that by Group, there is a significant difference in the average intercept. On average, there is an expected log-odds of 4.5 when comparing the likelihood of a Successful retrieval in Intact rats compared to the likelihood of a Successful retrieval in Ischemia rats.

disp(glm.result.PostOpDay);
        pVal: 0.0329
           F: 4.6257
         DF1: 1
         DF2: 169
      Intact: -0.0099
    Ischemia: 0.0666

The p-value for PostOpDay is also statistically significant at the alpha=0.05 level (p-value = 0.03), indicating that we reject the null hypothesis that the log-odds change by day is the same for both groups. The Ischemia group has an increased log-odds of 0.066 per day, indicating that there is a substantial improvement over the course of recovery for this group.

Fit model using PostOpDay nested in Animal as Grouping for Random Effect of Intercept

If we want to see when this effect is significant, we need a different model.

We can specify a Binomial Generalized Linear Mixed Effects model (as before), but this time we will change the model specification so that we include PostOpDay (Day) nested by Group as an ordinal grouping for a random Intercept term. We can then conduct hypothesis tests in a similar manner as the previous model, but using a grouping requirement for PostOpDay, which allows us to qualitatively split the model by "phases" of recovery (for example, by Week) so we can make a general statement about when this differentiation is the strongest.

First, update the data so that we have an Ordinal variable for day.

glm.day.data = Sf;
glm.day.data.Day = ordinal(glm.day.data.PostOpDay);
glm.day.exclusions = glm.outcome.exclusions;

Next, construct the model. Our model specification (using Wilkinson Notation) is "Successful_Trials~1+AnimalID+(1|Group:Day)", which gives us a fixed expected successful trial rate, a fixed expected successful trial rate by AnimalID, and a random intercept term for each PostOpDay nested in the two Group levels.

glm.day.model = fitglme(glm.day.data,...
   "Successful_Trials~1+AnimalID+(1|Group:Day)",...    
   'FitMethod',"Laplace",...
   'Distribution',"Binomial",...                 % Note that this is different than the first three
   'BinomialSize',glm.day.data.Total_Trials, ... % For Binomial distribution, we also need to specify the total number of trials
   'Link',"logit",...                            % We should use the corresponding canonical link function, since optimizer is best-suited for that case
   'DummyVarCoding',"effects",...
   'Exclude',glm.day.exclusions);
disp(glm.day.model);
Generalized linear mixed-effects model fit by ML

Model information:
    Number of observations             170
    Fixed effects coefficients          10
    Random effects coefficients         56
    Covariance parameters                1
    Distribution                    Binomial
    Link                            Logit 
    FitMethod                       Laplace

Formula:
    Successful_Trials ~ 1 + AnimalID + (1 | Group:Day)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    1155.1    1189.6    -566.56          1133.1  

Fixed effects coefficients (95% CIs):
    Name                      Estimate    SE          tStat      DF     pValue        Lower       Upper   
    {'(Intercept)'   }         -1.4402    0.047716    -30.184    160    6.0082e-68     -1.5345      -1.346
    {'AnimalID_RC-02'}         0.55755    0.076807      7.259    160    1.5982e-11     0.40586     0.70923
    {'AnimalID_RC-04'}        -0.42403    0.082042    -5.1684    160    6.9552e-07    -0.58605      -0.262
    {'AnimalID_RC-05'}        -0.25318    0.070364    -3.5981    160    0.00042683    -0.39214    -0.11422
    {'AnimalID_RC-08'}         -0.7756    0.073334    -10.576    160    3.7203e-20    -0.92043    -0.63078
    {'AnimalID_RC-14'}        0.021218    0.081959    0.25888    160       0.79606    -0.14064     0.18308
    {'AnimalID_RC-18'}         0.64686    0.082564     7.8347    160    6.1556e-13     0.48381     0.80992
    {'AnimalID_RC-21'}         0.60834    0.078297     7.7697    160     8.939e-13     0.45371     0.76297
    {'AnimalID_RC-26'}        -0.37388    0.077574    -4.8197    160    3.3176e-06    -0.52709    -0.22068
    {'AnimalID_RC-30'}         -0.3786     0.10616    -3.5663    160    0.00047779    -0.58826    -0.16895

Random effects covariance parameters:
Group: Group:Day (56 Levels)
    Name1                  Name2                  Type           Estimate
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.30014 

Group: Error
    Name                        Estimate
    {'sqrt(Dispersion)'}        1       

Since this model was fit using the "effects" specification for DummyVarCoding, it means that the model is fit such that the coefficient estimates for each member of AnimalID must sum to zero, excluding RC-43, which is set to zero by default to reduce the number of degrees of freedom. The low p-value for each coefficient value tells us that almost every Rat, with the exception of RC-14, has a substantially different expected performance level (as we might anticipate). To assess differences by Group, we need to look at the Random Effects. Before doing that, let's check to see if this model is a reasonable approximation of the data, or if the previous Random Effects model is a much better fit (in which case, we should simply stop here).

analyze.stat.panelized_residuals(glm.day.model);

R-squared Ordinary: 0.7370 Adjusted: 0.7222

This model has an even better fit of residuals than the previous model, and that's reflect in its R-squared ordinary and adjusted values. We can also compare the generalized linear mixed effects models based on model likelihood since both were estimated using the Laplace method.

glm.comparison = compare(glm.outcome.Success_Rate,glm.day.model);
disp(glm.comparison);
    Theoretical Likelihood Ratio Test

    Model      DF    AIC       BIC       LogLik     LRStat    deltaDF    pValue   
    GLME        4    1221.4      1234    -606.71                                  
    ALTGLME    11    1155.1    1189.6    -566.56    80.293    7          1.199e-14

The low p-value indicates that we should reject glm.outcome.Success_Rate in favor of the alternative model, glm.day.model.

[~,~,glm.day.randomStats] = randomEffects(glm.day.model);
disp(glm.day.randomStats);
    Random effect coefficients: DFMethod = 'residual', Alpha = 0.05

    Group                Level                  Name                   Estimate     SEPred     tStat        DF     pValue       Lower         Upper      
    {'Group:Day'}        {'Ischemia 3' }        {'(Intercept)'}        -0.069173     0.2138     -0.32353    160      0.74671      -0.49141        0.35307
    {'Group:Day'}        {'Ischemia 4' }        {'(Intercept)'}          0.49616    0.14992       3.3095    160    0.0011547       0.20009        0.79224
    {'Group:Day'}        {'Ischemia 5' }        {'(Intercept)'}           0.6846    0.20534        3.334    160    0.0010639       0.27907         1.0901
    {'Group:Day'}        {'Ischemia 6' }        {'(Intercept)'}          0.27018    0.13863        1.949    160     0.053044    -0.0035897        0.54396
    {'Group:Day'}        {'Ischemia 7' }        {'(Intercept)'}         -0.30694     0.1676      -1.8314    160       0.0689      -0.63793        0.02405
    {'Group:Day'}        {'Ischemia 8' }        {'(Intercept)'}         -0.10729     0.1406     -0.76307    160      0.44654      -0.38496        0.17038
    {'Group:Day'}        {'Ischemia 9' }        {'(Intercept)'}         -0.37546    0.16066       -2.337    160      0.02068      -0.69276       -0.05817
    {'Group:Day'}        {'Ischemia 10'}        {'(Intercept)'}         -0.16157    0.18515     -0.87265    160      0.38416      -0.52722        0.20408
    {'Group:Day'}        {'Ischemia 11'}        {'(Intercept)'}         -0.68209    0.18922      -3.6047    160     0.000417       -1.0558       -0.30839
    {'Group:Day'}        {'Ischemia 12'}        {'(Intercept)'}          0.08668    0.16537      0.52415    160       0.6009      -0.23992        0.41328
    {'Group:Day'}        {'Ischemia 13'}        {'(Intercept)'}         -0.14887    0.12186      -1.2217    160      0.22363      -0.38954       0.091789
    {'Group:Day'}        {'Ischemia 14'}        {'(Intercept)'}         -0.41209    0.14632      -2.8163    160    0.0054698      -0.70106       -0.12311
    {'Group:Day'}        {'Ischemia 15'}        {'(Intercept)'}         -0.19008    0.14246      -1.3343    160      0.18401      -0.47143       0.091265
    {'Group:Day'}        {'Ischemia 16'}        {'(Intercept)'}         -0.12395    0.12365      -1.0025    160      0.31763      -0.36814        0.12024
    {'Group:Day'}        {'Ischemia 17'}        {'(Intercept)'}         -0.37032     0.1432      -2.5861    160     0.010599      -0.65313      -0.087518
    {'Group:Day'}        {'Ischemia 18'}        {'(Intercept)'}         -0.13442    0.11241      -1.1957    160      0.23356      -0.35642       0.087587
    {'Group:Day'}        {'Ischemia 19'}        {'(Intercept)'}          0.16027     0.1153         1.39    160      0.16646     -0.067441        0.38799
    {'Group:Day'}        {'Ischemia 20'}        {'(Intercept)'}         0.064277    0.12129      0.52996    160      0.59688      -0.17525        0.30381
    {'Group:Day'}        {'Ischemia 21'}        {'(Intercept)'}          0.19479    0.14496       1.3437    160      0.18095       -0.0915        0.48107
    {'Group:Day'}        {'Ischemia 22'}        {'(Intercept)'}          0.12171    0.15153      0.80318    160      0.42306      -0.17756        0.42098
    {'Group:Day'}        {'Ischemia 23'}        {'(Intercept)'}          0.27139    0.17563       1.5452    160      0.12426      -0.07546        0.61824
    {'Group:Day'}        {'Ischemia 24'}        {'(Intercept)'}         0.027281    0.17577      0.15521    160      0.87685      -0.31984         0.3744
    {'Group:Day'}        {'Ischemia 25'}        {'(Intercept)'}           0.2781    0.14457       1.9237    160     0.056169     -0.007404         0.5636
    {'Group:Day'}        {'Ischemia 26'}        {'(Intercept)'}         0.064858    0.14382      0.45095    160      0.65264      -0.21918         0.3489
    {'Group:Day'}        {'Ischemia 27'}        {'(Intercept)'}          -0.0645    0.14683     -0.43927    160      0.66106      -0.35448        0.22548
    {'Group:Day'}        {'Ischemia 28'}        {'(Intercept)'}          0.11183    0.14563      0.76787    160      0.44369      -0.17578        0.39944
    {'Group:Day'}        {'Ischemia 29'}        {'(Intercept)'}          0.31429     0.1429       2.1994    160     0.029284      0.032081         0.5965
    {'Group:Day'}        {'Ischemia 30'}        {'(Intercept)'}          0.14677    0.23752      0.61793    160       0.5375      -0.32231        0.61585
    {'Group:Day'}        {'Intact 3'   }        {'(Intercept)'}         -0.11329    0.14481     -0.78237    160      0.43515      -0.39928        0.17269
    {'Group:Day'}        {'Intact 4'   }        {'(Intercept)'}         -0.23187    0.12806      -1.8106    160     0.072082      -0.48478       0.021044
    {'Group:Day'}        {'Intact 5'   }        {'(Intercept)'}        -0.099588    0.12052     -0.82633    160      0.40985       -0.3376        0.13842
    {'Group:Day'}        {'Intact 6'   }        {'(Intercept)'}          0.22387    0.14104       1.5873    160      0.11442     -0.054668        0.50241
    {'Group:Day'}        {'Intact 7'   }        {'(Intercept)'}         0.033481    0.15877      0.21088    160      0.83325      -0.28008        0.34704
    {'Group:Day'}        {'Intact 8'   }        {'(Intercept)'}         -0.15012    0.11125      -1.3494    160      0.17913      -0.36984       0.069594
    {'Group:Day'}        {'Intact 9'   }        {'(Intercept)'}        -0.069568    0.11772     -0.59096    160      0.55538      -0.30205        0.16292
    {'Group:Day'}        {'Intact 10'  }        {'(Intercept)'}          0.33698    0.12949       2.6024    160     0.010127      0.081253        0.59271
    {'Group:Day'}        {'Intact 11'  }        {'(Intercept)'}         0.034845    0.13287      0.26225    160      0.79347      -0.22756        0.29725
    {'Group:Day'}        {'Intact 12'  }        {'(Intercept)'}         0.010598    0.13443     0.078836    160      0.93726      -0.25489        0.27609
    {'Group:Day'}        {'Intact 13'  }        {'(Intercept)'}          0.36571    0.13366       2.7361    160    0.0069196       0.10174        0.62967
    {'Group:Day'}        {'Intact 14'  }        {'(Intercept)'}         0.096041    0.13981      0.68693    160      0.49312      -0.18007        0.37216
    {'Group:Day'}        {'Intact 15'  }        {'(Intercept)'}         -0.41323    0.15845       -2.608    160    0.0099684      -0.72615       -0.10032
    {'Group:Day'}        {'Intact 16'  }        {'(Intercept)'}         -0.30286    0.14829      -2.0424    160     0.042751      -0.59571      -0.010013
    {'Group:Day'}        {'Intact 17'  }        {'(Intercept)'}        -0.015227    0.16684    -0.091271    160      0.92739      -0.34471        0.31426
    {'Group:Day'}        {'Intact 18'  }        {'(Intercept)'}          0.15898    0.12433       1.2786    160      0.20287     -0.086567        0.40452
    {'Group:Day'}        {'Intact 19'  }        {'(Intercept)'}         -0.30174     0.1942      -1.5538    160      0.12222      -0.68526       0.081785
    {'Group:Day'}        {'Intact 20'  }        {'(Intercept)'}          0.17876    0.14096       1.2682    160      0.20659     -0.099623        0.45714
    {'Group:Day'}        {'Intact 21'  }        {'(Intercept)'}          0.25241    0.13252       1.9046    160     0.058621    -0.0093114        0.51414
    {'Group:Day'}        {'Intact 23'  }        {'(Intercept)'}          0.17963    0.16112       1.1149    160      0.26657      -0.13857        0.49782
    {'Group:Day'}        {'Intact 24'  }        {'(Intercept)'}          0.16625     0.1598       1.0403    160      0.29976      -0.14935        0.48184
    {'Group:Day'}        {'Intact 25'  }        {'(Intercept)'}         0.055324    0.15492      0.35711    160      0.72148      -0.25062        0.36127
    {'Group:Day'}        {'Intact 26'  }        {'(Intercept)'}          0.31224    0.17918       1.7426    160     0.083318     -0.041616        0.66609
    {'Group:Day'}        {'Intact 27'  }        {'(Intercept)'}          0.22551    0.17496       1.2889    160      0.19929      -0.12002        0.57103
    {'Group:Day'}        {'Intact 28'  }        {'(Intercept)'}           0.2138     0.1806       1.1839    160      0.23823      -0.14286        0.57047
    {'Group:Day'}        {'Intact 29'  }        {'(Intercept)'}          -0.3222    0.19728      -1.6332    160      0.10438      -0.71181       0.067401
    {'Group:Day'}        {'Intact 30'  }        {'(Intercept)'}          -0.3222    0.19728      -1.6332    160      0.10438      -0.71181       0.067401
    {'Group:Day'}        {'Intact 31'  }        {'(Intercept)'}         -0.40312    0.20404      -1.9757    160     0.049911      -0.80608    -0.00015801

What does the improved model tell us about these data?

This model can test a few different hypotheses. The first hypothesis to test is whether the two levels of Group are significantly different (we assume they are based on the previous model, but since we have updated the model we should check that the new model agrees with that result). To do that, we will identify the Random Effect terms that belong to each Group level, then create a design matrix for combining the terms with the null hypothesis that the grouped coefficients should sum to zero. Since there are an equal number of coefficient estimates for the two levels of Group, we can simply assign the elements of the group design matrix corresponding to Ischemia level elements values of -1, and those corresponding to Intact level elements values of 1.

To get an intuition of when the two groups diverge significantly on behavioral success rate, we need to conduct a different set of hypothesis tests. We will create three design matrices, for days in Week-1, days in Week-2, and days in Week-3. Depending on how many coefficients correspond to values at those levels, we may need to adjust the offsetting design matrix coefficients to account for differences in number by Group, but otherwise follow a similar strategy as outlined above. Last, we will test all of these hypotheses simultaneously using the GeneralizedLinearMixedModel.coefTest method.

K = zeros(4,length(glm.day.randomStats));
index.day.Intact   = startsWith(string(glm.day.randomStats.Level),"Intact");
index.day.Ischemia = startsWith(string(glm.day.randomStats.Level),"Ischemia");
index.day.Week1    = endsWith(string(glm.day.randomStats.Level)," 3") | endsWith(string(glm.day.randomStats.Level)," 4") | ...
   endsWith(string(glm.day.randomStats.Level)," 5") | endsWith(string(glm.day.randomStats.Level)," 6") | endsWith(string(glm.day.randomStats.Level)," 7") | ...
   endsWith(string(glm.day.randomStats.Level)," 8") | endsWith(string(glm.day.randomStats.Level)," 9");
index.day.Week2 = endsWith(string(glm.day.randomStats.Level),"10") | endsWith(string(glm.day.randomStats.Level),"11") | ...
   endsWith(string(glm.day.randomStats.Level),"12") | endsWith(string(glm.day.randomStats.Level),"13") | endsWith(string(glm.day.randomStats.Level),"14") | ...
   endsWith(string(glm.day.randomStats.Level),"15") | endsWith(string(glm.day.randomStats.Level),"16");
index.day.Week3 = endsWith(string(glm.day.randomStats.Level),"17") | endsWith(string(glm.day.randomStats.Level),"18") | ...
   endsWith(string(glm.day.randomStats.Level),"19") | endsWith(string(glm.day.randomStats.Level),"20") | endsWith(string(glm.day.randomStats.Level),"21") | ...
   endsWith(string(glm.day.randomStats.Level),"22") | endsWith(string(glm.day.randomStats.Level),"23");
K(1,index.day.Intact)   = 1;
K(1,index.day.Ischemia) = -1;
K(2,index.day.Intact & index.day.Week1)   =  1;
K(2,index.day.Ischemia & index.day.Week1) = -1;
K(3,index.day.Intact & index.day.Week2)   =  1;
K(3,index.day.Ischemia & index.day.Week2) = -1;
K(4,index.day.Intact & index.day.Week3)   =  1;
K(4,index.day.Ischemia & index.day.Week3) = -1;

Test all hypotheses simultaneously to see if there is any hypothesis that rejects the null hypothesis:

glm.day.result = struct('pVal',cell(5,1),'F',cell(5,1),'DF1',cell(5,1),'DF2',cell(5,1)); % Initialize struct array
[glm.day.result(1).pVal,glm.day.result(1).F,glm.day.result(1).DF1,glm.day.result(1).DF2] = coefTest(glm.day.model,zeros(4,10),zeros(4,1),'REContrast',K,'DFMethod','residual');
disp(glm.day.result(1));
    pVal: 0.0032
       F: 4.1454
     DF1: 4
     DF2: 160

The combined test indicates that yes, there is a significant difference somewhere among these hypothesis tests. We can conduct individual hypothesis tests.

Test the hypothesis that the coefficients grouped by the two levels of Group sum to zero

[glm.day.result(2).pVal,glm.day.result(2).F,glm.day.result(2).DF1,glm.day.result(2).DF2] = coefTest(glm.day.model,zeros(1,10),zeros(1,1),'REContrast',K(1,:),'DFMethod','residual');
disp(glm.day.result(2));
    pVal: 0.9833
       F: 4.3873e-04
     DF1: 1
     DF2: 160

Contrary to what was found on the basis of the previous model, we do not reject the null hypothesis that the intercept coefficient estimates grouped by Intact and Ischemia sum to zero, after accounting for the intercept offset to performanced that is introduced for the different levels of the Fixed Effect of AnimalID.

What does this mean? Since AnimalID is used as a fixed effect, we can think of this as removing the average pellet retrieval success rate of each rat from the picture, leaving us only to consider the "scale" of recovery for each animal. Therefore one way to think about this is that the overall dispersion of performance is not statistically significant in terms of the difference between Intact and Ischemia rats, when grouped over all the days.

Test the hypothesis that the coefficients grouped by the two levels of Group during Week-1 sum to zero

[glm.day.result(3).pVal,glm.day.result(3).F,glm.day.result(3).DF1,glm.day.result(3).DF2] = coefTest(glm.day.model,zeros(1,10),zeros(1,1),'REContrast',K(2,:),'DFMethod','residual');
disp(glm.day.result(3));
    pVal: 0.1749
       F: 1.8571
     DF1: 1
     DF2: 160

Similar to the by-group hypothesis test, we check the Week-1 test on a per-group basis. Again, we do not reject the null hypothesis that the random intercept coefficient estimates that correspond to a by-day difference in performance between the two groups during this time period sum to zero. This also makes some sense: both groups undergo a (relatively) massive surgical procedure and are put on buprenorphine to aid in their recovery. A side-effect of buprenorphine is that the rats lose their appetite, which anecdotally is due to them becoming constipated. So during the first week, in addition to the effects of buprenorphine altering how they might typically act, there is also the consideration that they are adjusting to being tethered and having an implant on top of their skull that wasn't there before. It makes sense that behavior could be comparable between groups at this time; generally, it might be the case that this period introduces an increased variance in performance that means that even if the Intact animals perform better during this time on average, the group size is underpowered to detect such a difference.

Test the hypothesis that the coefficients grouped by the two levels of Group during Week-2 sum to zero

[glm.day.result(4).pVal,glm.day.result(4).F,glm.day.result(4).DF1,glm.day.result(4).DF2] = coefTest(glm.day.model,zeros(1,10),zeros(1,1),'REContrast',K(3,:),'DFMethod','residual');
disp(glm.day.result(4));
    pVal: 0.0169
       F: 5.8300
     DF1: 1
     DF2: 160

In our Week-2 grouping, which starts on PostOpDay==10, we reject the null hypothesis that the random intercept coefficient estimates sum to zero on a day-to-day basis when split by Group.

glm.day.Week2.Intact = sum(glm.day.randomStats.Estimate(index.day.Intact & index.day.Week2));
glm.day.Week2.Ischemia = sum(glm.day.randomStats.Estimate(index.day.Ischemia & index.day.Week2));
disp(glm.day.Week2);
      Intact: 0.1281
    Ischemia: -1.6320

We note that the coefficient estimates indicate that the log-odds of a successful pellet retrieval during this time epoch are much higher for animals in the Intact level of Group. This corresponds to what we expect to see: that there is a time-dependent difference in recovery between the Intact and Ischemia groups.

Test the hypothesis that the coefficients grouped by the two levels of Group during Week-3 sum to zero

[glm.day.result(5).pVal,glm.day.result(5).F,glm.day.result(5).DF1,glm.day.result(5).DF2] = coefTest(glm.day.model,zeros(1,10),zeros(1,1),'REContrast',K(4,:),'DFMethod','residual');
disp(glm.day.result(5));
    pVal: 0.8328
       F: 0.0447
     DF1: 1
     DF2: 160

By the third week, we no longer reject the null hypothesis that the random intercept coefficient estimates sum to zero on a day-to-day basis when split by Group.

These results indicate that on the basis of behavioral differences, results focusing on the neurophysiological divergence between the two groups should focus on Post-Op Day-10 to Post-Op Day-16 (at least for the cohort of animals in this study). Similarly, it may be possible to compare activity during Week-1 (although due to the assumption that differences are not observed during this period due to an increased variance during this epoch, this may be unadvisable), or Week-3 (which seems like a more reasonable time to compare dynamics or multi-unit tuning).