Modeling Individual Recovery after Peripheral Nerve Injury in Rats and the Effects of Parturition

PREVENTING chronic pain after physical trauma, including that of major surgery, is an active area of research and drug development. Both research and development can be guided by clinical observations, most typically from genetic predictors of chronic pain. We recently observed that an environmental and biologic stimulus—childbirth—can affect the response to physical trauma. As such, chronic pain from delivery, including complicated vaginal delivery and surgical cesarean delivery, is remarkably rare.¹  In rats, spinal nerve ligation, a surgical model of neuropathic pain, results in less sustained mechanical hypersensitivity when performed at the time of parturition than in virgin females, suggesting a similar protective phenomenon occurs in rodents.² 

In this study, we advance our previous efforts by applying two novel approaches to the study of recovery from hypersensitivity behaviors in rats after surgical trauma. First, we use partial spinal nerve ligation (pSNL), a recently described model in which hypersensitivity resolves slowly over 2 to 3 months, with large interanimal variability, similar in many ways to the time course of and interpatient variability in resolution of pain after major surgery. We hypothesize that the postpartum period speeds recovery after pSNL and that this requires the presence of pups. In addition, we explored whether there is a sex difference in time course of recovery from pSNL, because women have a slightly higher incidence of chronic pain after surgery than men.

In laboratory studies, pain behavior over time is typically compared among experimental groups using repeated-measures ANOVA (or one-way analysis of covariance)³  for data that satisfy parametric assumptions or a combination of Friedman rank test and Kruskal–Wallis test for nonparametric data. Although quite useful, these methods are highly sensitive to unbalanced/unequal time points between subjects, missing values, and the violation of the underlying assumptions of these tests, especially in small laboratory samples. Traditionally applied fixed-effect ANOVA methods can elegantly examine pairwise contrasts (e.g., groups × time), or even polynomial contrasts, but cannot describe the degree of individual differences in the amount of change, or easily define the duration of treatment effects using traditional null hypothesis tests.

An alternative to ANOVA, growth curve modeling has been used extensively in large-scale clinical studies to analyze longitudinal data^4–7  but only rarely in small-n studies.⁸  In many situations, because of their flexible nature, growth curve models may offer several distinct advantages over that of ANOVA. Growth curve models allow the inspection of interindividual variability in the examination of within-subject change.⁹  Because the assumptions of growth curve modeling are more flexible than traditional methods, data may be partially missing, time points can be unequally spaced, and the distributions of outcomes and their repeated-measure covariances can be uniquely specified. Also, growth curve models typically have higher statistical power than traditional methods, which can lead to decreased expense and a reduction in the number of patients or animals needed to test a hypothesis.¹⁰  A secondary purpose of this study was to test the application of growth curve modeling to examine the trajectory of recovery of pSNL over time.

Forty-one Sprague–Dawley rats (15 to 16 weeks old; 29 females and 12 males) were acquired from Harlan Industries (Indianapolis, IN). Four groups were studied: postpartum-with-pups (n = 10), postpartum-separated-from-pups (n = 9), virgin females (n = 10), and males (n = 12). These four groups underwent pSNL surgery and were euthanized at 24 to 25 weeks of age. Pups were separated from their mothers after delivery in the postpartum-separated-from-pups group. Pups were housed with their dams after delivery in the postpartum-with-pups group until weaning at 21 days postsurgery (which is being described as the intervention). Animals were housed under a 12-h light–dark cycle with food and water ad libitum. All experiments were approved by the Animal Care and Use Committee at Wake Forest University (Winston Salem, North Carolina). Animals were studied in two cohorts—female versus male and postpartum with and postpartum separated from pups. Assignment to postpartum group was randomized, and investigators were blinded to group by having a different individual place the animals in the testing environment. No data were excluded from analysis.

