Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2753194

Formats

Article sections

- SUMMARY
- 1. Introduction
- 2. Model specification
- 3. Two-stage plug-in method
- 4. Joint Modeling of survival and TMA core-level data
- 5. Simulation
- 6. Case study in prostate cancer
- 7. Discussion
- REFERENCES

Authors

Related links

Stat Med. Author manuscript; available in PMC 2009 September 28.

Published in final edited form as:

Stat Med. 2008 May 20; 27(11): 1944–1959.

doi: 10.1002/sim.3217PMCID: PMC2753194

NIHMSID: NIHMS106719

See other articles in PMC that cite the published article.

Tissue Microarrays (TMAs) measure tumor-specific protein expression via high-density immunohistochemical staining assays. They provide a proteomic platform for validating cancer biomarkers emerging from large-scale DNA microarray studies. Repeated observations within each tumor result in substantial biological and experimental variability. This variability is usually ignored when associating the TMA expression data with patient survival outcome. It generates biased estimates of hazard ratio in proportional hazards models. We propose a Latent Expression Index (LEI) as a surrogate protein expression estimate in a two-stage analysis. Several estimators of LEI are compared: an Empirical Bayes (EB), a Full Bayes (FB), and a Varying Replicate Number (VRN) estimator. In addition, we jointly model survival and TMA expression data via a shared random effects model. Bayesian estimation is carried out using a Markov Chain Monte Carlo (MCMC) method. Simulation studies were conducted to compare the two-stage methods and the joint analysis in estimating the Cox regression coefficient. We show the two-stage methods reduce bias relative to the naive approach, but still lead to under-estimated hazard ratios. The joint model consistently outperforms the two-stage methods in terms of both bias and coverage property in various simulation scenarios. In case studies using prostate cancer TMA data sets, the two stage methods yields a good approximation in one data set while an insufficient one in the other. A general advice is to use the joint model inference whenever results differ between the two-stage methods and the joint analysis.

DNA microarray technology has enabled expression measurement of thousands of genes simultaneously, providing a platform for rapid screening of genomic biomarkers in cancer. Translating these discovery-type findings into clinical relevance is a more challenging task. Gene expression profiling studies using spotted cDNA arrays or Affymetrix GeneChip arrays often yield a few hundred candidate cancer genes that are associated with a phenotype. Only a small portion of these will be eventually validated for the corresponding expression changes at the protein level. The advent of Tissue Microarray (TMA) technology has provided a proteomic platform for such validation studies to find clinically useful biomarkers. TMA experiments measure tumor-specific protein expression via high-density immunohistochemical staining assays, allowing simultaneous evaluation of hundreds of patient samples in a single experiment [1]. Since their initial development, TMA-based expression studies have quickly become an integral part of cancer biomarker development [2-4].

A typical tissue array slide comprises up to 1000 tiny biopsy tissue elements, which we will refer to as cores, with multiple cores corresponding to repeated sampling from the same tumor. Expression measures on these replicate cores constitute the TMA core-level data. These can display substantial within-subject variability for both biological and experimental reasons. Biologically, for tumors that are highly infiltrative and heterogeneous in nature (e.g., prostate tumors), protein expression pattern can be quite variable. For example, cell proliferation genes often exhibit localized high expression within a tumor, indicating elevated aggressiveness and metastatic potential in the corresponding areas. Replicate sampling from various regions of the tumor is therefore important in capturing the underlying heterogeneity within a tumor. Experimental sources of the variability can come from a combination of probe affinity, measurement imprecision, and further missing data due to insufficient sampling. Without appropriately accounting for these variabilities, the noise-prone measurements tend to attenuate the prognostic value of a potential biomarker in predicting disease outcome. The lack of a model-based approach for TMA core-level expression data to effectively model the intra-tumor variation has motivated us to carry out a full investigation.

A good analogy for understanding TMA data structure is from probe-level data generated by the Affymetrix GeneChip arrays. GeneChip arrays measure gene expression at the mRNA transcripts level. The probe-level data refer to the replicate expression measures on a set of 16-20 small oligonucleotide probe sets derived for a target gene. The biological variation comes primarily from these oligonucleotide probe sequence variants, while experimental variations arise during the process of slide printing, hybridization and optical reading. Li *et al.* [6] reported that the variation of a specific probe across multiple arrays could be considerably smaller than the variance across probes within a probe set. Modeling Affymetrix probe-level data has generated much attention [6, 7] as the technology has become more mature and widely used.

