## Statistical Power and Confidence Regions

The significance level is the determining factor in the specification of the rejection region of a statistical test. Only the distribution under the null assumption of no signal plays a role in setting the level of the threshold, once the test statistic and the general form of test are decided upon. However, after setting that threshold, one can examine other statistical properties of the resulting test. A central property is statistical power of the test - the probability to reject the null hypothesis when a signal is present. Since this probability depends on the values of the parameters, one often speaks of the power function to emphasize this dependence. For a test at a single marker, this probability is obtained approximately from the normal distribution; it is a function of the noncentrality parameter given by (4.6). In this chapter we will examine the concept of power for a whole-genome scan.

The primary interest is now focused on the case where the null hypothesis is false. Statisticians define the power of an hypothesis test as the probability of concluding correctly the falsity of the null hypothesis. However, the case of a genome scan is more subtle than a simple test of hypothesis. There exists the possibility that due to random fluctuations the significance threshold is exceeded on a chromosome that does not actually contain a QTL. Unless the threshold is also crossed on a chromosome containing a QTL, one would correctly conclude that the simple null hypothesis of no QTL anywhere in the genome is false, but would identify the chromosomal location of the QTL incorrectly. We are particularly interested in the probability of exceeding the threshold on a chromosome containing a QTL, or perhaps even at some marker close to the QTL, say within 20 cM. Although any definition of power in this context is somewhat arbitrary, in this book the power to detect a particular QTL will refer to the probability of correctly identifying the chromosome inhabited by the QTL. This means that when there is more than one QTL, power refers to specific QTL and can vary from one to another, depending on the effect of the QTL on the phenotype. In the case that multiple QTL lie on the same chromosome, one might want to make more subtle distinctions.

While keeping this possibility in mind, we shall for the most part ignore it in our statistical analysis.

At one end of the spectrum, when there is no QTL or only very weak QTL on a given chromosome, the power function is essentially equal to the (chromosome specific) significance level. At the other extreme are parameter values which correspond to a signal that is so large as to make the power approximately equal to one. We will mainly be interested in interim parameter values, values for which the power function is in the range 50%-95%.

The main application of the power function is to help us choose an experimental design - especially the breeding design, marker density, and sample size. The significance level is set to be some fixed value - typically 5% - regardless of the design. The differences between designs will be reflected in their power to detect the signal. Thus, for example, a sample size which is too small, or a collection of markers which is not dense enough, may compromise the chances of successful detection of a QTL. On the other hand, it is neither economically efficient nor ethical to use more animals then needed. Moreover, since genotyping is an expensive component in a genome scan, using more markers than needed is a waste of time and money that can be used for other purposes. Careful planning of an experiment can ensure efficient distribution of resources, without a substantial reduction in power. In the body of the text we focus on issues of sample size determination and the selection of the inter-marker spacing in the context of the backcross design. The power of other breeding designs is left to exercises in the problem set at the end of the chapter.

In the first section we identify the terms that affect the power of detection in a whole-genome scan. In the second section we introduce analytic formulas for the power. These formulas, like the formulas used for the computation of the significance level, are given in the context of the Ornstein-Uhlenbeck process. They will allow us to analyze the power function and to examine the effect of changing the values of various parameters. Consequently, in the third section we apply these formulas to select a good experimental design for the detection of a QTL. In the last two sections we deal with issues related more to estimation. In the fourth section we consider the construction of confidence intervals for the location of a QTL and in the last section the construction of a lower confidence bound for the effect of the QTL.

### 6.1 The Power to Detect a QTL

In this section we identify the parameters that determine the statistical properties of the monitoring process in the presence of a QTL. The examination is carried out in the context of local alternatives. In the case of a single marker, which was analyzed in Chap. 4, this corresponded to considering a shifted normal distribution. Similarly, for the multi-marker process, where the null distribution corresponded to the Gaussian Ornstein-Uhlenbeck process, the computation under the alternative will involve the same process with a shifted mean function. Since we deal with local alternatives, the correlation structure is not affected. The power to detect a signal is the probability that the maximum absolute value of the process, with the shifted mean, exceeds the threshold.

Recall that for the single marker process the expected value under the alternative of the test statistic at a marker equals its expectation at the QTL multiplied by the correlation between the marker and the QTL (see 4.6). For the backcross design and the simple model of QTL, we use the expectation at the QTL itself, which equals £ = (a0 + S0)/(2ay). Here a0 corresponds to the (local) parameter of additive effect, S0 is the (local) parameter of dominance effect, and ay is the standard deviation of the phenotype. In terms of the original parameters of the model, a0 = n12a and So = n1/2S. The correlation between marker and QTL for a backcross design, under the Haldane model of recombination, is equal to 1 — 2 0 = exp{—0.02 \t — t|}. Here \t — t\ corresponds to the distance between a QTL located t cM from the telomere, and a marker located t cM from the telomere. This information is summarized by the formula