Partial spinal nerve ligation surgery was performed as previously described¹¹  within 24 h of delivery in the postpartum-with-pups and the postpartum-separated-from-pups groups. Surgery timing for the virgin females and males was age matched with the postpartum groups. On the day of pSNL surgery, rats were anesthetized with 2% isoflurane in oxygen and a 3-cm incision was placed along the right dorsal surface near the spine under aseptic conditions, penetrating underlying muscles. The sixth lumbar transverse process was removed and the dorsal half of the L5 nerve was ligated using 8-0 silk. Muscles and skin were closed in separate layers. After surgery, animals were housed individually in plastic cages in a climate-controlled room under a 12-h light–dark cycle with free access to food and water.

Withdrawal threshold to punctate mechanical stimulation was determined before and after pSNL surgery in all four groups by the application of calibrated von Frey filaments (Stoelting, Wood Dale, IL) to the hind paw. Animals were separated from pups, if present, and placed on a plastic mesh floor in individual clear plastic boxes and allowed to accommodate to their environment for at least 30 min. Filaments were applied to the bending point for 5 s, and a brisk paw withdrawal was considered a positive response. Withdrawal threshold was determined using an up–down statistical method.¹²  Behavioral testing was performed before surgery and commencing the day after surgery through postsurgery day 70. Because of the variability in day of the week of delivery of pups and behavioral testing largely restricted to weekdays, the timing of the withdrawal threshold measurements varied among animals in the postpartum groups. The person performing the behavioral testing was the same for all cohorts and was blinded to surgery and treatment, but not to sex.

In a growth curve model with experimental groups or interventions, the data can be viewed as a hierarchical structure where lower-level units consisting of within-subject repeated measurements (level 1) are nested in higher-level units such as animals or experimental group (level 2). Covariates that account for within-individual changes that occur across measurements are incorporated into level 1 (time-varying covariates). The level 2 equations model each of the level 1 parameters (e.g., intercepts and slopes as outcomes) using time-invariant covariates (e.g., intervention group).

The appendix details the formulation of these models, including references. Following the steps set forth for the building of a generic growth curve model outlined, plots for each subject with the time variable on the x-axis and withdrawal threshold on the y-axis, superimposed with regression lines, were studied to determine the model form. The quadratic growth model appeared by visual inspection to be a good fit. This fit was confirmed by comparing the Bayesian Information Criterion of the linear and quadratic models (see the appendix for more details). To examine the variance components, an unconditional means model with no predictors at either level was estimated. In addition, the results of an unconditional growth model with time as a level 1 predictor and no predictors at level 2 were used to evaluate baseline change. Because the effect of group and intervention and the group × intervention interactions were of interest, those predictors were added in level 2. The terms and parameters are as described in the generic quadratic model above and in the appendix. An examination of each residual separately with normal probability plots ensured that the error structure was independent and normally distributed with a mean of 0 and constant variance, σ² and that the growth parameters did in fact vary across subjects.

The model can be used to estimate interesting aspects of the experiment. For example, we can extrapolate when the postpartum groups are expected to have similar thresholds. This could be done using simulation methods from the model (i.e., assuming random variation as observed from the sample), or simply using a calculation of a point estimate (i.e., solving the model equations for no difference between the two groups). The latter approach is made more complicated when large individual differences in the fitted parameters are observed such that any point estimate is not representative of any individual in the sample. As an illustration, the model created in this study predicts the two groups to have the same withdrawal threshold on day 77. In addition, a formal statistical inference could be generated from a model of this type, but this would require an equivalence judgment (i.e., that the two groups do not likely differ by more than some critical threshold) that was not planned for here. Finally, extrapolations outside the model space are always prone to error (e.g., a poorly specified model may fit the data well within a certain range, but lead to large errors on extrapolation).

Independent samples Kruskal–Wallis testing was used to examine withdrawal threshold before surgery. Where appropriate, all hypothesis testing is two-tailed with a statistical significance threshold of P value less than 0.05.

