A maximum Likelihood Approach to Analyzing Incomplete Longitudinal Data in Mammary Tumor Development Experiments with Mice

Longitudinal mammary tumor development studies using mice as experimental units are affected by i) missing data towards the end of the study by natural death or euthanasia, and ii) the presence of censored data caused by the detection limits of instrumental sensitivity. To accommodate these characteristics, we investigate a test to carry out K-group comparisons based on maximum likelihood methodology. We derive a relevant likelihood ratio test based on general distributions, investigate its properties of based on theoretical propositions, and evaluate the performance of the test via a simulation study. We apply the results to data extracted from a study designed to investigate the development of breast cancer in mice.


Introduction
The objective of this article is to develop an inference method based on outcomes from mammary tumor studies containing different types of missing observations.We derive accurate likelihood functions and focus on the maximum likelihood comparison between populations that may have values below the limit of detection caused by low instrument sensitivity.Specifically, the method proposed in this article is motivated by a gene deficiency study performed at Roswell Park Cancer Institute, Buffalo, NY, which investigated the longitudinal growth of breast cancer in mice.The mouse tumor model was used to study the relationship between the loss of the gene Thoc1 and mammary tumor growth progression (Wang, 2008).In this mouse model, mammary hyperplasia can be detected as early as 4 weeks of age.In addition a large percentage of mice develop carcinoma by 14 weeks (Kelly- Spratt et al., 2008).Female mice with three genotypes, namely Thoc1F/FPyMT+/-, Thoc1F/FMMTVCre+/-PyMT+/-and Thoc1F/+MMTVCre+/-PyMT+/-(control) were followed for 35 weeks (Wang, 2008).The mice were palpated once weekly from 8 weeks of age to monitor mammary tumor development.
Breast tumors can be palpable when their diameters reach 0.5mm (Kuang et al., 2004); however, the tumor is often too small to be detected during certain periods (e.g., early stage) of its development (e.g., Tan et al., 2002).This problem can be considered in terms of the limit of detection (henceforth, referred as LOD) due to instrument sensitivity (e.g., Schisterman, 2006;Perkins et al., 2007;Vexler et al., 2008;Yu et al., 2009).Further, when primary tumors reach a certain size (approximately 2cm in diameter), mice are euthanized according to the study protocol, thus generating missing data.In general, a substantial proportion of data is missing as a function of time, leading to underestimated tumor growth rates.Another possible cause of missing data is a drop-out from natural death attributed to a treatment effect (e.g., toxicity, see Tan et al., 2002).In this setting, the missing pattern itself has valuable information related to treatment response.
When portions of the data are below the LOD, existing literature commonly suggests a parametric approach (e.g., Lynn, 2001;Schisterman, 2006;Perkins et al., 2006;Vexler, 2008), which we adapt in this article.The current literature recommendations are slanted towards the use of parametric distributions to accommodate the unobserved components of data distributions.The approaches based on normal assumptions have been frequently developed (e.g., Laird and Ware 1982;Hughes, 1999;Jacqmin-Gadda et al., 2000;Tan et al., 2002).To handle missing data in longitudinal studies, parametric methods also have been actively developed (e.g., Rubin, 1977Rubin, , 1987;;Rosenbaum and Rubin, 1983;Little, 1988;Diggle and Kenward, 1994).Hughes (1999) and Vaida et al. (2007) proposed an MCEM algorithm and hybrid EM algorithm (lmec package in R), respectively, to handle censored responses in the frame work of a linear mixedeffects model described in Laird and Ware (1982).
For mouse tumor models, most statistical studies have focused on measuring the tumor.Koziol et al. (1981) proposed a nonparametric method based on multivariate rank statistics to compare mice tumor growth curves, but did not handle the critical LOD issue.Tan et al. (2002) considered this in the tumor xenograft model by assuming normal distribution based models proposing two alternative parametric approaches, one based on EM algorithm and the other on a Bayesian methodology.Missing values were considered as missing at random while the death rate was not of specific interest in their study.Liang and Sha (2004) used a non-linear mixed effect model to analyze the data from Tan et al. (2002), where the observations were assumed to have normal distributions and the dropout pattern was considered as missing at random.The death rate was not taken into account in their analysis and the unavailable data caused by the LOD were simply imputed by a single value in their analysis.Albert and Shih (2003) modeled the relationship between clones of mouse cells with genetic mutations and tumor progression.Modeling tumor growth was suggested as an approach to overcome the bias problem caused by the missing data, where the drop-out mechanism was considered as missing at random.Mean tumor growth was defined as a function of the probability of the onset time in conjunction with a model of tumor growth after onset.The LOD was not of specific interest in their study.
In this article, we develop a test to compare K-groups utilizing the maximum likelihood ratio which is known to be locally most powerful taking into account the LOD and the drop-out mechanism jointly.Unlike other tumor studies with animals, for the mammary tumor experiments, the tumor is measurable throughout the study; however, the tumor measurement is often subject to the limit of detection due to insensitive tools as stated earlier.Besides, under a wellcontrolled study environment, the cause of drop-out is clear, thus often the dropout itself has information to be considered.We classify the causes of missing data into two categories, namely missing data due to euthanasia and missing data due to natural death (causes other than euthanasia).This distinction is important in that natural death provides extra information regarding treatment differences, while the decision of euthanasia directly depends on the size of the tumor.While we adapt the likelihood approach, in principle, we propose a method that is not limited by the normal distribution assumption but can work with a data set under more general distribution assumption (e.g., Bondon, 2009), so that our method has its own unique generalizability (more details in Section 4).This article has the following structure.In Section 2 we first specify the observation patterns for successive assessments of tumor volume, and derive the parametric likelihood function based on incomplete longitudinal data.A likelihood ratio test for examining differences among K -groups is developed and its asymptotic properties are considered.In Section 3 the performance of the proposed method is evaluated using an extensive Monte-Carlo study.In Section 4 we apply the proposed method to analyze actual data sets from mouse tumor experiments.Section 5 is devoted to concluding remarks.