Similarly in TMA experiments, modeling within-tumor protein expression heterogeneity is an important problem. In these tissue-based experiments, the variation across core samples within a tumor can be substantially larger than the variation observed across subjects. Etzioni *et al.* [5] used a compositional analysis to model such heterogeneity, and compared the proportion of cells stained at different intensity levels between normal and tumorous tissues. In this study, we focus on the effect of modeling intra-tumor variation in the context of predicting patient survival outcome. In a latent variable modeling framework, we assume that an underlying `true' expression value predicts survival. In real experiments, this `true' expression can not be precisely measured due to sampling variabilities and measurement imprecision. Instead, we observe the core-level expression measurements that are subject to these measurement errors. In a two-stage method, we adapt ideas from measurement error modeling and propose a latent expression index (LEI) to approximate the underlying true value, and focus on its behavior in proportional hazards models. We adapt an empirical Bayes estimator [8] to 1) incorporate important clinical parameters such as Gleason score and pathological stage of the tumor and 2) adjust for the varying number of cores. We further establish a joint model for TMA core-level data and survival outcome via a shared random effect. There is a large literature on joint modeling of longitudinal data and survival [9-16]. These methods have been developed predominantly for modeling survival and CD4 counts in AIDS patients; here their application to Tissue Microarray data in cancer biomarker studies is novel. Using both simulations and two published TMA data sets, we compare the performances of the naive, two-stage LEI, and the joint model approach in obtaining parameter estimates and associated inference.

The article is organized as follows. Section 2 specifies notation and models for the TMA core-level expression data and for the patient survival data. Section 3 introduces LEI and its use in a two-stage method. Section 4 presents the joint modeling approach and the Bayesian estimation framework. Simulation results to compare the performances of these methods are then discussed in Section 5. Case studies using two prostate cancer TMA data sets are presented in Section 6. Further discussion can be found in Section 7.

Let ${X}_{i}^{\ast}$ be the latent expression value for a biomarker in tumor *i*, *i* = 1, … , *n*. Assuming the observed TMA core-level measurement for the *j*th core in the *i*th tumor is

$${X}_{ij}={X}_{i}^{\ast}+{U}_{ij},\phantom{\rule{thickmathspace}{0ex}}j=1{,}_{\cdots},{r}_{i};i=1{,}_{\cdots},n,$$

(2.1)