Statistical analysis for the growth curve model was performed with SAS version 9.2 (SAS Institute, Inc., Cary, NC). Full SAS code is available upon request. PROC MIXED uses restricted maximum likelihood estimation to analyze multilevel models, and can incorporate both random and fixed effects. To fit an individual subject growth model, individual intercepts and time (slopes) are allowed to vary randomly with the remaining terms in the model treated as fixed effects. Because of the small sample size, the quadratic term was also estimated as fixed. The primary inferences of the study involve examining if experimental condition impacts some aspect of the change process. To estimate these effects, group and intervention (Int = weaning at 21 days postsurgery) were entered as level-2 predictors of intercepts, slopes, and quadratic change parameters (i.e., as group × parameter interaction). Because subject intercepts and slopes do not have the same variance because of the introduction of heteroscedasticity from allowing the slopes to vary by subject, an unstructured variance/covariance matrix has been specified in the random statement, but other covariance structures were also considered.⁴ 

For comparison with the growth curve modeling, the data were also examined using two-way repeated-measures ANOVA using SPSS version 22 (IBM, New York, NY). To minimize loss of data due to listwise deletion, measurements were assigned to nearest neighbor days where possible. To meet the assumptions of two-way repeated-measures ANOVA, the positively skewed response variable was logarithmically transformed. After transformation, there were no outliers and the distribution of the data satisfied parametric assumptions for each group.

To illustrate the efficiency of the growth curve model in comparison with the two-way repeated-measures ANOVA model, statistical power analyses were performed on data simulated with R version 3.0.2 (R Foundation for Statistical Computing, Vienna, Austria) and RStudio version 0.98.501 (RStudio, Inc., Boston, MA). Using MASS package (version 7.3–29; Oxford, United Kingdom), data sampling was done from a multivariate normal distribution so that the simulated data possessed the same means and SDs (within 0.005) by time point and similar correlations between time points as the actual observed data of the postpartum-with-pups and postpartum-separated-from-pups groups only. Five thousand datasets (i.e., runs) were simulated for each sample size per group of n = 6, 9, 12, 15, 18, 24, and 36, for a total of 35,000 datasets. Repeated-measures ANOVA and growth curve analyses were performed on each individually simulated untransformed dataset as was done with the actual observed experimental data.

In the repeated-measures ANOVA power simulations, the proportion of runs with statistically significant differences (P < 0.05) found in the group × time interactions was interpreted as the statistical power for each sample size. This was the effect of interest given that group differences over time were the primary focus of the analysis. In the growth curve power analyses, a smaller model with parameters of intercept, time and time2 was considered nested within a larger model with additional parameters of group, group × time and group × time2 (i.e., group impacts multiple aspects of change). In both the larger and smaller growth curve models, intercepts and slopes were allowed to vary between subjects. The chi-square statistic for the likelihood ratio test was obtained for each model by run and sample size. These 35,000 (5,000 per sample size) runs of both models were compared using chi-square difference tests that examine the difference in model fit between the larger (with group) and smaller (without group) models in the context of the number of parameters in each model.^13,14  The proportion of statistically significant differences (P < 0.05) found between the two models were interpreted as the statistical power for each sample size.

All the animals recovered from surgery without evidence of infection, and all deliveries occurred spontaneously with an average litter size of 11. Median (interquartile range) withdrawal thresholds before surgery did not differ among groups (males 29 [20 to 36] g, virgin females 22 [17 to 26] g, postpartum-with-pups 24 [19 to 26] g, postpartum-separated-from-pups 23 [18 to 26] g, P = 0.22 by Independent Samples Kruskal–Wallis Test), nor did they meaningfully differ on the lower bounds of the 95% CIs estimated from the model (13.2, 11.2, 13.4, and 14 g, respectively). A gradual increase in withdrawal threshold over 70 days after surgery, as previously noted in males¹¹  was present in all groups. In addition, a transient decrease in withdrawal threshold occurred after weaning on day 21 in the postpartum-with-pups group as previously described.² 