Model development
Consider a study with K treatment groups.For treatment group k ( there are k n experimental units.Each experimental unit is assumed to be examined at J equally spaced time points by a pre-specified schedule.Let ijk X denote the j th measurement of the tumor size of the i th unit in the k th group ( where   1    reflects the association between two consecutive measurements on the same experimental unit and ijk  are independent identically distributed random variables with the mean 0 and variance 2  .Model (1) is based on an autoregressive relationship between measurements.For example, in the particular case with One can remark that it is a very complicated problem to test for the correlation structure subject to the pattern of the missing data.In the context of the mammary tumor study, model ( 1) is reasonably supported by the sample correlations based on the observed data as shown as follows.For our logtransformed data, we were able to estimate the sample correlation between weeks 17 and 18 at 0.896 and that between weeks 17 and 19 at 0.797.Likewise, the sample correlation between weeks 19 and 20 was 0.930 and the correlation between weeks 18 and 20 was estimated at 0.866.Note that the selected example weeks are less affected by the LOD and death of animal mechanisms.
We note that higher order autoregression structures can be used by modifying (1), where the autoregression coefficient  can be considered as a function of j, which will require substantial mathematical details and be subject to future study.Note that, using (1), the variables   , where   (1 ) / (1 ) , indicating more variability into the variance/covariance structure when the time goes on.This covariance also indicates that, when the time passes, the values of the variance are more stabilized as j increases, much alike the typical autoregressive covariance structure with normally distributed data (e.g., Diggle, 1988;Tan et al., 2002).
Suppose that a mouse i in group k is euthanized if its tumor size is greater than Table 1 summarizes the possible cases of two consecutive observations.There are a total of 7 different scenarios.The values subject to the LOD are indicated by N/A and the check mark (✓) indicates that the tumor is detectable since . The dot indicates a missing value due to death of a mouse.Cases 5 and 6 both indicate missing at time j by reasons other than euthanasia.Case 7 indicates that an animal is euthanized at time j.

Defining the likelihood function
In this section, we limit our discussion within a single group (say, the k th group) for simplicity.In this case, the general form of the likelihood function is where θ indicates the parameters to define model (1) and f accounts for the appropriate contribution of ij S to the likelihood function.Let     where I is an indicator function that takes the value of 1 if the condition in the bracket is satisfied and 0 otherwise.According to the cases in Table 1, (4) can be expressed as , where 5) indicates the contribution from Cases 5-7.
Consider the conditional probabilities in (5).While euthanasia is directly a result of tumor growth, death by a natural cause is independent of the measurement process.Let 0 ij A  indicate that animal i died by a natural cause at time j , and . .
Similarly, for Cases 3 and 4, one can show that the probabilities of the measurement process ij X and ij A are separated.For Case 5, Similarly, for Case 6, .
Incorporating Formulas ( 6) and ( 7) into the likelihood function allows us to take into account the probabilities of natural deaths in hypothesis testing.
Differentiating (10) with respect to , 1 i j x  and ij x , we obtain since the denominator of ( 9) is   For Case 3 where , 1 i j X  is observed and , Note that, in a different context, a similar approach to these modifications of conditional probabilities is found in Park et al. (2007), where items such as j  and 2 j  are ignored.Note that, related to mammary tumor experiments, the values of j  and 2 j  may be relatively large since the observations are likely to be highly correlated.The density of   , if , and , if Combining ( 6)-( 8), and ( 11)-( 15), we derive the likelihood function with respect to a sequence of measurements for each animal.

The test statistic
be the conditional density defined within group k ( be the conditional density for all groups; that is, . The null and alternative hypotheses to test are : not where  for group k , and j p be jk p under the null hypothesis at time j .Under 0 H , the likelihood function has the form where   , 1, ijk i j k f x x  are obtained as described in Section 2.2 for different cases, , 1 . Finally, given the formal notation at (17) and ( 18), we complete the likelihood ratio statistic, to carry forth testing of hypothesis 0 H in (16).

Implementation to analyze the real data set
The estimation of j  for time j ( 1,..., ) may be biased because of the LOD mechanism and the large portion of missing values due to euthanasia as a consequence of large tumor sizes.Modeling tumor growth dependent upon time j should be a reasonable alternative to the current approaches.In general, the growth of tumors can be approximated by polynomial forms.For example, (20) assuming that the mean as a function of time has a   .Note that if an anti-cancer treatment is effective (e.g., Section 3), higher order polynomials may be desirable since the tumor size often decreases before it increases again.We also consider that jk p is constant over time, that is, 1 . The reasons for this consideration are two-fold.First, the rate of natural death incurred by toxicity can assumed to be constant throughout the study period.Second, estimating a death rate for each time j may not be feasible with small sample sizes.The robustness of this approach for an assumption of changing rates is discussed in Section 4. Using this adaptation and the likelihood function (17), the maximum likelihood estimator of the death rate, p ˆ, under 0 H , is The derivation of this estimator is easily carried out since the function of the measurement process and the function of the death rate in the likelihood ( 17) can be handled separately.Under 1 H , ˆk p is obtained for each group, k , in a similar manner to that shown above.By performing an appropriate transformation (e.g., Box-Cox power transformation), and assuming the tumor volumes to be normally distributed, the exact formulas of ( 11)-( 14) based on the normal distribution are provided in the Appendix A. In addition, for interested readers, we also outline an approach based on general distributions in the Appendix B.
Employing a linear trend for tumor growth and the same death rate throughout the study, θ in the likelihood ratio test (19) consists of , , ,     and p .With general underlying distributions, we can show the following proposition pertaining to the likelihood ratio test under some regularity conditions (Silvey, 1970).

Proposition 1
Under 0 H , (19) converges to 2  distribution with 3( 1) K  degrees of freedom when 1, , The degrees of the freedom reflect the difference between the number of parameters under 0 H (i.e., , , p   ) and the number of parameters under 1 H (i.e., different , , p   for each group).Since the development of ( 19) is based on the parametric likelihood as shown previously, the proof of this theorem is an immediate consequence of Wilk's Theorem pertaining to the likelihood ratio test.In Section 4, we evaluate the performance of the test statistic based on the extensive Monte Carlo simulations.Let θ indicate the maximum likelihood estimator of , , ,     and p Let N be . The following proposition is also immediate given the asymptotic normality of the maximum likelihood estimator under some regularity conditions (Silvey, 1970).Since the expectation is difficult to calculate for Fisher's information matrix, we can use the observed Fisher information evaluated at the maximum likelihood estimate (Park et al. 2007;Efron and Hinkley, 1978).We carried out the maximization of the likelihood function and the accompanying parameter estimations using an optimization function in R (optim, http://www.rproject.org).Optimization was carried out using the 'Nelder-Mead' algorithm and the likelihoods were calculated under 0 H and 1 H , respectively.Initial values of  ,  ,  and  were obtained by least square approach (e.g., the mixed effect model with the covariance structure of autoregression with order 1), based on the observed data (i.e., no consideration of the LOD or missing data).
Commonly, values of the -Log likelihood function around these initial values gradually decreasing toward an optimal point, thus the Nelder-Mead method provides a reasonable performance.The convergence is reasonably fast depending on a machine and the sample size (e.g., one optimization with 61 experimental units and 21 weeks takes approximately one minute using the IBM personal computer with Pentium D processor [3 GHz of CPU, 2.00 of RAM]).
The program is available from the authors on request.

Simulation study
To investigate the performance of the proposed method, a simulation study was conducted.We use the characteristics of generated samples that are empirically close to those of the mammary tumor data from the Thoc1 gene deficiency study.
For our simulation studies, three treatment groups ( 3 K  ) were investigated.We first carried out 5000 Monte-Carlo simulations based on a "naïve" practice in which the repeated measures analysis was performed based on a linear model approach via a generalized least squares method.A first order autoregressive process was also assumed for the correlation structure.The likelihood function was evaluated based on the generalized linear model with and without the group factor to construct the likelihood ratio test.The interaction between time and group was considered in the model along with the group factor.In the simulations, we supposed that the data were log-transformed prior to the analysis.Note that in the real data set, N/A cases due to the LOD were simply recorded as 0. In a naïve manner, unavailable observations due to the LOD were replaced with log(0.00001).We used intercept and slope values around -6 and 0.45, respectively, to emulate the log-transformed tumor volume growth in the mammary tumor study.The number of visits was 21 (J=21) and the total sample size was 63 with three balanced treatment groups.Different values of  and  in equation ( 1) were considered to evaluate the effect of the different magnitudes of them (e.g., 0.3   , 0.9   , or 0.9 ).The values of L d and ik R d were log(0.005)and 3, respectively.The natural death rate for each visit was 0.01.The target Type I error and 95 th percentile based on the 2  distribution ( df = 4) were 0.05 and 9.4877, respectively.Table 2 shows the Type I error and the 95 th percentile of the distribution for the test based on the naïve approach based on the Monte-Carlo simulations.The results in Table 2 show poor approximations to the 2  distribution.Among the examined cases, the Type I errors ranged from 0.0004 to 0.0870 and the 95 th percentiles ranged from 3.6480 to 11.6160, which indicate that the naïve test is not independent of the parameters.We also increased the sample size ( 90 N  ), but this did not improve the approximation (Table 2).It is probably because the LOD problem is in effect (e.g., Perkins, 2007Perkins, , 2009)).This result clearly demonstrates that the naïve approach may not be appropriate for analyzing tumor growth data with unfavorable characteristics.The method proposed in this article could be used for these kinds of data sets.
We next examined the Type I error and power for the proposed method given the various scenarios.Similarly to the scenario described above, the simulations were carried out in order to emulate the mammary tumor development data.
First, the values of log-transformed tumor volumes were generated based on the linear model with an intercept of -6 and a slope of 0.45 for all three treatment groups.The values of ij  were generated from Normal (0, 0.3 2 ), and  in (1)   was given 0.9.There were a total of 21 observation times.The natural death rate was 0.01 per week.The values of L d and ik R d were -log (0.005) and 3, respectively.The death rate for each time was 0.01.The sample was equally allocated to all three treatment groups.Toward the end of the study, the survival rate was approximately 40% on average.The simulated Type I errors for these scenarios are shown in Table 3.The survival rate in the table is the rate of mice alive at time J .The survival rate shows that data become scarce toward the end.The target Type I error is 0.05.The simulated Type I errors are somewhat greater than that when the sample size is small, but improve when the sample size increases.Table 3 also presents some other scenarios with different combinations of the sample size,  ,  and the number of visits.The simulated powers for various models. ,  and p are the intercept, slope for the linear function of j  and the death rate, respectively.Among three treatment groups, one group (indicated by the subscript 0 for the parameters) is different from the other two groups (indicated by the subscript 1 for the parameters). and  are 0.3 and 0.9, respectively.N indicates the total sample size.The survival rate is the rate of survived animals at time We also investigated the robustness of the model with constant death rate under the assumption of non-constant death rate.We generated data in which the death rates satisfied the model . To be comparable to the entries in Tables 3 and 4, we set up the overall mean death rate to be approximately equal to the constant death rate from Tables 3 and 4 (e.g.,   01 .0 / J p j ).
The proposed method demonstrated an overall robustness in terms of controlling Type I errors.For example, the Monte Carlo Type I error corresponding to the third entry in Table 3 ( 30 n  ) was 0.07 which is comparable to the Type I error for the correct model with a constant death rate in Table 3, namely, 0.065.The powers with the same parameters and similar death rates to the second ( 30 n  ) and the fifth ( 63 n  ) entries in Table 4 were 0.08 and 0.35 indicating that the power is somewhat compromised.Note that the power of the proposed method is not compared with any other approach (e.g., Koziol et al., 1981;Tan et al., 2002) as such tests do not use a fully correct model specification.
To demonstrate a generalizability of our method, i.e., an extension toward the other distributions, we also carried out the simulations based on the tumor measurements only using the error terms with the Laplace distribution.The simulations are based on ( , , , , , ) J K      ( 6, 2.5, 3, 0.9, 5, 2)


. The method to obtain the conditional densities is outlined in the Appendix B. Expectedly, the simulation results (based on 1000 simulations) showed a very good operating characteristic (Monte-Carlo Type I error = 0.050) even with a relatively small sample size ( 1 2 10 n n   ).Note that the likelihood function approach using the normal distribution showed a less accurate but viable performance (Monte-Carlo Type I error=0.055).

Applications
The Roswell Park data set includes observations for 21 weeks, from week 10 to week 30.A total of 63 mice were used (21, 24, and 18 mice for the three groups).No tumor less than 0.01cm 3 was detected in the actual data set.At week 10, 89% of the data were below the LOD even though it is known that mammary hyperplasia can be detected as early as 4 weeks of age for the transgenic polyoma middle T (PyMT) mouse model.Toward the end of the study (week 30), 89% of mice were euthanized so that the data were markedly scarce.
Table 5 shows the numbers of observations that are below the LOD and missing toward the end of the study.Note that none of the mice died naturally.Nevertheless, our proposed method and the data analysis strategy can easily be adapted to these.Tumor volumes are log-transformed before comparison.Histograms for the control group (Thoc1F/+MMTVCre+/-PyMT+/-) compared with the other two genotypes (Thoc1F/FPyMT+/-and Thoc1F/FMMTVCre+/-PyMT+/-) are shown in Figure 1.In the early period (weeks 13 and 17), the control group shows markedly slower tumor development than the other two groups.Toward the end of the experiment, the differences are less distinct, as mice with large tumor volumes were euthanized.Throughout the study period, mice in the control group were euthanized less than the other two genotypes (Table 6).We also performed the analysis for two additional longitudinal mouse experiments in order to demonstrate the applicability and performance of the proposed method across similar experiments.First, the data for the animal tumor immunotherapy experiment are taken from Koziol et al. (1981).The data consist of three treatment groups of ten BALB/c mice (total of 30 mice) producing repeated measurements across 11 time points (from week 7 to week 21).
Towards the end of the study, 37% of the total mice dropped out.When mice died, tumor sizes ranged from 320.1mm 3 to 937.5mm 3 .One mouse survived with a tumor size of 1014.0mm 3 .Note that Koziol et al. (1981) did not consider the probability of the natural death or LOD in their analysis as opposed to our approach.There was one case of an LOD observation.Since the study examines simple tumor growth, we assume a linear trend is appropriate (refer to  Second, we examined our model using the data set from Tan et al. (2002).In this dataset, they analyzed a set of log-transformed data from an experiment in which mice were allocated to three anti-cancer treatments in a randomized experiment.Twenty mice were randomly assigned to one of three treatments with 7, 7, and 6 mice.Tan et al. (2002) used a pair-wise comparison approach to investigate whether the differences in the drug dosages altered the antitumor activities manifested by the tumor volumes.Note that their approach does not provide a closed form for the test statistic.In the available data set, whether mice died due to natural cause or euthanasia is not specified.Since Tan et al. (2002) indicated that mice were sacrificed when their tumor grew to four times its initial volume, we considered mice deaths occurring before the threshold as natural.In the data set they provided (Tan et al., 2002), we observed a total of 3 mice that died naturally and 1 mouse that was euthanized prior to the end of the experiment.Unlike the previous two examples, in this example, an antitumor regimen was employed to suppress tumor growth such that the tumor size actually reduced for a certain period and then re-grew (refer to Figure 1 in Tan et al., 2002).Thus, a model with cubic polynomials was used adapting (20).Using the proposed method, the ), strongly indicating a difference between treatments.

Concluding remarks
We have proposed an efficient method to carry out a group comparison for animal studies investigating longitudinal mammary tumor measurements.The proposed method employs a maximum likelihood methodology that accounts for the dependent data structure, the natural death rate, and the limit of detection.It is demonstrated that the proposed method allows us to carry out accurate inference (in terms of Type I error).The proposed method is sensitive in detecting differences between groups that are manifested by different death rates, even if no difference is shown in the observed tumor volumes.Our approach is easily implemented using R software packages for the various scenarios demonstrated in Section 4.

Appendix A. The conditional density based on the normal distribution
Let   be the distribution function and density function of the standard normal distribution, respectively.Based on the normal distribution, the conditional densities corresponding to the cases in Table 1, ( 11), ( 12), (13), and ( 14) are as follows.
(i) Case 1: The density of ij X is obtained using the inversion formula,  In (A-1), a line integration over t is carried out only on the real part of where Im(.) indicates the imaginary part of the complex number.The calculations of (A-1) and (A-2) are easily implemented and quite fast using available software such as R (http://www.r-project.org),which provides complex number arithmetic in basic functions.The numerator at the right-sided equation in ( is obtained by numeric integration, where

d
is an individual threshold tumor size indicating the euthanasia criterion.
(20) uses ideas presented in the context of the local maximum likelihood methodology (e.g.,Fan et al., 1998).For the purpose of illustration and in view of Figure2, we primarily consider a simple linear trend of the form j j

Table 2 :
Type I errors and 95th percentiles using the Monte-Carlo simulations for the test based on the linear model using the generalized least squares approach. and  are the intercept and slope for the linear function of j  . and  are defined in Error!Reference source not found.. N indicates the total sample size.The number of visits is 21 for all models to emulate the gene (Thoc1) deficiency study in Section 3. The target type I error and 95th percentile are 0.05 and 9

Table 3 :
The simulated Type I errors for various models.The intercept and slope for the linear function of j  were -6 and 0.45, respectively. and  are defined in Error!Reference source not found.. N and J indicate the total sample size and the number of visits.The death rate for each visit is 0.01.The survival rate is the rate of survived animals at time J.The cases indicated by * and ** utilized tdistributions with the degrees of freedom of 5 and 10, respectively, for the underlying distribution of ijk  .

2 
test statistic was 17.59 (degrees of freedom [ df ] = 4, -value p = 0.0015).This result strongly indicates the difference in tumor growth between groups.

Figure 1 :
Figure 1: The histograms of the log of tumor volumes at each week in the mammary tumor development data.The left panels are the control group and the right panels are the other two genotypes combined (Thoc1F/FPyMT+/-and Thoc1F/FMMTVCre+/-PyM T+/).

Figure 2 :
Figure 2: The top plot is the sample means of log-transformed tumor volumes (solid line: control group, dashed line: Thoc1F/FPyMT+/-, dotted line: Thoc1F/FMMTVCre+/-PyMT+/-).N/A and missing data were removed to make the plot.Before the transformation, the unit is cm3.The bottom plot is observations minus the sample means at each time point.

1 (
integration of the imaginary part is typically 0. The probability can also be obtained conveniently using the inversion theorem (Gil-Pelaez, 1951), is obtained by (A-1).

Table 1 : The possible cases for two consecutive observations. Unavailable observations due to the LOD is indicated by N/A whereas the check mark (✓) indicates a correctly observed value. The dot indicates a missing observation.
is independent of  , the conditional density of the measurement process for each case in Table1.For Case 1, as well to correctly incorporate the absence of natural death.Similarly, the likelihood under 1H is

Table 6 : Comparison of rates of euthanized animals between the control group and the other two groups for the data in Section 3. Difference at each time indicates the euthanasia rate in control group subtracted by the euthanasia rate in the other two genotypes. Thus, the negative values indicate that the control group died less.
(19)idering all these facts, the control group seems to have a slower tumor growth overall.Figure2describes mean tumor volumes (the top plot) and the observation minus the sample mean tumor volumes at each time point (the bottom plot) using the log-transformed data.The top plot demonstrates the gradual growth of the tumors.Note that the plot may not correctly depict the tumor growth trend, as it was obtained by removing missing cases.Notice also that the mice with a high tumor volume died toward the end of the study so that the rate of growth of the mean tumor volume tends to be slower; however, the real tumor growth should be steeper than what is depicted in the plot.An increasing trend of the variance toward the end of the experiment is well-demonstrated in the bottom plot demonstrating that Model (1) and the corresponding covariance structure are reasonable.It also shows that some sizable losses of experimental units started around week 20 due to sacrifices, indicating that the simple sample variance of the tumor size can underestimate the actual variance of the tumor size.Overall, the linear trend looks reasonable for modeling the growth of tumor.The likelihood ratio test(19)did not consider the natural death rate since all mice died by euthanasia.The