E(Zt) = a0±S0(1 — 20) = £exp(—0.02 \t — t\) . (6.1)

2 ay

If the marker and the QTL are not on the same chromosome, then 0 = 1/2, the genetic distance from the QTL is defined to be infinite, and the expectation is equal to 0.

Equation (6.1) gives a complete description of the mean function of the multi-marker process under the model of a single QTL. The QTL is located on some chromosome, t cM from the telomere. The mean function for the markers on the same chromosome is a function of their distance from the QTL. It decreases exponentially fast, on both sides of the QTL, as the distance from the QTL increases. The mean function of the multi-marker process over the other 19 chromosomes, which do not contain the QTL, is identically equal to 0. The covariance structure of the process under the alternative is identical to the covariance structure under the null assumption. Thus, the standard deviation of the test statistics Zt is equal to one, regardless of the location of the marker and its distance from the QTL. The correlation between any pair of markers is a function of the genetic distance between them. If the two markers are located on different chromosomes, then the genetic distance between them is infinite, and the two markers are uncorrelated. If the markers are located on the same chromosome, s and t cM from the telomere, then the correlation between them equals exp{—0.02 \t — s\}.

The distribution of the multi-marker process can be generated under the alternative similarly to the way it is generated under the null. The basic random process in both cases is the Ornstein-Uhlenbeck process. This processes describes the stochastic element in the behavior of markers for each chromosome. The deterministic element in the behavior is the mean function. This deterministic element is the difference between the null distribution and the distribution under the alternative. For the latter case a non-zero mean function is added for any chromosome carrying a QTL. Consequently, one can use the R function "OU.sim" that we wrote in the previous chapter in order to simulate the multi-marker process on a given chromosome. The vector of means may be added to the simulated process of marker-specific test statistics. We implement this approach in the function "add.qtl". This function takes as an input the matrix produced by "OU.sim", the location of markers, the coefficient of recombination 3, and two new parameters: "q", the location on the first chromosome of the QTL (measured in cM from the telomere); and "xi", the noncentrality parameter. The mean vector is added and the altered scanning process is returned:

> add.qtl <- function(Z,beta,markers,q,xi) + {

+ stop("Number of columns of simulated matrix

+ does not match the number of markers")

The function "stop" may be used in order to stop a function in the case of a fatal error. The argument of the function is printed out if the error occurs. Similarly, a warning may be produced, in the case of a nonfatal errors, with the function "warning".

The function "sweep" returns a matrix obtained from an input matrix by sweeping the elements of a vector. The first argument is the input matrix. The second argument is the margin over which the elements of the vector should be applied. The third argument is the vector, and the fourth argument is the binary function that produces the elements of the output matrix from the application of the binary function to the element of the input matrix and the appropriate element of the vector. More generally, this function may be applied to arrays, which are higher-dimension extension of matrices. Let us generate some paths of the resulting multi-marker process:

> for (i in 1:20) Z <- cbind(Z,OU.sim(beta,markers,n.iter=5)) Loading required package: MASS

> Z[,chr1] <- add.qtl(Z,beta,markers,q,xi) Error in add.qtl(Z, beta, markers, q, xi) :

Fig. 6.1. Sample paths of the absolute value of the Ornstein-Uhlenbeck process under the alternative.

Number of columns of simulated matrix does not match the number of markers

Note an erroneous application of the function "add.qtl" resulted in an error message, which helped us to detect the source of the error and debug the mistake. The function "OU.sim" did not produce an error message since we added to its definition the expression "require(MASS)".

The paths of the scanning process are plotted in Fig. 6.1. The paths for chromosome 1 are shown on the left plot and the paths of chromosomes 2-20 are shown on the right plot. Observe that the values of the test statistics in the middle of chromosome 1 have levels which are consistently large. However, occasional extreme values can occur also in other chromosomes. The code that produces the two plots is:

> plot(c(0,80),range(abs(Z)),type="n",xlab="chr 1",ylab="")

> for(i in 1:5) lines(markers,abs(Z[i,chr1]),col=gray(i/7))

> plot(c(0,80),range(abs(Z)),type="n",xlab="chrs 2-20",

+ for(i in 1:5) lines(markers,abs(Z[i,chr]),col=gray(i/7))

We detect a QTL if the maximum absolute value of the multi-marker process exceeds the significance threshold. Let us examine the distribution of this maximum, both when the signal is absent and when it is present. However, for the sake of determining the significance level we simulate these distributions for a whole-genome scan with markers at 0, 10, 20, .80 cM. The power is considered only in the context of the chromosome that contains a QTL. Specifically, the QTL is located 40 cM from the telomere on chromosome 1. Note that a marker happens to be located right at that spot. The noncentrality parameter at the QTL equals 4 in our simulations.

> for (i in 1:20) Z0 <- cbind(Z0,OU.sim(beta,markers))

> d0 <- density(apply(abs(Z0),1,max),from=1,to=7)

> d1 <- density(apply(abs(Z1),1,max),from=1,to=7)

> plot(d0,main="Densities of maximal statistics",

> legend(4.7,1,legend=c("Under H0","Under H1"),lty=1:2)

Examine the distributions of the test statistic under the two scenarios (Fig. 6.2). Note that although the distribution of the test statistic under the alternative tends to get higher values, still the two distributions cannot be separated perfectly. Any reasonable threshold that eliminates most of the ex-ceedences of the test statistic under the null distribution must eliminate some occurrences of the test statistic under the alternative distribution as well. The traditional way to resolve this difficulty is to set a threshold with an acceptable proportion of the null distribution above it. This proportion is the significance level of the test - the probability of falsely rejecting the null hypothesis - and is typically set to be equal to 5%. The proportion of the distribution under the alternative that is below the threshold corresponds to the probability of falsely accepting the null hypothesis, and is called the probability of a type II error. Under the alternative the statistical power corresponds in the plot to the portion above the threshold, i.e., one minus the probability of a type II error. The larger this probability the better.

In order to determine the power, one must specify the threshold. In the previous chapter the value of 3.56 was suggested as a threshold for a genome scan in the backcross design with inter-marker spacing of 10 cM. This threshold was derived from an analytical formula for the significance value. Let us see the actual significance level and power via simulations.

Densities of maximal statistics

Densities of maximal statistics

Fig. 6.2. The distribution of the test statistic under the null and under the alternative hypotheses.

Fig. 6.2. The distribution of the test statistic under the null and under the alternative hypotheses.

> mean(apply(abs(Z0),1,max) >= 3.56) [1] 0.0482

> mean(apply(abs(Z1),1,max) >= 3.56) [1] 0.7044

It follows that the power to detect a QTL, located 40 cM from the telomere and with an (asymptotic) noncentrality parameter of 4, is about 70%.

What would be the power if the QTL happened to be between two markers? near the end of the chromosome, rather than near the center? for a smaller value of the noncentrality parameter? a larger value? We can obtain answers to these questions using simulations. For example:

> mean(apply(abs(Z1),1,max) >= 3.56) [1] 0.6482

The task of exploring these question can be carried out much more efficiently once we have formulas for the statistical power, similarly to the formulas we have for the significance level. In the next section we describe such formulas.

### 6.2 An Analytic Approximation of the Power

As we saw, the power can be affected quite heavily by the location of the QTL on the chromosome. The probability of detecting a QTL may be substantially reduced if the QTL is located midway between markers, compared to the case where the QTL is located at the same location of a marker in the middle of the chromosome, and even more so if it is near an end of the chromosome. Since the formulas for approximating the power for the case when the QTL is between markers are substantially more complex than the formula when it is located exactly at a marker, we will present here only the formula for the latter case. For the former case, however, we do provide an R function, but not an explicit display of the formula. The interested reader may extract the mathematical expression from the code of the function.

When a marker is located at a QTL, detection occurs on that chromosome if either (i) the test statistic associated with the QTL/marker or (ii) the process associated with the flanking markers exceeds the threshold. The test statistic at the QTL has a normal distribution with mean £. Thus, the probability of the first case is simply the probability that such normal variable exceeds a threshold z. The mathematical derivation of the second case proceeds by conditioning on the value of the test statistic at the QTL, and analyzing the asymptotic conditional distribution of the process at the other markers. The resulting formula for a QTL not near either end of a chromosome is:

Pr( max\> z) « 1 - \$(z - |£|) + \$(z - |£|) [2 - v2/(z + |£|)] , (6.2)

i where v = v(z{2 (3A}l/2). The first term, 1 - \$(z - |£|), is the probability that the test statistic associated with the QTL/marker exceeds the threshold z. The second term corresponds to the probability of crossing the threshold by one or the other of the two flanking processes when the value of the statistic at the QTL is below the threshold.

Remark 6.1. Recall that for A « 0, i.e., when the distribution of markers on the chromosome is very dense, the correction term v is close to one.

Remark 6.2. When the QTL is located at the first or last marker on a chromosome, there are flanking markers only to one side. Then the approximation becomes

0 0