Figure 1 shows the scatterplot of actual withdrawal threshold measurements and the predicted trajectory fits of two representative subjects, one from the postpartum-with-pups group and the other from the virgin females group. The predicted group trajectories for all four groups are displayed in figure 2, where the postpartum are shown in figure 2A and the nonpostpartum in figure 2B. Individual subject modeled trajectories detailing missing data by breaks in the line plots and demonstrating individual variability by group are shown in the “spaghetti plots” in figure 3. The model coefficients predicting withdrawal threshold are reported in the appendix.

The virgin females group exhibited higher thresholds after surgery than the other groups (see intercepts in fig. 2, A and B). The postpartum-with-pups (pp w/pups) group exhibited an increased trajectory before removal of pups compared with the postpartum-separated-from-pups (pp w/o pups) group (fig. 2A). The postpartum groups had nonoverlapping 95% CIs of modeled trajectories between days 6 and 66 after surgery demonstrating difference during the recovery period (fig. 2A). Immediately after surgery, the trajectory of the males group was not different from the virgin females (fig. 2B) and postpartum-separated-from-pups groups, but was different from the postpartum-with-pups group with nonoverlapping 95% CIs of modeled trajectories between days 3 and 64 postsurgery (not shown). In addition, substantial individual differences were observed within groups as evidenced by the varying patterns of change within each group in figure 3.

Regarding the random effects, the variance component (standard error) for intercepts (2.60 [0.96], P = 0.004) and slopes (0.002 [0.0007], P = 0.005), were statistically significantly different than 0, indicating that substantial individual differences were observed for these parameters across individuals (fig. 3). There was no association between intercepts and slopes, (−0.02 [0.02], P = 0.381), indicating that initial withdrawal thresholds are not good predictors of linear trajectories. Finally, the residual variance estimate (8.43 [0.53], P < 0.001) remained significantly different from 0, suggesting additional variation in intercepts and slopes not explained by the current factors and interactions (fig. 4). This leads to the assumption that there are additional factors that could explain the variation not considered in the present model.

Main effects for group (P < 0.001) and time (P < 0.001) were both observed in the ANOVA model, a statistically significant group × time interaction was also observed, F(30,370) = 3.768, P < 0.001. This indicates that group differences change over time. This interpretation is similar to what was observed in the growth curve modeling. To visualize the intersubject variability within group, figure 5A shows a scatterplot of actual withdrawal threshold measurements for the postpartum-with-pups group. This within subject variability not accounted for by group and time is attributed to error in the two-way repeated-measures ANOVA model, whereas this variability is assigned to the distribution of random effects in the growth curve model.⁴  To illustrate the difference in the nature of the estimates of each model and their corresponding precision, an example plot (fig. 5B) of the postpartum-with-pups group was analyzed using both approaches. Close inspection reveals narrower 95% CI widths for the growth curve model which demonstrates greater precision (i.e., the point estimates of withdrawal threshold are associated with less uncertainty). This greater precision would translate into greater statistical power for the growth curve approach. Power analyses conducted on the simulated data, collected as set forth in the Materials and Methods, show that 80% statistical power to detect group difference is obtained at 6 per group in growth curve modeling, whereas 14 per group are required to obtain 80% power in the repeated-measures ANOVA modeling (fig. 6). However, the two simulations estimate two very different parameters from the approaches and are not simply different estimates of the same construct. For instance, in the ANOVA simulations, the group × time interaction effect was examined, whereas a model fit was the focus of the growth curve simulations.

These results carry both biologic and methodologic implications for the study of recovery from peripheral injury. As regards biology, we confirm previous observations that peripheral injury during the postpartum period results in a lesser degree of persistent hypersensitivity than in virgin females, that this effect requires the presence of pups, and that there is a transient worsening of hypersensitivity when pups are weaned.²  We also confirm that recovery from hypersensitivity after pSNL injury occurs over several weeks, in distinction to spinal nerve ligation, which is permanent in most animals with partial recovery over many months in a minority of subjects.¹⁵ 