where we assume ${X}_{i}^{\ast}~N({\mu}_{x\ast},{\sigma}_{x\ast}^{2})$. The mean μ_{x*} is a linear function of clinical covariates: ${\mu}_{x\ast}={\theta}_{0}+\theta {\mathit{Z}}_{i}^{\prime}$, where **Z*** _{i}* = (

Let the observed survival time for patient *i* (*i* = 1, …, *n*) be *T _{i}* = min(

$$\lambda \left(t\right)={\lambda}_{0}\left(t\right){e}^{{\beta}^{\ast}{X}_{i}^{\ast}+{\mathrm{bZ}}_{i}^{\prime}},$$

(2.2)

where λ_{0}(*t*) is the baseline hazard function, β* is the regression coefficient associated with the `true' protein expression, and **b** = (*b*_{1}, … , *b _{p}*) is the row vector of coefficients associated with the vector of clinical variables ${\mathit{Z}}_{i}^{\prime}$. We also consider a parametric Weibull regression model with the following form for the hazard function:

$$\lambda \left(t\right)=\gamma {t}^{\gamma -1}{e}^{{\beta}_{0}+{\beta}^{\ast}{X}_{i}^{\ast}+{\mathrm{bZ}}_{i}^{\prime}},$$

(2.3)

Given the rare disease assumption and the assumption that the measurement error *U _{ij}* has no predictive value, i.e., $\lambda (t\mid {\mathbf{X}}_{i},{X}_{i}^{\ast})=\lambda (t\mid {X}_{i}^{\ast})$ Prentice (1982) introduced the induced hazard rate function for Cox regression model with covariates measured with error. When normality is assumed for the error distribution, a first-order bias-correction method is to replace ${X}_{i}^{\ast}$ with its conditional mean

$$\lambda (t\mid {\mathrm{X}}_{i})={\lambda}_{0}\left(t\right){e}^{{\beta}^{\ast}E[{X}_{i}^{\ast}\mid {\mathrm{X}}_{i},{\mathrm{Z}}_{i}]+{\mathrm{bZ}}_{i}^{\prime}}$$

(3.1)

and estimate β* by maximizing the corresponding partial likelihood. Here ${\mathbf{X}}_{i}=({X}_{i1},{X}_{i2},\cdots ,{X}_{i{r}_{i}})$ denotes the vector of observed expression measures for tumor *i*. Note that the regression calibration model in (3.1) is an approximation to the true model in (2.2) due to the rare disease assumption and the omission of a higher order term implicit in (2.2). Details of such first-order bias-correction methods can be found in [17, 18]. Define the Latent Expression Index (LEI) to be an estimate of the conditional mean, ${\mathrm{LEI}}_{i}=\widehat{E}[{X}_{i}^{\ast}\mid {\mathbf{X}}_{i},{\mathbf{Z}}_{i}]$, for each subject *i*. A two-stage plug-in method can be described by the following algorithm: 1) Compute LEI_{i} (*i* = 1, …, *n*) as a surrogate expression estimate that adjusts for measurement error; and 2) Apply the Cox or Weibull regression model using LEI_{i} to obtain an estimate of β* and the associated standard error.

In the next section, we describe methods for computing LEI for tissue microarray data. These include an empirical Bayes estimator conditional on clinical covariates, a full Bayes approach and a Varying Replicate Number (VRN) method as an extension to adjust for the number of cores per tumor.

Express (2.1) as a mixed effects model

$${X}_{ij}={\theta}_{0}+\theta {\mathrm{Z}}_{i}^{\prime}+{v}_{i}+{U}_{ij},$$

(3.2)

where ${\nu}_{i}~N(0,{\sigma}_{x\ast}^{2})$. To make a connection between (3.2) and (2.1), ${X}_{i}^{\ast}$ is the unobserved “true” protein expression level, which can be expressed as ${X}_{i}^{\ast}={\theta}_{0}+\theta {\mathbf{Z}}_{i}^{\prime}+{\nu}_{i}$. Subsequently, ν_{i}, which is normally distributed with mean zero and variance, has the interpretation of a mean centered “true” protein expression level. The empirical Bayes estimator can then be derived as

$$LE{I}_{i}^{eb}={\widehat{\gamma}}_{i}{\stackrel{\u2012}{X}}_{i}+(1-{\widehat{\gamma}}_{i})({\widehat{\theta}}_{0}+\widehat{\theta}{\mathrm{Z}}_{i}^{\prime})$$

(3.3)

where ${\widehat{\gamma}}_{i}\equiv {\widehat{\sigma}}_{x\ast}^{2}{({\widehat{\sigma}}_{x\ast}^{2}+{\widehat{\sigma}}_{u}^{2}{r}_{i}^{-1})}^{-1}$ is the attenuation factor [18]. Parameter estimates $\{{\widehat{\theta}}_{0},\widehat{\theta},{\widehat{\sigma}}_{u}^{2},{\widehat{\sigma}}_{x\ast}^{2}\}$ can be obtained by fitting a mixed effects model as described in (3.2), using a restricted maximum likelihood (REML) approach [19, 20].

The empirical Bayes estimator conditions on the set of parameter estimates derived from the data. The uncertainty of these estimates are not accounted for in LEI^{eb}. For this reason, we also consider a full Bayesian estimator, ${\mathit{LEI}}_{i}^{fb}$, where standard conjugate hyperprior distributions are adopted: ${\sigma}_{u}^{-2},{\sigma}_{x\ast}^{-2}~\Gamma ({r}_{0},{\gamma}_{0})$. The full Bayes estimator

$$LE{I}_{i}^{fb}={\stackrel{~}{\theta}}_{0}+\stackrel{~}{\theta}{\mathrm{Z}}_{i}^{\prime}+{\stackrel{~}{v}}_{i},$$

is then based on the posterior inference from model (3.2) where $\{{\stackrel{~}{\theta}}_{0},\stackrel{~}{\theta},{\stackrel{~}{\nu}}_{i}\}$ are the posterior means given the data.

In a typical TMA construction, *K _{i}* cores are placed on the array for each tumor

$$\text{logit}P({M}_{ij}=1)={\psi}_{0i}+\psi {\mathrm{Z}}_{i}^{\prime},$$

(3.4)

where ${\psi}_{0i}~N({\psi}_{0},{\sigma}_{\psi}^{2})$, **Z*** _{i}* = (

$${r}_{i}\sim \text{Binomial}\left({K}_{i},{p}_{i}=\frac{{e}^{{\psi}_{oi}+\psi {\mathrm{Z}}_{i}^{\prime}}}{1+{e}^{{\psi}_{oi}+\psi {\mathrm{Z}}_{i}^{\prime}}}\right).$$

(3.5)

The expression index under the VRN model is then derived by averaging over all the possible values of (*r _{i}*,

$$\begin{array}{cc}\hfill LE{I}_{i}^{vrn}& ={E}_{({r}_{i},{K}_{i})}E[{X}_{i}^{\ast}\mid {\mathbf{X}}_{i},{\mathbf{Z}}_{i},{r}_{i},{K}_{i}]\hfill \\ \hfill & =\sum _{s=1}^{R}\sum _{m=0}^{s}\left\{\frac{{\widehat{\sigma}}_{{x}^{\ast}}^{2}m}{{\widehat{\sigma}}_{{x}^{\ast}}^{2}m+{\widehat{\sigma}}_{u}^{2}}{\stackrel{\u2012}{X}}_{i}+\frac{{\widehat{\sigma}}_{u}^{2}}{{\widehat{\sigma}}_{{x}^{\ast}}^{2}m+{\widehat{\sigma}}_{u}^{2}}({\widehat{\theta}}_{0}+\widehat{\theta}{Z}_{i}^{\prime})\right\}\times \left({\scriptstyle \begin{array}{c}\hfill s\hfill \\ \hfill m\hfill \end{array}}\right){\left(\frac{{e}^{{\widehat{\psi}}_{0i}+\widehat{\psi}{\mathbf{Z}}_{i}^{\prime}}}{1+{e}^{{\widehat{\psi}}_{0i}+\widehat{\psi}{\mathbf{Z}}_{i}^{\prime}}}\right)}^{m}{\left(\frac{1}{1+{e}^{{\widehat{\psi}}_{0i}+\widehat{\psi}{\mathbf{Z}}_{i}^{\prime}}}\right)}^{s-m}\widehat{P}({K}_{i}=s)\hfill \end{array}$$

(3.6)

Here *R* is taken to be max_{i }*K _{i}*. An additional assumption for the above is that the expression measures do not correlate with

In a relatively balanced TMA array where the number of replicate measures *r _{i}* does not vary much across subjects, ${\widehat{\gamma}}_{i}\equiv {\widehat{\sigma}}_{x\ast}^{2}{({\widehat{\sigma}}_{x\ast}^{2}+{\widehat{\sigma}}_{u}^{2}{r}_{i}^{-1})}^{-1}$ is an approximately constant adjustment factor. The amount of shrinkage in LEI

The two-stage approaches described above are attractive for their simplicity and straightforward interpretation. They require minimal computation and can be easily implemented using existing statistical packages. However, there are major limitations for the two-stage method [23]. First, the two-stage method involves a first order approximation and ignores the second-order term ${\beta}^{{\ast}^{2}}{\sigma}^{2}({X}_{l}^{\ast}\mid \mathbf{X})$ in the induced hazard rate function (3.1). As we will illustrate in the simulation study, such approximation works well when β* is close to zero, but otherwise lead to sizeable bias in ${\widehat{\beta}}^{\ast}$. Second, parameter estimates in the second stage do not account for the uncertainty in estimating LEI in the first stage. The associated standard error for ${\widehat{\beta}}^{\ast}$ will be over-optimistic. Given these considerations, it is desirable to make inference based on the joint likelihood of the failure time and TMA expression data. In this article, we adopt a shared random effect model to induce correlation between the TMA data and the survival outcome.

Given the measurement model specified in (3.2) for the TMA data, we write the proportional hazards model for the survival outcome as

$$\lambda \left(t\right)={\lambda}_{0}\left(t\right)\phantom{\rule{thickmathspace}{0ex}}\mathrm{exp}({\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{v}_{i}).$$

(4.1)

The parameter ν* _{i}* constitutes a shared random effect that connects the measurement model (3.2) and the survival outcome model (4.1). The expression data and survival times are then assumed to be independent given ν

$${L}_{joint}=\prod _{i=1}^{n}f({t}_{i},{\delta}_{i}\mid {v}_{i},{\mathbf{z}}_{i})\times \prod _{i=1}^{n}f({X}_{i}\mid {v}_{i},{\mathbf{Z}}_{i})f\left({v}_{i}\right)={L}_{SU\phantom{\rule{thickmathspace}{0ex}}RV}\times {L}_{M\phantom{\rule{thickmathspace}{0ex}}E},$$

(4.2)

where

$$\begin{array}{cc}\hfill & {L}_{M\phantom{\rule{thickmathspace}{0ex}}E}=\prod _{i=1}^{n}\left(\prod _{j=1}^{{r}_{i}}\frac{1}{\sqrt{2\pi {\sigma}_{u}^{2}}}\mathrm{exp}\left\{-\frac{{({x}_{ij}-{\theta}_{0}-\theta {\mathbf{Z}}_{i}^{\prime}-{v}_{i})}^{2}}{2{\sigma}_{u}^{2}}\right\}\right)\frac{1}{\sqrt{2\pi {\sigma}_{{x}^{\ast}}^{2}}}\mathrm{exp}\left\{-\frac{{v}_{i}^{2}}{2{\sigma}_{{x}^{\ast}}^{2}}\right\};\hfill \\ \hfill & {L}_{SU\phantom{\rule{thickmathspace}{0ex}}RV}=\prod _{i=1}^{n}\prod _{l=1}^{L}{\lambda}_{l}^{dl}\mathrm{exp}(\sum _{i\in {D}_{l}}{\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{v}_{i})\mathrm{exp}(-{\lambda}_{t}\sum _{i\in {R}_{l}}{\Delta}_{il}{e}^{{\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{v}_{i}})\hfill \end{array}$$

(4.3)

We used a piecewise constant hazards model in which the time axis is partitioned into *L* disjoint intervals, *I*_{1}, …, *I _{L}*, where

Alternatively, a parametric Weibull model can be assumed for the survival outcome using the following hazard function:

$$\lambda \left(t\right)=\gamma {t}^{\gamma -1}\mathrm{exp}({b}_{0}+{\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{v}_{i}).$$

(4.4)

When γ = 1, the above reduces to exponential distribution with constant failure rate $\mathrm{exp}({b}_{0}+{\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{\nu}_{i})$. The survival time component of the joint likelihood in (4.3) is then replaced by

$${L}_{SU\phantom{\rule{thickmathspace}{0ex}}RV}=\prod _{i=1}^{n}{\left\{\gamma {t}_{i}^{\gamma -1}{e}^{\gamma ({b}_{0}+{\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{v}_{i})}\right\}}^{{\delta}_{i}}\mathrm{exp}(-{e}^{\gamma ({b}_{0}+{\mathbf{bZ}}_{i}^{\prime}+{\beta}^{\ast}{v}_{i})}{t}_{i}^{\gamma}).$$

(4.5)

In a Bayesian estimation framework, the following prior distributions are specified for the model parameters:

$$\begin{array}{cc}\hfill & ({\beta}^{\ast},{\theta}_{0},\theta ,{b}_{0},b)\phantom{\rule{thickmathspace}{0ex}}\sim \phantom{\rule{thickmathspace}{0ex}}N\phantom{\rule{thickmathspace}{0ex}}({\mu}_{0},{\sigma}_{0}^{2});\hfill \\ \hfill & ({\sigma}_{u}^{-2},{\sigma}_{{x}^{\ast}}^{2})\phantom{\rule{thickmathspace}{0ex}}\sim \phantom{\rule{thickmathspace}{0ex}}\Gamma ({r}_{0},{\gamma}_{0});\hfill \\ \hfill & (\gamma ,{\lambda}_{l},\phantom{\rule{thickmathspace}{0ex}}l=1{,}_{\cdots},L)\phantom{\rule{thickmathspace}{0ex}}\sim \phantom{\rule{thickmathspace}{0ex}}\Gamma ({r}_{0},{\gamma}_{0}).\hfill \end{array}$$

(4.6)

The prior choices in (4.6) render proper noninformative prior distributions with the hyperparameters chosen to achieve large variances (wide spread). In particular, ${\mu}_{0}=0,{\sigma}_{0}^{2}=10000,{r}_{0}=0.001,{\gamma}_{0}=0.001$. The values of *r*_{0} and γ_{0} are set to be small to impose noninformativeness (although there have been criticisms that the resulting posterior distributions remain highly sensitive to the choice of *r*_{0} and γ_{0}). We had tested different sets of values for *r*_{0} and γ_{0}: (0.1, 0.1), (0.01, 0.01), and (0.001, 0.001). When the sample size is moderate to large, the variance estimates, ${\widehat{\sigma}}_{u}$ and ${\widehat{\sigma}}_{{x}^{\ast}}$, are well-behaved and not significantly affected by the choice of *r*_{0} and γ_{0}. Samples from the posterior distribution are obtained using Markov Chain Monte Carlo (MCMC) methods. In the likelihood function of (4.2), the Bayesian estimation procedure treats ν* _{i}* as missing data and imputes it at each MCMC iteration via posterior sampling from the conditional distribution of ν

The additive measurement error model in (3.2) with one covariate *Z*_{1i} is used to simulate the expression measure *X _{ij}*,

Computation of *LEI ^{eb}*,

The simulation results are summarized in Table I. For β* = 1 in the survival models, the naive approach (using ${\stackrel{\u2012}{X}}_{i}$ as a surrogate expression) attenuates the true effect size by around 25%. The coverage probability of a nominal 95% confidence interval of ${\widehat{\beta}}^{\ast}$ is 0.10 at best. The two stage methods (LEI) achieved a considerable bias correction by adjusting for the measurement error in the LEI imputation. The joint modeling approach gives the best estimate ${\widehat{\beta}}^{\ast}=1.03$ and a coverage probability of 95% compared to the truth.

Next a larger effect size is simulated (β* = 2). The bias in the two-stage approaches due to the first-order approximation is evident. Overall the two-stage methods generate less biased ${\widehat{\beta}}^{\ast}\u2019\mathrm{s}$ compared to the naive estimate. Nevertheless, these are 15-25% smaller than the true β*. The coverages are poor for the two-stage approaches. The joint modeling approach should be advocated in this scenario for inference. Notice that the joint model estimates of β* are slightly bigger than 2, especially under the Weibull distribution. This may be driven by the prior distributions adopted for the parameters under the Bayesian estimation scheme. We have observed that such a difference disappears when the sample size gets larger.

In this study, we consider two prostate tumor tissue microarray data sets. The α-Methylacyl CoA racemase (AMACR) is a peroxisomal and mitochondrial enzyme that plays an important role in fatty acid metabolism. AMACR has been shown to consistently overexpress in prostate tumors [26]. Rubin *et al.* [4] profiled AMACR protein expression using a TMA constructed on 203 prostate tumors from a surgical cohort who underwent radical prostatectomy at the University of Michigan as a primary therapy for clinically localized prostate cancer diagnosed between 1994 and 1998. They found AMACR is a significant predictor of the Prostate Specific Antigen (PSA) failure in these 203 patients.

The second biomarker evaluated in this study is BM28. BM28 encodes a highly conserved mini-chromosome maintenance protein (MCM) that is involved in the initiation of genome replication. Bismar *et al.* [27] profiled a total of 41 genes (including BM28) in a TMA-based proteomic study. They identified a 12-gene model showing the expression combination of the twelve genes significantly associates with tumor progression and PSA failure in a set of 79 men following surgery for clinically localized prostate cancer. The expression level of BM28 however did not show significant prognostic value in their analysis. We chose BM28 to evaluate the possibility of its being a false negative biomarker due to measurement error.

For the AMACR data, an average of *K _{i}* =5.5 (range: 2 to 12) tissue core specimens were taken from each tumor sample and put on a tissue array for immunohistochemical staining. After initial diagnostic evaluation of each core, an average of 29% (range: 0-86%) of the cores were excluded due to reasons discussed earlier, leading to 5.5% missing subjects. The BM28 data has a similar tissue array design. A quantitative imaging analysis of the staining intensity was obtained using the ACIS II (Chromavision, San Juan Capistrano, CA) system. The intensity level ranges from 0 to 255 chromogen intensity units, and is transformed using the natural logarithm (one unit added to avoid taking logarithm of 0) and normalized to have mean zero and standard deviation one. Disease recurrence is defined as a serum PSA increase >0.2ng/mL after radical prostatectomy. Censored observations are those free of the recurrence at the time of last follow-up.

Figure 1 illustrates a considerable amount of measurement error in the core-level expression data from the two TMA experiments. In the AMACR data, methods-of-moments estimates of the variance components are ${\widehat{\sigma}}_{u}^{2}=0.46$ and ${\widehat{\sigma}}_{{x}^{\ast}}^{2}=0.54$. In the BM28 data, the methods-of-moments estimates of the variance components are ${\widehat{\sigma}}_{u}^{2}=0.62$ and ${\widehat{\sigma}}_{{x}^{\ast}}^{2}=0.21$. The within-subject variation is almost three times the between-subject variation.

Variance plots to represent the within-subject variation in the TMA core-level expression data. A) The AMACR data. Estimates of the variance components are: ${\widehat{\sigma}}_{u}^{2}=0.46$ and ${\widehat{\sigma}}_{x\ast}^{2}=0.54$. B)

In a Cox proportional hazards model context, we did a simple simulation where ${X}_{i}^{\ast}~N(0,1)$, β* = 1. A naive estimator ${\stackrel{\u2012}{X}}_{i}={r}_{i}^{-1}{\sum}_{j=1}^{{r}_{i}}{X}_{ij}$ —the average core-level expression for tumor *i*— is used as the surrogate expression to replace ${X}_{i}^{\ast}$ in (2.2). We simulate situations where measurement error is small $({\sigma}_{u}^{2}=0.1)$, moderate $({\sigma}_{u}^{2}=0.5)$, and large $({\sigma}_{u}^{2}=1)$. Figure 2 shows various degrees of regression attenuation in the estimate of β* as a function of replicate number *r _{i}* and the amount of error ${\sigma}_{u}^{2}$. When the parameter values are set to resemble the AMACR data, the naive estimate of β* is approximately 30% smaller than the true value. With the current TMA construction protocol specifying three cores per subject due to economic and tissue-preservation reasons, and a great amount of within-subject variability routinely observed in the core-level expression data, Figure 2 effectively conveys the importance of modeling measurement error in TMA data. In the following two sections, we implement the measurement error models to demonstrate how statistical inference differs from the previous results.

In prostate cancer, Gleason score, pathologic stage and tumor size are among the most important clinical parameters. We include these as clinical covariates **Z**_{i} to adjust in the measurement model (3.2), the replicate number model (3.4), and the survival outcome model (4.1). In the measurement model,${\widehat{\theta}}_{TumorSize}=.32$ with an associated standard error of 0.13, indicating a marginal association of tumor size with AMACR expression level. In the replicate number model,${\widehat{\psi}}_{TumorSize}=.72$ with an associated standard error of 0.16, which is consistent with our expectation that a larger tumor sample provides more abundant number of cores.

Table II lists the estimates and associated standard errors (posterior standard deviation) of β* in the outcome model. The measurement error adjustment has significantly improved upon the naive estimate. The error-adjusted $\widehat{\beta}\ast $ is around 0.75 $\left(\widehat{se}(\widehat{\beta}\ast )=.26\right)$, approximately 31% larger in absolute value than the naive estimate which is 0.57 $\left(\widehat{se}(\widehat{\beta}\ast )=.20\right)$. The amount of attenuation in β* is quite consistent with what we conclude from the simulated datasets in the previous section. In this dataset, the two-stage methods (LEI^{eb}, LEI^{fb}, LEI^{vrn}) perform equally well as the joint modeling approach. The simplicity and computational efficiency of LEI serves as a satisfactory core-level expression index for AMACR. However as mentioned earlier, the two-stage methods are based on a first-order approximation, the accuracy of which is largely driven by the size of β* and the ratio of within- and between-subject variation. As will be shown in the other data example, the two-stage methods will not be always a suitable approach.

Kaplan-Meier curves are useful as a graphical representation of the prognostic value of a biomarker. We examined these plots by dividing the subjects into different risk groups based on the values of AMACR expression estimates derived under each method. In Figure 3(a), subjects with AMACR high, median, and low expression groups based on *LEI ^{vrn}* and the joint model estimates (C and D respectively) are significantly better separated in terms of probability of recurrence-free survival, when compared to that using the naive mean estimates (A).

The measurement model indicates a marginal association of pathologic stage of the tumor with BM28 expression: ${\widehat{\theta}}_{PathStage}=0.59\left(\widehat{se}\left(\theta \right)=0.22\right)$. The replicate number model again suggests a strong dependence on the size of the tumor: ${\widehat{\psi}}_{Tumoresize}=1.12\left(\widehat{se}\left(\psi \right)=0.43\right)$.

In Table II, the differences under various models are more discernable in this datasets. First, the Weibull and piecewise exponential model overall generate slightly different results given a small sample size (n=52). Second, the empirical Bayes estimate differs substantially from the full Bayes LEI estimate. It is likely due to the large uncertainties in the parameter estimates that determine LEI^{eb}. Finally, the bias introduced by the first-order approximation is prominent here. Both the coefficient β* and the noise ratio in this dataset are much larger in magnitude compared to the AMACR data example. In this case, the two-stage methods alleviate regression attenuation, only to a limited extent. The joint model should be used for parameter estimation and associated inference.

Figure 3(b) plots the Kaplan-Meier curves using different expression estimates. The 5-year PSA recurrence-free survival probability is 0.95 $(\widehat{se}=0.05)$ versus 0.58 $(\widehat{se}=0.10)$ for low and high BM28 expression estimated by the joint model. Adjusting for measurement error in this dataset has made a dramatic change in the conclusion about the prognostic value of BM28, compared to the naive method.

Figure 4 compares the naive and the joint model expression estimates, plotted against the survival time on the x-axis. The mean expression ± two standard deviations is plotted for each individual. Two improvements under the joint model are clear: 1) the noise, represented by the error bars, is greatly reduced via the joint modeling, and 2) the mean expression levels are distributed more tightly around a regression line, accentuating the relationship of the TMA expression data and survival time.

In TMA data analysis, statistical methods often focus on downstream models in predicting disease outcome assuming ${\stackrel{\u2012}{X}}_{i}$ is a suffcient expression summary measure [29, 30]. Relatively little attention has been given to the modeling of within-tumor variation in these TMA experiments. As we have shown in this paper with real data examples, analysis ignoring intra-tumor variation can lead to false negative results, causing a tremendous waste of valuable tissue resource and experimental costs. In this study, we propose and compare a two-stage approach and a joint model to analyze tissue microarray data in a measurement error framework. The two-stage approach computes a Latent Expression Index (LEI) in the first stage and then associate the LEI with patient survival information using a proportional hazards model in the second stage. We further adjust for clinical/pathological covariates *Z _{i}* and the number of repeated measures (

Although the two-stage approach is attractive for its simplicity, interpretability and straightforward computation, practitioners need to be aware of the associated statistical limitations. As we discussed earlier, the regression calibration approach is an approximation of the true model which reduces bias relative to the naive method but still leads to underestimated β*. In addition, the survival information is not used in the first stage to compute LEI, which can cause bias and effciency loss in estimating β* in the second stage. Finally, the uncertainty of estimating the LEI quantity is not assimilated in the second stage, leading to over-optimistic standard error estimates. For a review of this topic, see Tsiatis and Davidian [23]. In simulation studies, we show the two-stage methods reduce bias relative to the naive approach in estimating the Cox regression coefficient, but still lead to under-estimated β*. The joint analysis, however, outperforms the two-stage methods in terms of the bias and coverage property in various simulated scenarios. In case studies, we applied the error methods to prostate TMA data sets where the two-stage approach and the joint analysis generate 1) very similar results in the AMACR data set and 2) disparate results in the BM28 data set. While the first case study demonstrates a scenario where the two-stage approach is a good approximation to the true model, the second suggests one where such approximation is not so accurate (the higher order term can not be ignored when β* is a large value). Inference should be based on the joint analysis whenever results differ.

A final notion is regarding the impact of the technology advancement on TMA data analysis. Recent advances in quantitative assessment of the immunohistochemical staining provide precise, objective, and reproducible protein expression measurements. Compared to the conventional pathologist scoring on an ordinal scale, the Chromavision system used in our data examples enables quantification of the antigen level on a continuous scale, free of the subjectivity associated with pathologist-based visual scoring system. We should point out that the Chromavision system not only quantifies the amount of protein expression, but also provides a more accurate measure of the percentage of tissue that express the protein. A statistical problem of interest is to extend our model to incorporate the percentage information when summarizing the expression profile of a tumor. Another quantitative system, AQUA [28], which stands for Automated Quantitative Analysis, measures fluorescence signals, leading to higher sensitivity to very low antibody concentrations. In addition, it allows the separation of tumor from stromal elements and the sub-cellular localization of signals for a co-localization of the antigens in different cell compartments. As the technology is becoming more and more refined, statistical models underpinning both biological and experimental issues are in great need.

We thank Dr. Arul Chinnaiyan at the Pathology and Urology Department at the University of Michigan Medical School and Dr. Mark Rubin at the Pathology Department of Weill Medical College of Cornell University for making the prostate tissue microarray datasets available.

1. Kononen J, Bubendorf L, Kallioniemi A, Barlund M, Schraml P, et al. Tissue microarrays for high-throughput molecular profiling of tumor specimens. Nature Medicine. 1998;4:844–7. [PubMed]

2. Divito KA, Berger AJ, Camp RL, Dolled-Filhart M, Rimm DL, Kluger HM. Automated quantitative analysis of tissue microarrays reveals an association between high bcl-2 expression and improved outcome in melanoma. Cancer Research. 2004;64:8773–7. [PubMed]

3. Seligson DB, Horvath S, Shi T, Yu H, Tze S, Grunstein M, Kurdistani SK. Global histone modification patterns predict risk of prostate cancer recurrence. Nature. 2005;435:1262–6. [PubMed]

4. Rubin MA, Bismar TA, Andrén O, Mucci L, Kim R, Shen R, Ghosh D, Wei J, Chinnaiyan AM, Adami H, Kantoff PM, J. J. Decreased *α*-Methylacyl CoA racemase expression in localized prostate cancer is associated with an increased rate of biochemical recurrence and cancer-specific death. Cancer Epidemiology and Biomarkers Prevention. 2005;14:1424–31. [PubMed]

5. Etzioni R, Hawley S, Billheimer D, True LD, Knudsen B. Analyzing patterns of staining in immunohistochemical studies: application to a study of prostate cancer recurrence. Cancer Epidemiology and Biomarkers Prevention. 2005;14:1040–6. [PubMed]

6. Li C, Wang WH. Model-based analysis of oligonucleotide arrays: Expression index computation and outlier detection. Proceedings of the National Academy of Sciences. 2001;98:31–6. [PubMed]

7. Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP. Summaries of Affymetrix genechip probe level data. Nucleic Acids Research. 2003;31:1–8. [PMC free article] [PubMed]

8. Tsiatis AA, DeGruttola V, Wulfsohn MS. Modeling the relationship of survival to longitudinal data measured with error. applications to survival and CD4 counts in patients with AIDS. Journal of the American Statistical Association. 1995;90:27–37.

9. Wulfsohn MS, Tsiatis AA. A joint modeling for survival and longitudinal data measured with error. Biometrics. 1997;53:330–39. [PubMed]

10. Faucett CJ, Thomas DC. Simultaneously modeling censored survival data and repeatedly measured covariates: a Gibbs sampling approach. Statistics in Medicine. 1996;15:1663–85. [PubMed]

11. Henderson R, Diggle P, Dobson A. Joint modeling of longitudinal measurements and event time data. Biostatistics. 2000;1:465–80. [PubMed]

12. Wang Y, Taylor JMG. Joint modeling longitudinal and event time data with application to acquired immunodeficiency syndrome. Journal of the American Statistical Association. 2001;96:895–905.

13. Xu J, Zeger SL. Joint analysis of longitudinal data comprising repeated measures and times to events. Applied Statistics. 2001;50:375–87.

14. Brown ER, Ibrahim JG. A Bayesian semiparametric joint hierarchical model for longitudinal and survival data. Biometrics. 2003;59:221–28. [PubMed]

15. Guo X, Carlin B. Separate and joint modeling of longitudinal and event time data using standard computer packages. The American Statistician. 2004;58:1–9.

16. Tadesse MG, Ibrahim JG, Gentleman R, Chiaretti S, Ritz J, Foa R. Bayesian error-in-variable survival model for the analysis of genechip arrays. Biometrics. 2005;61:488–97. [PubMed]

17. Prentice RL. Covariate measurement errors and parameter estimation in a failure time regression model. Biometrika. 1982;69:331–42.

18. Carroll RJ, Ruppert D, Stefanski LA. Measurement error in nonlinear models. Chapman and Hall; 1995.

19. Harville DA. Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association. 1977;72:320–38.

20. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38:963–74. [PubMed]

21. Lindstrom MJ, Bates DM. Nonlinear mixed effects models for repeated measures data. Biometrics. 1990;46:673–87. [PubMed]

22. McCulloch CE. Maximum likelihood variance components estimation for binary data. Journal of the American Statistical Association. 1994;89:330–35.

23. Tsiatis AA, Davidian M. Joint modeling of longitudinal and time-to-event data: an overview. Statistica Sinica. 2004;14:809–34.

24. Spiegelhalter D, Thomas A, Best N, Lunn D. WinBUGS 1.4 Manual. 2003 Http://www.mrcbsu.cam.ac.uk/bugs/winbugs/contents.shtml.

25. Thomas A. BRugs: An R interface to OpenBUGS. Version 1.0 User's Manual. 2004 Http://mathstat.helsinki.fi/openbugs/

26. Rhodes DR, Barrette TR, Rubin MA, Ghosh D, Chinnaiyan AM. Meta-analysis of microarrays: interstudy validation of gene expression profiles reveals pathway dysregulation in prostate cancer. Cancer Research. 2002;62:4427C33. [PubMed]

27. Bismar T, Demichelis F, Riva A, Kim R, Varambally S, He L. Defining aggressive prostate cancer using a 12-gene model. Neoplasia. 2006;8:59–68. [PMC free article] [PubMed]

28. Camp RL, Chung GG, Rimm DL. Automated subcellular localization and quantification of protein expression in tissue microarrays. Nature Medicine. 2002;8:1323–28. [PubMed]

29. Horvath S, Liu X, Huang Y, Kim HL, Seligson DB. Forest-based methods for analyzing tissue microarray. ASA Proceedings of the Joint Statistical Meetings, 1863-1869 American Statistical Association (Alexandria, VA) 2003

30. Liu X, Minin V, Huang Y, Seligson DB, Horvath S. Statistical methods for analyzing tissue microarray data. Journal of Biopharmaceutical Statistics. 2004;14:671–85. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |