flowchart LR Absorption(Absorption) --> A(Circulatory System) A -- Distribution --> B(Gastrointestinal System) B -- Metabolism --> A B --> Excretion(Excretion) style Absorption fill:transparent,stroke:transparent style Excretion fill:transparent,stroke:transparent style A fill:#E7D6AE,stroke:#C9B27D,color:#000 style B fill:#E7D6AE,stroke:#C9B27D,color:#000
7 Early Phase Studies
Early Phase Studies
Dose-Finding Studies
Bayesian Methods
Pharmacokinetics
Correlation and Agreement
Overview
Early phase studies are designed to characterize the dose-response relationship of a new treatment, assess safety and pharmacokinetics, and identify an appropriate dose for later phase trials. This lesson covers the design and analysis considerations for dose-finding (Phase I) studies, Bayesian and likelihood methods for dose escalation, and methods for assessing correlation and agreement.
Objectives
Upon completion of this lesson, you should be able to:
- State the components of a typical dose-finding design.
- State the objectives of a pharmacokinetic model.
- Describe Bayesian and likelihood approaches used in Phase I studies.
- Explain the continual reassessment method (CRM) and Fibonacci dose-ranging designs.
- Interpret measures of correlation and agreement.
- Recognize appropriate use of Pearson correlation, Spearman correlation, Kendall’s tau-b and Cohen’s Kappa statistics.
- Use a SAS program to produce confidence intervals for correlation coefficients and interpret the results.
- Adapt a SAS program to produce the correlation coefficients, their confidence intervals and Kendall’s tau-b.
- Recognize situations that call for the use of a statistic measuring concordance.
- Distinguish between a concordance correlation coefficient and a Kappa statistic based on the type of data used for each.
- Interpret a concordance correlation coefficient and a Kappa statistic.
References: 1 2 3 4
7.1 Considerations for Dose Finding Studies
The terms describing several types of early clinical studies are given below.
Definition 7.1 (Treatment Mechanism) Early developmental trial that investigates mechanism of treatment effect, e.g., a pharmacokinetics study of absorption and elimination of the drug from the human body
Definition 7.2 (Phase I) Imprecise term for dose-ranging studies
Definition 7.3 (Dose-escalation) Design or component of a design that specifies methods for increases in dose for subsequent subjects
Definition 7.4 (Dose-ranging) Design that tests some or all of a prespecified set of doses (fixed design points)
Definition 7.5 (Dose-finding) Design that titrates dose to a prespecified optimum based on biological or clinical considerations
Definitions from Piantodosi (2005)
Dose-finding (DF) trials are Phase I studies with the objective of determining the optimal biological dose (OBD) of a drug. In order to determine the dose with the highest potential for efficacy in the patient population that still meets safety criteria, dose-finding studies are typically conducted by administering sequentially rising doses to successive groups of individuals. Such studies may be conducted in healthy volunteers or in patients with the disease.
A question the investigator must answer in designing a dose-finding study is how to characterize an optimum dose. Should the optimum dose be selected on the basis of the highest therapeutic index (the maximal separation between risk and benefit)? Or is the optimal dose the level that maximizes therapeutic benefit while maintaining risk below a predetermined threshold? What measures will denote risk and benefit?
An optimal dose can be selected on the basis of efficacy alone, such as when a minimum effective dose (MED) is chosen for a pain-relieving medication and defined as the dose which eliminates mild-to-moderate pain in 80% of trial participants. In another case, the optimal dose might be selected as the highest dose that is associated with serious side effects in no more than 1 of 20 patients. This would be a maximum nontoxic dose (MND). In cancer therapeutics, the optimal dose for a cytotoxic drug designed to shrink tumors could be defined as the level that yields serious but reversible toxicity in no more than 30% of the patients. This is a maximum tolerated dose (MTD). Care in defining the conditions for optimality is critical to a dose-finding study.
Most DF trials are sequential studies such that the number of subjects is itself an outcome of the trial. Convincing evidence characterizing the relationship between dose and safety can be obtained after studying a small set of patients. Hence, sample size is not a major concern DF trials.
An idealized DF study would be similar to an animal bioassay design, with K fixed doses at increasing levels, d1, d2, … , dK. The hypothesized optimal dose would lie between d1 and dK. The n participants would be randomized to each of the K dose groups and the binary response of toxicity would be noted for each participant. A mathematical model could then be fit to the proportional responses over the doses such that the optimal dose could be determined.
Would you agree to participate in such a study? Think carefully …
Most likely, your answer is no because you would not want to risk being assigned to the highest dose level of this unproven drug as your first treatment. There is a principle here: it is unethical to treat humans at high doses of a drug without any prior knowledge of their responses at lower levels. Furthermore, ethics compel a design that minimizes the number of patients treated with both low ineffective doses and high toxic doses.
Thus, along with defining optimality, a DF study design usually includes a method for determining the starting dose for the patient, specification of dose increments and cohort sizes, definition of dose-limiting toxicities as well as the decision rules for escalation and de-escalation of the dose.
Continual Reassessment Method
The continual reassessment method (CRM) allows fitting a mathematical model to observed data during the study from which it estimates an optimal dose via extrapolation or interpolation. The next cohort of patients is assigned to the estimated optimal dose. A study using CRM would not have an a priori defined set of doses; thus is dose-finding study. The CRM itself can be thought of an algorithm for updating the best guess regarding the optimal dose. Bayesian approaches have also been incorporated into the CRM and the method is applicable for many types of responses.
In contrast to a dose-finding study, a dose-ranging study uses pre-specified design points.
Fibonacci Dose-Ranging Designs
Fibonacci, a thirteenth-century Italian mathematician, popularized the number sequence 1, 1, 2, 3, 5, 8, 13, 21, 34 … (a number in the sequence is the sum of the two previous numbers).
A fixed dosing scheme can be based on the Fibonacci sequence. For example, the first cohort of n participants is assigned dose D, the initial dose. If they tolerate this dose D well, the next cohort of n participants is assigned dose 2D. If all goes well with the second cohort, then a third cohort is assigned dose 3D, the fourth cohort is assigned does 5D, etc. The process is discontinued when one of the cohorts exhibits toxicity. Numerous modifications have been proposed to the Fibonacci scheme, such as allowing for de-escalation as well as escalation.
7.2 Dose-Finding (Phase I) Studies
One of the principal objectives of dose-finding (DF) studies is to assess the distribution and elimination of drug in the human system. What level of the drug is appropriate?
Pharmacokinetic (PK) models, also known as compartmental models, provide useful analytical approaches for DF studies. The objective of a PK model is to account for the absorption, distribution, metabolism, and excretion of a drug in the human system. An example of a two-compartment PK model is as follows and the objective is to estimate the rates into, between, and out of the two compartments.
Estimates are made for each of these areas, i.e., absorption rate, distribution rate, etc. for the drug in question.
7.3 Bayesian and Likelihood Approaches
Definition 7.6 (Data-dependent stopping) general term to describe any statistical or administrative reason for stopping a trial
Consideration of the reasons given earlier may lead you to stop the trial at an early stage, or at least change the protocol.
The review, interpretation, and decision-making aspects of clinical trials based on interim data are necessary but prone to error. If the investigators learn of interim results, then it could affect objectivity during the remainder of the trial, or if statistical tests are performed repeatedly on the accumulating data, then the Type I error rate is increased.
There is a natural conflict. On one hand, terminating the trial as early as possible will save costs and labor, expose as few patients as possible to inferior treatments, and allow disseminating information about the treatments quickly. On the other hand, there are pressures to continue the trial for as long as possible in order to increase precision, reduce errors of inference, obtain sufficient statistical power to account for prognostic factors and examine subgroups of interest, and gather information on secondary endpoints.
All the available statistical methods for interim analyses have some similar characteristics. They
- require the investigators to state the objectives clearly and in advance of the study,
- assume that the basic study design is sound,
- require some structure to the problem beyond the data being observed
- impose some penalty for early termination.
No statistical method is a substitute for judgment. The statistical criteria provide guidelines for terminating the trial because the decision to stop a trial is not based just on statistical information collected on one endpoint.
7.4 Likelihood Methods
It may be possible to assess treatment effects after each patient is accrued, treated, and evaluated. Such an approach is impractical in most circumstances, especially for trials that require lengthy follow-up to determine outcomes.
The first classical likelihood method proposed for this situation is called the sequential probability ratio test (SPRT) and it is based on the likelihood function. (This method is very rarely implemented because it is impractical in the clinical setting, but is important for historical reasons.) Let’s review this method in general terms here.
A likelihood function is constructed from a probability model for a sequence of random variables which correspond to the outcome measurements on the experimental units. In the likelihood function, however, the observed data points replace the random variables. Suppose we have a binary response (success/failure) from each patient which is determined immediately after a treatment is administered. (Again, not very practical.) However, for the situation discussed, we are examining one treatment which is administered to every patient. If there are N patients with K successes, and p represents the probability of success within each patient, then the likelihood function is based on the binomial probability function:
\[ L(p, K)=p^K(1-p)^{N-K} \]
This is a very simple likelihood function for a very simple example.
If the investigator is trying to decide whether \(p_0\) or \(p_1\) is the more appropriate value of p, then the likelihood ratio can be constructed to assess the evidence:
\[ R=\dfrac{L(p_0, K)}{L(p_1, K}=\left(\dfrac{p_0}{p_1} \right)^K \left(\dfrac{1-p_0}{1-p_1} \right)^{N-K} \]
This is a ratio of two different likelihood functions. If R is large, then the evidence is going to favor \(p_0\). If R is small, then the evidence is going to favor \(p_1\). Therefore, when analyzing interim data, we can calculate the likelihood ratio and stop the trial only if we have the amount of evidence that is expected for the target sample size.
Suppose that N is the target sample size and that after n patients there are k successes. After each treatment we will stop and analyze the data to determine whether to continue the trial or not. Under this scenario, we stop the trial if:
\[ R=\dfrac{L(p_0, K)}{L(p_1, K}=\left(\dfrac{p_0}{p_1} \right)^k \left(\dfrac{1-p_0}{1-p_1} \right)^{n-k} \le R_L \text{ or }\ge R_U \]
where \(R_L\) and \(R_U\) are prespecified constants. Let’s not worry about the details of the statistical calculation here. The values of \(R_L\) and \(R_U\) that correspond to testing \(H_0\colon p = p_0\) versus \(H_1 \colon p = p_1\) are \(R_L = \dfrac{\alpha}{(1 - \beta)}\) and \(R_U = \dfrac{(1 - beta)}{\alpha}\).
A sample schematic of the SPRT in practice is shown below. Here you would calculate R after the treatment of each patient. As you accumulate patients you can see that R is moving around as the trial proceeds. Before we had accrued all the patients that we wanted we hit the upper boundary and would not recruit the remaining patients.
Here is another example…
The SPRT might be useful in a phase II SE trial in which a treatment is to be monitored closely to determine if it reaches a certain level of success or failure. For example, suppose the investigator considers the treatment successful if \(p = 0.4\) (40% or greater), but considers it a failure if \(p = 0.2\) (20 % or less). Thus, the hypothesis testing problem is \(H_0 \colon p = 0.2\) vs. \(H_1 \colon p = 0.4\). Suppose we take \(\alpha = 0.05\) and \(\beta = 0.05\). Then the bounds would be calculated as \(R_L = \dfrac{1}{19}\) and \(R_U = 19\). We would reject \(H_0\) in favor of \(H_1\), and claim success, as soon as R gets small enough, \(R = (0.5)^k (1.33)^{n-k} \leq \dfrac{1}{19}\). On the other hand, we would stop the trial and accept \(H_0\) and reject \(H_1\), and claim failure, as soon as \(R \geq 19\).
The statistical formulation for the SPRT is relatively straightforward, but it is more commonly used in a quality control setting than in clinical trials. The obvious criticism is that each patient’s outcome must be observed quickly before you recruit the next patient. The SPRT also has the statistical property that it has a positive probability of never reaching the boundaries \(R_L\) and \(R_U\). If this is the case after the target sample size, N, is reached, then the trial is inconclusive.
7.5 Bayesian Methods
First, let’s review the Bayesian approach in general and then apply it to our current topic of likelihood methods.
The Bayesian approach to statistical design and inference is very different from the classical approach (the frequentist approach).
Before a trial begins, a Bayesian statistician summarizes the current knowledge or belief about the treatment effect, say we call it \(\theta\), in the form of a probability distribution. This is known as the prior distribution for \(\theta\). These assumptions are made prior to conducting the study and collecting any data.
Next, the data from the trial are observed, say we call it X, and the likelihood function of X given \(\theta\) is constructed. Finally, the posterior distribution for \(\theta\) given X is constructed. In essence, the prior distribution for \(\theta\) is revised into the posterior distribution based on the data X. The data collection in the study informs or revises the earlier assumptions.
The following schematic describes this Bayesian approach:
The development of the posterior distribution may be very difficult mathematically and it may be necessary to approximate it through computer algorithms.
The Bayesian statistician performs all inference for the treatment effect by formulating probability statements based on the posterior distribution. This is a very different approach and is not always accepted by the more traditional frequentist oriented statisticians.
In the Bayesian approach, \(\theta\) is regarded as a random variable, about which probability statements can be made. This is the appealing aspect of the Bayesian approach. In contrast, the frequentist approach regards \(\theta\) as a fixed but unknown quantity (called a parameter) that can be estimated from the data.
As an example of the contrasting philosophies, consider the frequentist description and the Bayesian description of a 95% confidence interval for \(\theta\).
Frequentist: “If a very large number of samples, each with the same sample size as the original sample, were taken from the same population as the original sample, and a 95% confidence interval constructed for each sample, then 95% of those confidence intervals would contain the true value of \(\theta\).” This is an extremely awkward and dissatisfying definition but technically represents the frequentist’s approach.
Bayesian: “The 95% confidence interval defines a region that covers 95% of the possible values of \(\theta\).” This is much more simple and straightforward. (As a matter of fact, most people when they first take a statistics course believe that this is the definition of a confidence interval.)
In a Bayesian analysis, if \(\theta\) is a parameter of interest, the analysis results in a probability distribution for \(\theta\). Using the probability distribution, many statements can be made. For example, if \(\theta\) represents a probability of success for a treatment, a statement can be made about the probability that \(\theta > 0.90\) (or any other value).
7.6 Bayesian approach in Clinical Trials
With respect to clinical trials, a Bayesian approach can cause some difficulties for investigators because they are not accustomed to representing their prior beliefs about a treatment effect in the form of a probability distribution. In addition, there may be very little prior knowledge about a new experimental therapy, so investigators may be reluctant to or not be able to quantify their prior beliefs. In the business world, the Bayesian approach is used quite often because of the availability of prior information. In the medical field, more often than not, this is not the case.
The choice of a prior distribution can be very controversial. Different investigators may select different priors for the same situation, which could lead to different conclusions about the trial. This is especially true when the data, X, are based on a small sample size because in such situations the prior distributions are modified only slightly to form the posterior distributions. Small sample sizes only modify the prior slightly. This tends to weight the posterior distribution very closely to the prior, therefore you are basing your results almost entirely on your prior assumptions.
When there is little prior information to base your assumptions of the distribution on, Bayesians employ a reference (or vague or non-informative) prior. These are intended to represent a minimal amount of prior information. Although vague priors may yield results similar to those of a frequentist approach, the priors may be unrealistic because they attempt to assign equal weight to all values of θ. Below you can see a very flat distribution, very spread out over a wide range of values.
Similarly, skeptical prior distributions are those that quantify the belief that large treatment effects are unlikely. Enthusiastic prior distributions are those that quantify large treatment effects. Let’s not worry about the calculations, but focus instead on the concepts here…
An example of a Bayesian approach for interim monitoring is as follows. Suppose an investigator plans a trial to detect a hazard ratio of 2 \(\left(\Lambda = 2\right)\) with 90% statistical power \(\left(\beta = 0.10\right)\) using at least a sample size of 90 events. The investigator plans one interim analysis, approximately halfway through trial, and a final analysis. (This is the more standard approach, as opposed to the SPRT where R was calculated after each treatment.)
The estimated logarithm of the hazard ratio is approximately normally distributed with variance \(\left(\dfrac{1}{d_1}\right) + \left(\dfrac{1}{d_2}\right)\), where \(d_1\) and \(d_2\) are the numbers of events in the two treatment groups. The null hypothesis is that the treatment groups are the same, i.e., \(H_0\colon \Lambda = 1\). Note that the \(log_e\) hazard ratio is 0 under the null hypothesis and the \(log_e\) hazard ratio is 0.693 when \(\Lambda = 2\), the proposed effect size.
Suppose the investigator has access to some pilot data or the published report of another investigator, in which there appeared to be a very small treatment effect with 16 events occurring within each of the two treatment groups. The investigator decides that this preliminary study will form the basis of a skeptical prior distribution for the \(log_e\) hazard ratio with a mean of 0 and a standard deviation of \(0.35 = {\dfrac{1}{16} + \dfrac{1}{16}}^{\frac{1}{2}}\). This is called a skeptical prior because it expresses skepticism that the treatment is beneficial.
Next, suppose that at the time of the interim analysis, (45 events have occurred), there are 31 events in one group and 14 events in the other group, such that the estimated hazard ratio is 2.25 (calculations not shown). These values are incorporated into the likelihood function, which modifies the prior distribution to yield the posterior distribution for the estimated \(log_e\) hazard ratio that has a mean = 0.474 and standard deviation = 0.228 (calculations not shown). Therefore, we can calculate the probability that \(\Lambda\) is \(> 2\). From the posterior distribution we construct the following probability statement:
\[ Pr[\Lambda \ge 2]=1-\Phi \left(\dfrac{log_e(2)-0.474}{0.228} \right)=1-\Phi(0.961)=0.168 \]
where \(\Phi\) represents the cumulative distribution function for the standard normal and is the true hazard ratio.
Conclusion: Based on the results from the interim analysis with a skeptical prior, there is not strong evidence that the treatment is effective because the posterior probability of the hazard ratio exceeding 2 is relatively small. Therefore, there is not enough evidence here to suggest that the study be stopped. What is too large? A reasonable value should be specified in your protocol before these values are determined.
In contrast, suppose that before the onset of the trial the investigator is very excited about the potential benefit of the treatment. Therefore, the investigator wants to use an enthusiastic prior for the \(log_e\) hazard ratio, i.e., a normal distribution with mean \(= log_e(2) = 0.693\) and standard deviation = 0.35 (same as the skeptical prior).
Suppose the interim data results are the same as those described above. This time, the posterior distribution for the \(log_e\) hazard ratio is normal with mean = 0.762 and standard deviation = 0.228. Then the probability for the posterior distribution is:
\[ Pr[\Lambda \ge 2]=1-\Phi \left(\dfrac{log_e(2)-0.762}{.228} \right)=1-\Phi(-0.302)=0.682 \]
This is a drastic change in the probability based on the assumptions that were made ahead of time. In this case, the investigator still may not consider this to be strong evidence that the trial should terminate because the posterior probability of the hazard ratio exceeding 2 does not exceed 0.90.
Nevertheless, the example demonstrates the controversy that can arise with a Bayesian analysis when the amount of experimental data is small, i.e., the selection of the prior distribution drives the decision-making process. For this reason, many investigators prefer to use non-informative priors. Using the Bayesian methods, you can make probability statements about your expected results.
7.7 Correlation and Agreement
Overview
Many biostatistical analyses are conducted to study the relationship between two continuous or ordinal scale variables within a group of patients.
Purposes of these analyses include:
- assessing correlation between the two variables, i.e., identifying whether values of one variable tend to be higher (or possibly lower) for higher values of the other variable;
- assessing the amount of agreement between the values of the two variables, i.e., comparing alternative ways of measuring or assessing the same response;
- assessing the ability of one variable to predict values of the other variable, i.e., formulating predictive models via regression analyses.
This lesson will focus only on correlation and agreement, (issues numbered 1 and 2 listed above).
7.8 Pearson Correlation Coefficient
Correlation is a general method of analysis useful when studying possible association between two continuous or ordinal scale variables. Several measures of correlation exist. The appropriate type for a particular situation depends on the distribution and measurement scale of the data. Three measures of correlation are commonly applied in biostatistics and these will be discussed below.
Suppose that we have two variables of interest, denoted as X and Y, and suppose that we have a bivariate sample of size n:
\[ \left(X_{1} , Y_{1} \right), \left(X_{2} , Y_{2} \right), \dots , \left(X_{n} , Y_{n} \right) \]
and we define the following statistics:
\[ \bar{X}=\dfrac{1}{n}\sum_{i=1}^{n}X_i , S_{XX}=\dfrac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar{X})^2 \]
\[ \bar{Y}=\dfrac{1}{n}\sum_{i=1}^{n}Y_i , S_{YY}=\dfrac{1}{n-1}\sum_{i=1}^{n}(Y_i-\bar{Y})^2 \]
\[ S_{XY}=\dfrac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar{X})(Y_i-\bar{Y}) \]
These statistics above represent the sample mean for X, the sample variance for X, the sample mean for Y, the sample variance for Y, and the sample covariance between X and Y, respectively. These should be very familiar to you.
The sample Pearson correlation coefficient (also called the sample product-moment correlation coefficient) for measuring the association between variables X and Y is given by the following formula:
\[ r_p=\dfrac{S_{XY}}{\sqrt{S_{XX}S_{YY}}} \]
The sample Pearson correlation coefficient, \(r_{p}\) , is the point estimate of the population Pearson correlation coefficient
\[ \rho_p=\dfrac{\sigma_{XY}}{\sqrt{\sigma_{XX}\sigma_{YY}}} \]
The Pearson correlation coefficient measures the degree of linear relationship between X and Y and \(-1 ≤ r_{p} ≤ +1\), so that \(r_{p}\) is a “unitless” quantity, i.e., when you construct the correlation coefficient the units of measurement that are used cancel out. A value of +1 reflects perfect positive correlation and a value of -1 reflects perfect negative correlation.
For the Pearson correlation coefficient, we assume that both X and Y are measured on a continuous scale and that each is approximately normally distributed.
The Pearson correlation coefficient is invariant to location and scale transformations. This means that if every \(X_{i}\) is transformed to
\[ X_{i} * = aX_{i} + b \]
and every \(Y_{i}\) is transformed to
\[ Y_{i} * = cY_{i} + d \]
where \(a > 0, b, c > 0\), and d are constants, then the correlation between X and Y is the same as the correlation between \(X^*\) and \(Y^*\).
With SAS, PROC CORR is used to calculate \(r_{p}\). The output from PROC CORR includes summary statistics for both variables and the computed value of \(r_{p}\). The output also contains a p-value corresponding to the test of:
\[ H_{0} : \rho_{p} = 0 \text{ versus } H_{1} : \rho_{p} ≠ 0 \]
It should be noted that this statistical test generally is not very useful, and the associated p-value, therefore, should not be emphasized. What is more important is to construct a confidence interval.
The sampling distribution for Pearson’s \(r_{p}\) is not normal. In order to attain confidence limits for \(r_{p}\) based on a standard normal distribution, we transform \(r_{p}\) using Fisher’s Z transformation to get a quantity, \(z_{p}\), that has an approximate normal distribution. Then we can work with this value. Here is what is involved in the transformation.
Fisher’s Z transformation is defined as
\[ z_p=\dfrac{1}{2}log_e\left( \dfrac{1+r_p}{1-r_p} \right) \sim N\left( \zeta_p , sd=\dfrac{1}{\sqrt{n-3}} \right) \]
where
\[ \zeta_p=\dfrac{1}{2}log_e\left( \dfrac{1+\rho_p}{1-\rho_p} \right) \]
We will use this to get the usual confidence interval, so, an approximate \(100(1 - \alpha)\%\) confidence interval for \(\zeta_{p}\) is given by \(z_{p, \frac{\alpha}{2}} , z_{p, 1-\frac{\alpha}{2}}\), where
\[ z_{p , \alpha/2}=z_p-\left( t_{n-3 , 1-\alpha/2}/\sqrt{n-3} \right) , z_{p , 1-\alpha/2}=z_p+\left( t_{n-3 , 1-\alpha/2}/\sqrt{n-3} \right) \]
But really what we want is an approximate \(100(1 - \alpha)\%\) confidence interval for \(\rho_{p}\) is given by \(r_{p, \frac{\alpha}{2}} , r_{p, 1- \frac{\alpha}{2}}\), where
\[ r_{p , \frac{\alpha}{2}}=\dfrac{exp(2z_{p , \alpha/2})-1}{exp(2z_{p , \alpha/2})+1},r_{p , 1-\alpha/2}=\dfrac{exp(2z_{p , 1-\alpha/2})-1}{exp(2z_{p , 1-\alpha/2})+1} \]
Again, you do not have to do this by hand. PROC CORR in SAS will do this for you, but it is important to have an idea of what is going on.
7.9 Spearman Correlation Coefficient
The Spearman rank correlation coefficient, \(r_s\), is a nonparametric measure of correlation based on data ranks. It is obtained by ranking the values of the two variables (X and Y) and calculating the Pearson \(r_p\) on the resulting ranks, not the data itself. Again, PROC CORR will do all of these actual calculations for you.
The Spearman rank correlation coefficient has properties similar to those of the Pearson correlation coefficient, although the Spearman rank correlation coefficient quantifies the degree of linear association between the ranks of X and the ranks of Y. Also, \(r_s\) does not estimate a natural population parameter (unlike Pearson’s \(r_p\) which estimates \(\rho_p\) ).
An advantage of the Spearman rank correlation coefficient is that the X and Y values can be continuous or ordinal, and approximate normal distributions for X and Y are not required. Similar to the Pearson \(r_p\), Fisher’s Z transformation can be applied to the Spearman \(r_s\) to get a statistic, \(z_s\), that has an asymptotic normal distribution for calculating an asymptotic confidence interval. Again, PROC CORR will do this as well.
7.10 Kendall Tau-b Correlation Coefficient
The Kendall tau-b correlation coefficient, \(\tau_b\), is a nonparametric measure of association based on the number of concordances and discordances in paired observations.
Suppose two observations \(\left(X_i , Y_i \right)\) and \(\left(X_j , Y_j \right)\) are concordant if they are in the same order with respect to each variable. That is, if
- \(X_i < X_j\) and \(Y_i < Y_j\) , or if
- \(X_i > X_j\) and \(Y_i > Y_j\)
They are discordant if they are in the reverse ordering for X and Y, or the values are arranged in opposite directions. That is, if
- \(X_i < X_j\) and \(Y_i > Y_j\) , or if
- \(X_i > X_j\) and \(Y_i < Y_j\)
The two observations are tied if \(X_i = X_j\) and/or \(Y_i = Y_j\) .
The total number of pairs that can be constructed for a sample size of n is
\[ N=\binom{n}{2}=\dfrac{1}{2}n(n-1) \]
N can be decomposed into these five quantities:
\[ N = P + Q + X_0 + Y_0 + (XY)_0 \]
where P is the number of concordant pairs, Q is the number of discordant pairs, \(X_0\) is the number of pairs tied only on the X variable, \(Y_0\) is the number of pairs tied only on the Y variable, and \(\left(XY\right)_0\) is the number of pairs tied on both X and Y.
The Kendall tau-b for measuring order association between variables X and Y is given by the following formula:
\[ t_b=\dfrac{P-Q}{\sqrt{(P+Q+X_0)(P+Q+Y_0)}} \]
This value becomes scaled and ranges between -1 and +1. Unlike Spearman it does estimate a population variance as:
\[ t_b \text{ is the sample estimate of } t_b = Pr[\text{concordance}] - Pr[\text{discordance}] \]
The Kendall tau-b has properties similar to the properties of the Spearman \(r_s\). Because the sample estimate, \(t_b\) does estimate a population parameter, \(t_b\), many statisticians prefer the Kendall tau-b to the Spearman rank correlation coefficient.
Example 7.1 (Correlation Coefficients) This example provides an IML module for calculating point and interval estimates of the Pearson correlation coefficient and the concordance correlation coefficient.
Download the SAS code: sas_correlation.sas
Age and percentage body fat were measured in 18 adults. SAS PROC CORR provides estimates of the Pearson, Spearman, and Kendall correlation coefficients. It also calculates Fisher’s Z transformation for the Pearson and Spearman correlation coefficients in order to get 95% confidence intervals.
*******************************************************************************
* This program indicates how to construct a bivariate scatterplot with an *
* overlay of the least squares regression line. *
* *
* This program also provides an example for calculating point and *
* interval estimates of the Pearson, Spearman, and Kendall correlation *
* coefficients. *
*******************************************************************************;
data bodyfat;
input subject age bodyfat_perc;
cards;
01 23 9.5
02 23 27.9
03 27 7.8
04 27 17.8
05 39 31.4
06 41 25.9
07 45 27.4
08 49 25.2
09 50 31.1
10 53 34.7
11 53 42.0
12 54 29.1
13 56 32.5
14 57 30.3
15 58 33.0
16 58 33.8
17 60 41.1
18 61 34.5
;
run;
proc gplot data=bodyfat;
plot bodyfat_perc*age/vaxis=axis1 haxis=axis2 nolegend frame;
axis1 label=(a=90 '% Body Fat') minor=none;
axis2 label=('Age') minor=none;
symbol1 value=star color=black interpol=r;
title "Scatterplot";
run;
proc corr data=bodyfat Pearson Spearman Kendall Fisher(biasadj=no);
var age;
with bodyfat_perc;
title "Correlation Coefficients";
run;The resulting estimates for this example are 0.7921, 0.7539, and 0.5762, respectively for the Pearson, Spearman, and Kendall correlation coefficients. The Kendall tau-b correlation typically is smaller in magnitude than the Pearson and Spearman correlation coefficients.
The 95% confidence intervals are (0.5161, 0.9191) and (0.4429, 0.9029), respectively for the Pearson and Spearman correlation coefficients. Because the Kendall correlation typically is applied to binary or ordinal data, its 95% confidence interval can be calculated via SAS PROC FREQ (this is not shown in the SAS program above).
7.11 Use and Misuse of Correlation Coefficients
Correlation is a widely-used analysis tool that sometimes is applied inappropriately. Some caveats regarding the use of correlation methods follow.
- The correlation methods discussed in this chapter should be used only with independent data; they should not be applied to repeated measures data where the data are not independent. For example, it would not be appropriate to use these measures of correlation to describe the relationship between Week 4 and Week 8 blood pressures in the same patients.
- Caution should be used in interpreting results of correlation analysis when large numbers of variables have been examined, resulting in a large number of correlation coefficients.
- The correlation of two variables that both have been recorded repeatedly over time can be misleading and spurious. Time trends should be removed from such data before attempting to measure correlation.
- To extend correlation results to a given population, the subjects under study must form a representative (i.e., random) sample from that population. The Pearson correlation coefficient can be very sensitive to outlying observations and all correlation coefficients are susceptible to sample selection biases.
- Care should be taken when attempting to correlate two variables where one is a part and one represents the total. For example, we would expect to find a positive correlation between height at age ten and adult height because the second quantity “contains” the first quantity.
- Correlation should not be used to study the relation between an initial measurement, X, and the change in that measurement over time, Y - X. X will be correlated with Y - X due to the regression to the mean phenomenon.
- Small correlation values do not necessarily indicate that two variables are unassociated. For example, Pearson’s \(r_p\) will underestimate the association between two variables that show a quadratic relationship. Scatterplots should always be examined.
- Correlation does not imply causation. If a strong correlation is observed between two variables A and B, there are several possible explanations:
- A influences B
- B influences A
- A and B are influenced by one or more additional variables
- the relationship observed between A and B was a chance error.
- “Regular” correlation coefficients are often published when the researcher really intends to compare two methods of measuring the same quantity with respect to their agreement. This is a misguided analysis because correlation measures only the degree of association; it does not measure agreement. The next section of this lesson will present a measure of agreement.
7.12 Concordance Correlation Coefficient for Measuring Agreement
How well do two diagnostic measurements agree? Many times continuous units of measurement are used in the diagnostic test. We may not be interested in correlation or linear relationship between the two measures, but in a measure of agreement.
The concordance correlation coefficient, \(r_c\) , for measuring agreement between continuous variables X and Y (both approximately normally distributed), is calculated as follows:
\[ r_c=\dfrac{2S_{XY}}{S_{XX}+S_{YY}+(\bar{X}-\bar{Y})^2} \]
Similar to the other correlation coefficient, the concordance correlation satisfies \(-1 ≤ r_c ≤ +1\). A value of \(r_c = +1\) corresponds to perfect agreement. A value of \(r_c = - 1\) corresponds to perfect negative agreement, and a value of \(r_c = 0\) corresponds to no agreement. The sample estimate, \(r_c\) , is an estimate of the population concordance correlation coefficient:
\[ \rho_c=\dfrac{2\sigma_{XY}}{\sigma_{XX}+\sigma_{YY}+(\mu_{X}-\mu_{Y})^2} \]
Let’s look at an example that will help to make this concept clearer.
Example 7.2 (SAS PROC FREQ option for constructing Cohen’s kappa and weighted kappa statistics) Download the SAS code: sas_agreement_concordanc.sas:
*******************************************************************************
* This program indicates how to construct a bivariate scatterplot with an *
* overlay of the line of identity. *
* *
* This program also provides an IML module for calculating point and *
* interval estimates of the Pearson correlation coefficient and the *
* concordance correlation coefficient. *
*******************************************************************************;
data dice_baseline;
input subject cort_auc1 cort_auc2;
cards;
61001 5.28170 5.37851
61002 5.58796 5.33628
61003 6.47607 6.59770
61005 6.36019 6.39746
61007 5.81121 5.82528
61008 6.03036 6.21147
61009 5.84549 6.10434
61014 6.80349 6.90689
61015 6.28977 6.27369
61023 5.88446 5.94352
61024 5.79701 6.04876
61026 5.01302 4.72154
61027 6.48824 6.37891
61028 5.30862 5.53405
61033 5.70905 5.77803
61034 5.98545 5.77613
61035 6.19924 6.20880
61036 6.00639 5.98313
61037 6.27793 6.44342
61038 6.57390 6.62936
61039 5.69639 5.43509
61042 5.96588 6.01282
61044 6.04803 6.11529
61046 6.39423 6.03876
62002 7.44584 7.58421
62003 5.90813 5.87230
62004 6.05483 5.94695
62005 5.65735 5.64983
62006 6.44815 6.44280
62007 6.28611 6.45374
62009 6.40863 6.14994
62012 5.62564 5.58142
62013 6.68375 6.58815
62014 5.76951 5.84802
62015 5.94383 5.95489
62016 5.66024 5.73711
62017 4.77492 4.57465
62018 5.60468 5.43495
62019 6.82819 6.86652
62020 5.18986 5.05725
62021 6.48810 6.59655
62022 6.08867 5.76965
62023 5.91400 5.89672
62024 5.58217 5.57651
62026 6.32857 6.48921
62027 7.67703 7.76541
62028 5.92411 5.66689
62029 6.15313 6.16558
62030 5.20392 5.42481
62032 6.43962 6.46171
62033 6.20661 6.28542
63001 6.04767 5.82528
63002 6.46923 6.51728
63003 5.68370 5.79701
63004 5.11719 5.47363
63005 6.10993 6.13541
63006 4.91744 5.06968
63008 5.35972 5.56605
63010 6.73016 6.76285
63011 5.93700 5.94092
63012 6.07716 5.92548
63014 6.58185 6.52781
63015 5.84317 5.85030
63016 5.98144 6.22389
63017 5.77452 5.81662
63018 5.46142 5.84898
63021 5.44920 5.43688
63022 5.96519 6.01302
63026 6.29258 6.42339
63027 6.86608 6.92806
63028 5.47875 5.72634
63030 6.16190 6.16608
63032 5.66707 5.97114
63033 5.80634 5.63640
63034 6.37256 6.24416
63035 5.65755 5.98070
64001 6.07596 6.06010
64002 6.57898 6.54552
64003 6.60733 6.89724
64004 5.69611 5.82963
64005 6.60331 6.62972
64007 5.89963 5.83195
64008 6.19731 6.07044
64009 5.88875 6.03427
64010 5.64912 5.46135
64011 6.99962 7.10324
64012 6.61282 6.68520
64015 6.60477 6.76468
64016 5.35161 5.55307
64017 6.63249 6.77868
64020 5.37717 5.19760
64023 5.58781 5.87044
64025 5.74499 5.77090
64027 5.92655 6.09265
64028 5.01060 5.17112
64030 7.12400 7.17368
64031 5.79909 5.37603
64032 5.75609 5.85174
64033 6.79275 6.78255
64035 6.02198 6.03915
64036 5.43960 5.82229
64037 5.20163 5.12928
65001 6.18838 6.66593
65002 6.13860 6.26224
65003 6.98807 6.97658
65005 5.54628 5.50547
65006 4.47249 4.59256
65007 5.04034 5.07775
65008 5.42025 5.46227
65009 5.26772 5.41463
65011 6.43019 6.38438
65013 6.56323 6.46607
65014 5.06134 4.70619
65016 6.32666 6.31564
65017 5.90235 6.05890
65018 6.05800 5.99467
65020 5.96388 6.01204
65021 5.57324 5.61324
65022 6.21017 6.15262
65025 6.01934 6.00318
65026 5.52082 5.54575
65027 5.89237 5.67469
65028 6.24592 6.37106
65031 6.31524 6.44334
65032 6.19602 6.29576
65033 6.07305 6.07966
65037 5.60960 5.38492
65038 5.39806 5.18466
65040 5.95464 6.20802
66003 5.96880 5.86128
66004 5.89707 5.69116
66005 6.27067 6.35294
66006 5.39744 5.23236
66010 6.64408 6.64990
66011 6.29430 6.32724
66012 5.22507 5.28754
66013 5.15840 5.05580
66014 5.85385 5.54974
66018 6.49611 5.87795
66019 5.43285 5.50871
66021 6.32416 6.31612
66023 5.77253 5.74469
66024 5.89920 5.95774
;
run;
proc means data=dice_baseline noprint;
var cort_auc1 cort_auc2;
output out=dice_baselinemin min=var1 var2;
run;
proc means data=dice_baseline noprint;
var cort_auc1 cort_auc2;
output out=dice_baselinemax max=var1 var2;
run;
data dice_baselineall;
set dice_baselinemin dice_baselinemax;
drop _type_ _freq_;
if _N_=1 then var1=floor(min(var1,var2));
if _N_=1 then var2=floor(min(var1,var2));
if _N_=2 then var1=ceil(max(var1,var2));
if _N_=2 then var2=ceil(max(var1,var2));
run;
data dice_baselineall;
set dice_baseline dice_baselineall;
run;
proc gplot uniform data=dice_baselineall;
plot cort_auc1*cort_auc2 var1*var2/overlay vaxis=axis1 haxis=axis2 nolegend frame;
axis1 label=(a=90 'Cortisol Every Hour') minor=none;
axis2 label=('Cortisol Every Two Hours') minor=none;
symbol1 value=star color=black interpol=none;
symbol2 value=none color=black interpol=join;
title "Concordance Correlation Coefficient";
run;
proc iml;
*******************************************************************************
* Enter the appropriate SAS data set name in the use statement and enter the *
* appropriate variable names in the read statements. *
*******************************************************************************;
use dice_baseline;
read all var {cort_auc1} into var1;
read all var {cort_auc2} into var2;
*******************************************************************************
* The IML module, labeled concorr, starts next. *
*******************************************************************************;
start concorr;
nonmiss=loc(var1#var2^=.);
var1=var1[nonmiss];
var2=var2[nonmiss];
free nonmiss;
n=nrow(var1);
mu1=sum(var1)/n;
mu1=round(mu1,0.0001);
mu2=sum(var2)/n;
mu2=round(mu2,0.0001);
sigma11=ssq(var1-mu1)/(n-1);
sigma11=round(sigma11,0.0001);
sigma22=ssq(var2-mu2)/(n-1);
sigma22=round(sigma22,0.0001);
sigma12=sum((var1-mu1)#(var2-mu2))/(n-1);
sigma12=round(sigma12,0.0001);
lshift=(mu1-mu2)/((sigma11#sigma22)##0.25);
rho=sigma12/sqrt(sigma11#sigma22);
rho=round(rho,0.0001);
z=log((1+rho)/(1-rho))/2;
se_z=1/sqrt(n-3);
t=tinv(0.975,n-3);
z_low=z-(se_z#t);
z_upp=z+(se_z#t);
rho_low=(exp(2#z_low)-1)/(exp(2#z_low)+1);
rho_low=round(rho_low,0.0001);
rho_upp=(exp(2#z_upp)-1)/(exp(2#z_upp)+1);
rho_upp=round(rho_upp,0.0001);
crho=(2#sigma12)/((sigma11+sigma22)+((mu1-mu2)##2));
crho=round(crho,0.0001);
z=log((1+crho)/(1-crho))/2;
if sigma12^=0 then do;
t1=((1-(rho##2))#(crho##2))/((1-(crho##2))#(rho##2));
t2=(2#(crho##3)#(1-crho)#(lshift##2))/(rho#((1-(crho##2))##2));
t3=((crho##4)#(lshift##4))/(2#(rho##2)#((1-(crho##2))##2));
se_z=sqrt((t1+t2-t3)/(n-2));
end;
else se_z=sqrt(2#sigma11#sigma22)/((sigma11+sigma22+((mu1-mu2)##2))#(n-2));
t=tinv(0.975,n-2);
z_low=z-(se_z#t);
z_upp=z+(se_z#t);
crho_low=(exp(2#z_low)-1)/(exp(2#z_low)+1);
crho_low=round(crho_low,0.0001);
crho_upp=(exp(2#z_upp)-1)/(exp(2#z_upp)+1);
crho_upp=round(crho_upp,0.0001);
Results=n//mu1//mu2//sigma11//sigma22//sigma12//rho_low//rho//rho_upp//
crho_low//crho//crho_upp;
r_name={'SampleSize' 'Mean_1' 'Mean_2' 'Variance_1' 'Variance_2' 'Covariance' 'Corr LowerCL'
'Corr' 'Corr UpperCL' 'ConcCorr LowerCL' 'ConcCorr' 'ConcCorr UpperCL'};
print 'The Estimated Correlation and Concordance Correlation (and 95% Confidence Limits)';
print Results [rowname=r_name];
finish concorr;
*******************************************************************************
* The IML module, labeled concorr, is finished. *
*******************************************************************************;
run concorr;The ACRN DICE trial was discussed earlier in this course. In that trial, participants underwent hourly blood draws between 08:00 PM and 08:00 AM once a week in order to determine the cortisol area-under-the-curve (AUC). The participants hated this! They complained about the sleep disruption every hour when the nurses came by to draw blood, so the ACRN wanted to determine for future studies if the cortisol AUC calculated on measurements every two hours was in good agreement with the cortisol AUC calculated on hourly measurements. The baseline data were used to investigate how well these two measurements agreed. If there is good agreement, the protocol could be changed to take blood every two hours.
The SAS program yielded \(r_c = 0.95\) and a 95% confidence interval = (0.93, 0.96). The ACRN judged this to be excellent agreement, so it will use two-hourly measurements in future studies.
What about binary or ordinal data? Cohen’s Kappa Statistic will handle this…
7.13 Cohen’s Kappa Statistic for Measuring Agreement
Cohen’s kappa statistic, \(\kappa\) , is a measure of agreement between categorical variables X and Y. For example, kappa can be used to compare the ability of different raters to classify subjects into one of several groups. Kappa also can be used to assess the agreement between alternative methods of categorical assessment when new techniques are under study.
Kappa is calculated from the observed and expected frequencies on the diagonal of a square contingency table. Suppose that there are n subjects on whom X and Y are measured, and suppose that there are g distinct categorical outcomes for both X and Y. Let \(f_{ij}\) denote the frequency of the number of subjects with the \(i^{th}\) categorical response for variable X and the \(j^{th}\) categorical response for variable Y.
Then the frequencies can be arranged in the following g × g table:
| Y = 1 | Y = 2 | … | Y = g | |
|---|---|---|---|---|
| X = 1 | \(f_{11}\) | \(f_{12}\) | … | \(f_{1g}\) |
| X = 2 | \(f_{21}\) | \(f_{22}\) | … | \(f_{2g}\) |
| … | … | … | … | …… |
| X = g | \(f_{g1}\) | \(f_{g2}\) | … | \(f_{gg}\) |
The observed proportional agreement between X and Y is defined as:
\[ p_0=\dfrac{1}{n}\sum_{i=1}^{g}f_{ii} \]
and the expected agreement by chance is:
\[ p_e=\dfrac{1}{n^2}\sum_{i=1}^{g}f_{i+}f_{+i} \]
where \(f_{i+}\) is the total for the \(i^{th}\) row and \(f_{+i}\) is the total for the \(i^{th}\) column. The kappa statistic is:
\[ \hat{\kappa}=\dfrac{p_0-p_e}{1-p_e} \]
Cohen’s kappa statistic is an estimate of the population coefficient:
\[ \kappa=\dfrac{Pr[X=Y]-Pr[X=Y|X \text{ and }Y \text{ independent}]}{1-Pr[X=Y|X \text{ and }Y \text{ independent}]} \]
Generally, \(0 ≤ \kappa ≤ 1\), although negative values do occur on occasion. Cohen’s kappa is ideally suited for nominal (non-ordinal) categories. Weighted kappa can be calculated for tables with ordinal categories.
Example 7.3 (Multiple Regression) Download the SAS code: sas_agreement_Cohen.sas
Two radiologists rated 85 patients with respect to liver lesions. The ratings were designated on an ordinal scale as:
0 =‘Normal’
1 =‘Benign’
2 =‘Suspected’
3 =‘Cancer’
SAS PROC FREQ provides an option for constructing Cohen’s kappa and weighted kappa statistics.
*******************************************************************************
* This program indicates how to calculate Cohen's kappa statistic for *
* evaluating the level of agreement between two variables. *
*******************************************************************************;
proc format;
value raterfmt 0='Normal' 1='Benign' 2='Suspected' 3='Cancer';
run;
data radiology;
input rater1 rater2 count;
format rater1 rater2 raterfmt.;
cards;
0 0 21
0 1 12
0 2 0
0 3 0
1 0 4
1 1 17
1 2 1
1 3 0
2 0 3
2 1 9
2 2 15
2 3 2
3 0 0
3 1 0
3 2 0
3 3 1
;
run;
proc freq data=radiology;
tables rater1*rater2/agree;
weight count;
test kappa;
exact kappa;
title "Cohen's Kappa Coefficients";
run;The weighted kappa coefficient is 0.57 and the asymptotic 95% confidence interval is (0.44, 0.70). This indicates that the amount of agreement between the two radiologists is modest (and not as strong as the researchers had hoped it would be).
7.14 Summary
In this lesson, among other things, we learned:
- To state the components of a typical dose-finding design.
- To state the objectives of a pharmacokinetic model.
- To describe Bayesian and likelihood approaches used in Phase I studies.
- To explain the continual reassessment method (CRM) and Fibonacci dose-ranging designs.
- To interpret measures of correlation and agreement.
- recognize appropriate use of Pearson correlation, Spearman correlation, Kendall’s tau-b and Cohen’s Kappa statistics.
- use a SAS program to produce confidence intervals for correlation coefficients and interpret the results.
- adapt a SAS program to produce the correlation coefficients, their confidence intervals and Kendall’s tau-b.
- recognize situations that call for the use of a statistic measuring concordance.
- distinguish between a concordance correlation coefficient and a Kappa statistic based on the type of data used for each.
- interpret a concordance correlation coefficient and a Kappa statistic.
DeMets DL, Lan KK, 1994, Interim analysis: The alpha spending function approach, Statistics in Medicine 13: 1341-1352.↩︎
Piantadosi Steven. (2005) Clinical trials as research, Why clinical trials are ethical, Contexts for clinical trials. In: Piantadosi Steven. Clinical Trials: A Methodologic Perspective. 2nd ed. Hobaken, NJ: John Wiley and Sons, Inc.↩︎
Ellenberg, SS. Fleming, TR. DeMets, DL. 2002, Data Monitoring Committees in Clinical Trials, New York, NY: Wiley.↩︎
Pocock, S.J. 1983 Clinical Trials: A Practical Approach. Chichester: John Wiley and Sons.↩︎