We extend previous biologic work in two ways. First, we show that recovery from hypersensitivity occurs more rapidly after pSNL when this injury occurs in the postpartum period than in virgin females. The similar degree of initial hypersensitivity and clear divergence of trajectory within days of injury (fig. 2A) strongly suggests an active process in the postpartum period hastening recovery. These data agree with observational data in women after cesarean delivery, and a very low prevalence (<1%) of surgery-induced pain 1 yr after surgery.¹  Second, we demonstrate a lack of sex difference in the trajectory of recovery from hypersensitivity after pSNL. Although this result superficially disagrees with clinical observations that chronic pain after surgery is more common in women than in men,¹⁶  this difference is small (10 to 30% increased risk in women) and could not be tested with adequate power in traditionally small laboratory studies in animals.

Withdrawal threshold decreased after removal of pups 21 days after delivery, the time of normal weaning. The cause of this transient hypersensitivity, which also occurs in the absence of nerve injury,²  may reflect maternal stress at this time, because stress itself is associated with hypersensitivity in rats.¹⁷  Alternatively, hypersensitivity at this time may reflect an acute reduction in oxytocin release with cessation of lactation, wince withdrawal threshold is temporarily reduced in postpartum rats with pups by intrathecal injection of an oxytocin receptor antagonist.² 

The current study also illustrates the potential value of growth curve modeling to laboratory studies of recovery after injury and provides much of the information needed for researchers to apply this tool. Strengths of individual growth curve modeling relative to more traditional methods of analysis include flexible specification of how the outcome responds over time (i.e., the model form), and explicit modeling of both group-level and subject-level growth curves. In this approach, change is modeled at the individual level using change parameters, allowing for the examination of individual variability in intercepts and rates of change in the outcome under investigation (random effects). The ANOVA approach focuses on change at the aggregate level and assumes variation around a fixed effect for each group. In addition, the growth curve model has the ability to accommodate embedded missing data, whereas listwise deletion is used in ANOVA. Due to its flexibility, growth curve modeling easily handles nonnormally distributed outcomes (e.g., binary, rates, counts), where an ANOVA model requires assumptions about the distribution under study.

Growth curve modeling is not without its weaknesses and is not advocated for all laboratory studies of change. Model misspecification can become an issue as there are many choices of model forms that can be used. For example, is a quadratic or a (log)linear model the correct choice? Models are evaluated to ensure they describe the data well as set forth in this article, but all models are wrong to varying degrees.¹⁸  The number of repeated measures is an important consideration, with small numbers of measurement occasions (e.g., t <4) greatly limiting the value of the approach. Furthermore, the estimation of complex random effects can be problematic in small-n studies. In the current study, because of the estimation difficulties, we opted to only allow intercepts and slopes to vary randomly across subjects. In addition, because there is a higher level of technical sophistication needed, interpretation can be difficult for persons familiar with the ANOVA approach, although this can be overcome with plots and practice. Finally, at the time of this writing, there are no simple power calculations for most applications. Published tables do exist, but statistical simulations are generally helpful.¹³ 

In summary, using a model of surgical nerve injury with recovery of hypersensitivity over a couple of months, we demonstrate a similar time course of recovery in males and females, but an accelerated rate of recovery in females when surgery occurs in the immediate postpartum period. These data agree with a high incidence of acute pain, but a low incidence of chronic pain in women undergoing traumatic or surgical delivery, as observed clinically. We describe the advantages of growth curve modeling to better understand and probe individual-level and group-level differences in the time course of recovery after injury in animal studies. Because recovery is modeled in each animal, this approach provides a unique tool to study the broader question of promoting recovery from injury rather than the dichotomous question of presence of pain/hypersensitivity at an arbitrarily defined “chronic” time after injury.

Supported in part by grant R37-GM48085 from the National Institutes of Health, Bethesda, Maryland.

Dr. Eisenach has received payments for consultation from Adynxx (San Francisco, California) and Aerial Biopharma (Morrisville, North Carolina) on topics unrelated to this article. Dr. Houle has received payments for consultation and research from Merck, Sharp & Dohme (White House Station, New Jersey) on topics unrelated to this article. The other authors declare no competing interests.

In a growth curve model with differing groups and interventions, the hierarchical structure where lower-level units are nested in higher-level units consists of repeated measurements/individual subject growth trajectories (level 1) and variation in between-subject growth (level 2). Covariates that account for within-person changes that occur across measurements are incorporated into level 1 (time-varying covariates). The level 2 equations model each of the level 1 parameters. Time-invariate covariates are added in level 2. These are the predictors that do not change over time, such as group and intervention.

The first step to building the model is to plot and study the data using a line plot with the time variable on the x-axis and the dependent variable on the y-axis for each subject and superimposing a regression line. Investigate whether the growth is linear, quadratic, or cubic and determine the best fit. Many times the linear, quadratic, and cubic models will be tested by the likelihood ratio test (chi-square difference test) or by comparing the Bayesian Information Criterion of the models to determine whether one is a better fit than another.¹⁹ 

For the following model build, the equation structure follows that of Raudenbush and Bryk.²⁰ Y_ti is the observed measurement at time t for subject i, i = 1,…,n subjects, where a_ti is time t for person i and π_pi is the growth trajectory parameter p for subject i associated with the polynomial of degree P. Therefore, π_0i is the intercept, π_1i is the coefficient for linear slope (i’s true rate of linear change per unit of time) and π_2i is the coefficient for quadratic slope (i’s true rate of quadratic change per unit of time which can also be understood as acceleration or deceleration). The portion of subject i’s outcome that is unexplained by the specified growth parameters on occasion t is e_ti. The parameters for level 2 are set forth in table 1. The level 2 residuals, which are deviations of individual change trajectories around the predicted averages, are denoted r_0i, r_1i, and r_2i.

The next step is to run two unconditional models which aid in examination of variance. The first is an unconditional means model with no predictors at either level. This model partitions the outcome variation into within-subjects (variation within subjects over time) and between-subject (variation between subjects with no time effect).

From this model, we determine whether the outcome is worth exploring further by examining where the systematic variation lies.

The second is an unconditional growth model with time as the level 1 predictor and no predictors at level 2. This model allows evaluation of baseline change.

This random slopes and intercepts model shows how much total variation there is within and between-subjects. This information aids in the decision of whether to add level 2 predictors.

If there is significant variation in initial status and rate of change, exploration of time-invariate predictors to account for variation in growth parameters across individuals is warranted. In a growth curve model where the effect of group and interventio n are of interest, group and intervention and their interaction can be chosen as predictors. Finally, the quadratic model with these level 2 predictors is:

When substituting level 2 effects into level 1, a composite growth model is formed:

This formula shows how the outcome depends concurrently on the level 1 and level 2 predictors and the cross-level interactions with time. The final step is to ensure that the error structure is independent and normally distributed with a mean of 0 and constant variance, σ² and that the growth parameters vary across subjects. Nonnormality of the level 1 error, e_ti, will bias the standard errors at both levels one and two. An examination of each residual separately (e_ti, r_0i, r_1i, r_2i) with normal probability plots is a method commonly used to check the error structure. If the residuals are not independently and normally distributed, a data transformation procedure can be used. These procedures are set forth in detail in the study by Judd and McClelland.²¹  The estimation method for the covariance parameters is restricted maximum likelihood, which is a form of maximum likelihood estimation which does not base estimates on the fit of all of the data, but instead uses a likelihood function calculated from a transformed set of data, so that nuisance parameters have no effect. The maximum likelihood method selects values of the parameters that maximize the agreement between the model and the observed data. In the residual maximum likelihood method, the likelihood function is partitioned into components representing information about the parameters of interest and other nuisance parameters and then completes the estimation.

A quadratic model has differing rates of change at each time point. These rates of change are estimated by the slopes of the tangent lines to the growth curve at each point. These are called simple-slopes and they change across the span of the growth curve.¹⁹  The rates of change are the first derivatives with respect to time of the level 1 equation evaluated at each particular time point. In the composite specification above, the rate of change is represented by π_1i + 2π_2ia_ti, where π_1i is the instantaneous rate of growth at initial status.

Growth curve modeling allows for incomplete or missing data if the data are missing at random. But, the models do require one time point more than the number of growth parameters in the level 1 model. A quadratic model with three growth parameters in the level 1 equation, π_0i, π_1i, π_2i, requires at least four time points. The more time points that are obtained, the more precise the model.¹⁹  Before analysis, data must be organized in long form where there are multiple rows per subject, one for each time point at which the subject has a measurement. The number of rows equals the number of measurements taken. The variables may include subject ID, group, time, intervention, and outcome. The time variable will need to be created and represents the passing of time for which measurements were taken. The first row for each subject should be coded as the initial status with Time = 0. Therefore, the duration that is subtracted from the first time point to start the new variable “Time” at 0 should be subtracted from each time point thereafter. For example, if the first measurement was taken on day 3, the Time variable for the first measurement is 3 – 3 = 0. The Time variable for day 4 is 4 – 3 = 1, etc. This process of coding is called “centering” and allows for easier interpretation of the estimated parameters. If working with a quadratic model, a variable for time squared should also be created. These data conditioning procedures will be set forth in more detail in the data obtained in the current study.

Following the steps set forth for the building of a generic growth curve model outlined above, plots for each subject with the Time variable on the x-axis and withdrawal threshold on the y-axis, superimposed with regression lines, were studied to determine the model form. The quadratic model appeared by visual inspection to be a good fit. This fit was confirmed by comparing the Bayesian Information Criterion of the linear and quadratic models. To examine the variance, an unconditional means model with no predictors at either level was run. In addition, the results of an unconditional growth model with time as a level 1 predictor and no predictors at level 2 were used to evaluate baseline change. Because the effect of group and intervention and the group × intervention interactions were of interest, those predictors were added in level 2. The final quadratic model to predict behavioral recovery, that is, withdrawal threshold over time, is:

The terms and parameters are as described in the generic quadratic model above and in table 1. The model coefficients predicting withdrawal threshold in our small animal study are reported in table 2.

Y_ti is the predicted withdrawal threshold measurement for animal i at time t, a_ti is time t for person i. Intervention is a dummy variable (0/1), 0 for up to and including day 21 and 1 for after day 21. Group is postpartum-with-pups, postpartum-separated-from-pups, males or virgins. An examination of each residual separately with normal probability plots was done to ensure the error structure was independent and normally distributed with a mean of 0 and constant variance, σ² and that the growth parameters vary across subjects.

To condition data, we set the intervention dummy variable “Int” to either 0 or 1. If day is less than or equal to 21, Int = 0. If day is greater than 21, Int = 1. As we are modeling recovery, we removed the presurgery baseline measurement and let Time = Time − 1 so that the initial status is Time = 0 [Postoperative day = 1]. Create a new variable “Time squared” which is computed Time × Time (table 3). In addition, data were allowed to be unbalanced.

Statistical analysis for this study was performed with SAS version 9.2. Full code is available upon request. PROC MIXED uses restricted maximum likelihood estimation as discussed above to analyze multilevel models and can incorporate both random and fixed effects. To fit an individual subject growth model, intercepts and time (slopes) are allowed to be random as denoted in the random statement, with the remaining terms on the model statement being fixed. Because of the small sample size, the quadratic term, Time squared was not allowed to be a random effect. The subject = subject option of the random statement specifies that the random terms should vary across subjects. Because subject intercepts and slopes do not have the same variance due to the introduction of heteroscedasticity from allowing the slopes to vary by subject, an unstructured variance/covariance matrix has been specified in the random statement, but other covariance structures should be considered.⁴