engineering

The engineering essay below has been submitted to us by a student in order to help you with your studies.

Numerical Comparisons Of Bioassay Methods Engineering Essay

The efficiency of pesticide to pests is widely studied. The level of a stimulus that results in a response by 50% of individuals in a population under study is an important characterizing parameter and it is denoted by median lethal dose (LD50) or median lethal concentration (LC50). Estimation of LC50 is a type of quantal response assays that belongs to qualitative indirect bioassays. In this report, seven methods of estimating LC50 are reviewed. The efficiency of them is discussed with reference to several normal distributions of tolerance in four different cases. Some modified methods are also discussed. Maximum likelihood method outperforms the other methods. The comparing result indicates that the modified Dragstedt-Behrens method and modified Reed-Muench method are good substitutes for the original ones in most of conditions.

Table of Contents

List of Figures 3

List of Tables 4

List of Tables 4

Acknowledgements 6

Acknowledgements 6

CHAPTER 1 - Introduction 1

1.1 Bioassay 1

1.2 LC50 2

1.3 Methods of estimating LC 50 2

1.4 Other Methods of estimating LC 50 3

CHAPTER 2 - Methods of Estimating LC50 4

2.1 The Dragstedt-Behrens Method 4

2.1.1 Algorithm 4

2.1.2 Confidence Interval 5

2.1.3 Modification and Improvement 5

2.1.4 Simulation Study 7

2.1.5 R-code 10

2.2 The Spearman-Kärber Method 15

2.2.1 Algorithm 16

2.2.2 Confidence Interval 16

2.2.3 Simulation study 17

2.2.4 R-code 19

2.3 The Reed-Muench Method 22

2.3.1 Algorithm 23

2.3.2 Modification and Improvement 24

2.3.3 Simulation study 25

2.3.4 R-code 28

2.4 The Thompson Moving Average Method 32

2.4.1 Algorithm 32

2.4.2 Modification and Improvement 33

2.4.3 Simulation study 35

2.4.4 R-code 36

2.5 The Shuster-Dietrich Method 41

The method was recommended by Shuster and Dietrich (1976) for estimating the dose response curves in quantal response assays. It is a general inverse regression procedure. 41

2.5.1 Algorithm 41

2.5.2 Simulation study 42

2.5.3 R-code 44

2.6 The Shuster-Yang Method 47

2.6.1 Algorithm 48

2.6.2 Simulation study 48

2.6.3 R-code 50

2.7 Maximum Likelihood Method 53

2.7.1 Algorithm 53

2.7.2 Simulation Study 54

2.7.3 R-code 55

CHAPTER 3 - Conclusion and Discussion 58

References Or Bibliography 60

References Or Bibliography 60

List of Figures

Figure 2.1 Plot of p* versus x, Dragstedt-Behrens method and Modified Dragstedt-Behrens method. 6

Figure 2.2 Plot of T versus x for Reed-Muench method and Modified Reed-Muench method. 24

Figure 2.3 Plot of P* versus x for moving average method 34

List of Tables

Table 2.1 Case 1, normal distribution. 7

Table 2.2 Case 1, Cauchy distribution. 7

Table 2.3 Case 2, normal distribution. 7

Table 2.4 Case 2, Cauchy distribution. 8

Table 2.5 Case 3, normal distribution. 8

Table 2.6 Case 3, Cauchy distribution. 8

Table 2.7 Case 4, normal distribution. 9

Table 2.8 Case 4, Cauchy distribution. 9

Table 2.9 Spearman-Kärber method, Case 1, normal distribution 17

Table 2.10 Spearman-Kärber method, Case 1, Cauchy distribution. 17

Table 2.11 Spearman-Kärber method, Case 2, normal distribution. 17

Table 2.12 Spearman-Kärber method, Case 2, Cauchy distribution. 18

Table 2.13 Spearman-Kärber method, Case 3, normal distribution. 18

Table 2.14 Spearman-Kärber method, Case 3, Cauchy distribution. 18

Table 2.15 Spearman-Kärber method, Case 4, normal distribution. 19

Table 2.16 Spearman-Kärber method, Case 4, Cauchy distribution. 19

Table 2.17 Reed- Muench method, Case 1, normal distribution. 25

Table 2.18 Reed- Muench method, Case 1, Cauchy distribution. 26

Table 2.19 Reed- Muench method, Case 2, normal distribution. 26

Table 2.20 Reed- Muench method, Case 2, Cauchy distribution. 26

Table 2.21 Reed- Muench method, Case 3, normal distribution. 26

Table 2.22 Reed- Muench method, Case 3, Cauchy distribution. 27

Table 2.23 Reed- Muench method, Case 4, normal distribution. 27

Table 2.24 Reed- Muench method, Case 4, Cauchy distribution. 27

Table 2.25 Case 1 35

Table 2.26 Case 2 35

Table 2.27 Case 3 36

Table 2.28 Case 4 36

Table 2.29 Case 1 42

Table 2.30 Case 2 43

Table 2.31 Case 3 43

Table 2.32 Case 4 44

Table 2.33 Case 1 48

Table 2.34 Case 2 49

Table 2.35 Case 3 49

Table 2.36 Case 4 49

Table 2.37 Case 1 54

Table 2.38 Case 2 54

Table 2.39 Case 3 54

Table 2.40 Case 4 55

Acknowledgements

Enter text for your Acknowledgements here. This section is optional. If you do not include an Acknowledgements, delete this section. Be sure to retain the “Page Break” that occurs after the “List of Tables” section above.

Introduction

Bioassay

Finney defined a biological assay as “an experiment for estimating the nature, constitution, or potency of a material (or of a process), by means of the reaction that follows its application to living matter.” (Finney, 1978) That is to say bioassay includes a lot of procedures. In those procedures, the amount or strength of an agent (usually is a drug) or stimulus is determined by a response (particular characteristic changes such as death) of a subject (usually is an animal or animal tissue). According to Webster's International Dictionary, bioassay is defined as “the estimation of the strength of a drug, etc., by comparing its effect on biological material, as animals or animal tissue, with those of a standard product.” For example, Xt units of drug perform as Xs units of the standard. When=2, it means one unit of the test drug is equivalent to 2 units of the standard with respect to biological activity.

Bioassay includes stochastic assay and non-stochastic assay. For stochastic assay, is affected by other factors besides the preparations. It’s very difficult to control those factors. Non-stochastic assay assume that is a constant that is independent of the subject. For majority of cases, the preparations are influenced by multiple factors such as species of the animal and environmental differences. So biometricians always talk about the stochastic assay.

Stochastic assay can be performed in two ways, directly or indirectly. In direct assays, the amount of stimulus is measured by holding the response fixed. In this situation, the dose is not the response but the variable of interest. However, the most common assay is indirect assay in which people observe a variable related to an event.

Indirect assay includes qualitative assay and quantitative assay. In qualitative assay, response always is only death. The number of response is the variable of interest. In quantitative indirect assays, both the response and the level of dosage are variables of interest.

LC50

Quantal assays belong to the qualitative indirect bioassays. In quantal response assays, a stimulus (e.g., dose of a drug) is applied to n experimental units, then there will be r units that are responding (e.g., death) and n-r that are not. Through these assays, we can obtain information about the amount of stimulus that produces a response. The quantal response assays are used to estimate the tolerance of individual, which refers as the individual effective dose (IED). The result will be the level of a stimulus corresponding to a particular response. A very important parameter, median lethal dose (LD50), denotes to the level of a stimulus that causes a response of 50% of individuals. There are some similar parameters such as median lethal concentration (LC50), median effective dose (ED50) and median effective concentration (EC50). We’ll concentrate on LC50 from the next chapter.

Why is LD50 so important? Why not use the minimal lethal dose under which all individuals response or the maximum lethal dose under which no one responses? Theoretically, these two definitions seem to be reasonable, but practically extreme members always exist in most samples no matter how many experimental units they have. Hence, the 100% or zero of lethal dose won’t be accurate. In some condition, LD20, LD90, etc. are also used, but LD50 are usually estimated more precisely.

Methods of estimating LC 50

From next chapter, the following seven methods of estimating LC50 will be studied:

(1). Dragster-Behrens method

(2). Spearman-Kärber method

(3). Reed-Muench method

(4). Thompson moving average method

(5). Shuster-Dietrich method

(6). Shuster-Yang method

(7). Maximum Likelihood method

Each method is simulated in 4 cases. Two cases have equally spaced log dose levels, which are 100.1, 100.2, 100.3, 100.4, 100.5, 100.6 and 100.7. The other two cases are unequally spaced dose levels, which are 1, 2, 3, 3.5, 5, 7 and 9.

When it comes to the number of subject at each dose level, two cases have same number of subject at each dose level (n= 10, 20, 30, 40, 50). In other two cases, numbers of subject at each dose level are different (n is generated randomly based on n=m±1, where m=10, 20, 30, 40, 50).

Other Methods of estimating LC 50

There are still some methods that we won’t discuss in this report such as the Dixon-Mood method (1948) and Litchfield-Wilcoxon method (1949). Dixon-Mood method is also known as the “up-and-down method” because of its unique procedure. In the procedure, if the subject survives, the dose for the next one will be increased; if it dies, the dose will be decreased. This method needs much less experiment subjects than others, so it’s very popular in many expensive and time-consuming experiments. Litchfield-Wilcoxon method uses a line drawn by eye to fit points for each dose and response data. This is a rapid graphic method for estimating LD50 and confidence limits. For the Litchfield-Wilcoxon method, accuracy is a conspicuous merit. Finney (1971) states, “The results are often very close to those from maximum likelihood estimation”.

In addition, many researchers such as Ramsey (1972), Chmiel (1976), Freeman (1980) and Davis (1972) studied other methods that require more statistical background.

Methods of Estimating LC50

2.1 The Dragstedt-Behrens Method

Dragstedt-Behrens method was independently introduced by Dragstedt and Lang (1928) and Behrens (1929). The algorithm of Dradstedt-Behrens method is very simple and is very easy to use. However its numerical simplicity is considered to be the “only merit” and it has many limitations. For example, it “may behave reasonably well if the data are equally spaced in X have equal n at each X, and are moderately symmetric” (Finney 1971). Besides, Dragstedt-Behrens method has no sound theoretical basis.

2.1.1 Algorithm

Suppose an experiment looking for LC50 produces the following data set:

Dose Level(X)

Size(n)

Response(r)

x1

n1

r1

x2

n2

r2

xk

nk

rk

where k is the number of dose levels used in the experiment, and X is the log-transformed dose level. Usually, xi’s are arranged in the ascending order.

Dradstedt-Behrens method uses the following procedure to find out the estimate of LC50.

(1). At each dose level, calculate

T1(x)= Total of all values of r’s for doses equal to or less than x,

T2(x)= Total of all values of r’s for doses equal to or greater than x,

and .

(2). If there is a dose level, say xi, such that p*(x)=0.5, then let m=xi,

(3). If there is no value of p* being equal to 0.5, then find the biggest dose level, say xi, which p* is less than 0.5. Then let

.

(4). The LC50 is estimated by LC50=10m

2.1.2 Confidence Interval

If the X values are equally spaced and the numbers of subject in each dose level are the same, then a crude estimate of the standard error of m can be estimated by

,

where h is the difference between two adjacent dose levels, n is the number of subject in each dose level and IR is an estimate of the inter-quartile range of log(LC75) – log(LC25), which can be obtained by using similar algorithm as finding the LC50. Then a confidence interval of LC50 with confidence level 1- α can be constructed as

,

where zα/2 is the upper 100 α/2 percentile of the standard normal distribution.

2.1.3 Modification and Improvement

(1). Non-equally spaced X or different number of subjects at each dose level.

If the values of X are roughly equally spaced or the numbers of subjects at the chosen dose levels are roughly same, the above procedure still works. The algorithm of find LC50 is stay unchanged. To construct the confidence interval using the above mentioned SE(m), one can use the following formulas to calculate h and n.

, and .

(2). Qubic-Fitting method of finding LC50.

Dragstedt-Behrens method uses linear interpolation to find the estimation of LC50. It only uses the information from two points (xi, p*(xi)), and (xi, p*(xi+1)), where xi is the biggest dose level when P*(xi) ≤0.5. The following is a typical plot of p* versus x.

The general pattern of the plot shows a sigmoid trend. Therefore, a cubic polynomial, estimated from all data, can be used to find a better estimation of LC50. One can easily show that p*(x) is an increasing function of the dose level x.

Figure 2.1 Plot of p* versus x, Dragstedt-Behrens method and Modified Dragstedt-Behrens method.

2.1.4 Simulation Study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.1 Case 1, normal distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

2.8179

0.0379

0.968

1.3836

2.7727

0.0322

n=20

2.7920

0.0161

0.972

1.2512

2.7662

0.0158

n=30

2.7890

0.0104

0.986

1.1988

2.7618

0.0106

n=40

2.7892

0.0067

0.984

1.1680

2.7617

0.0073

n=50

2.7877

0.0059

0.986

1.1492

2.7597

0.0065

Table 2.2 Case 1, Cauchy distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

2.7766

0.0563

0.966

1.4581

2.7263

0.0539

n=20

2.7506

0.0275

0.960

1.3033

2.7231

0.0282

n=30

2.7470

0.0195

0.948

1.2377

2.7160

0.0208

n=40

2.7458

0.0143

0.954

1.2033

2.7155

0.0161

n=50

2.7388

0.0125

0.956

1.1804

2.7137

0.0147

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.3 Case 2, normal distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

2.7914

0.0392

0.968

1.3840

2.7522

0.0341

n=20

2.7996

0.0145

0.978

1.2458

2.7719

0.0142

n=30

2.7902

0.0095

0.986

1.1973

2.7613

0.0098

n=40

2.7906

0.0063

0.990

1.1691

2.7617

0.0070

n=50

2.7861

0.0063

0.984

1.1503

2.7581

0.0069

Table 2.4 Case 2, Cauchy distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

2.7014

0.0668

0.948

1.4651

2.6528

0.0647

n=20

2.7659

0.0230

0.974

1.2960

2.7413

0.0256

n=30

2.7491

0.0171

0.966

1.2353

2.7222

0.0190

n=40

2.7423

0.0127

0.964

1.2038

2.7118

0.0157

n=50

2.7292

0.0137

0.946

1.1826

2.7056

0.0161

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.5 Case 3, normal distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

4.1028

0.4093

0.914

1.7975

4.2918

0.2894

n=20

4.1246

0.2331

0.918

1.5073

4.2754

0.1478

n=30

4.1150

0.1805

0.884

1.3967

4.2596

0.1012

n=40

4.1177

0.1391

0.890

1.3362

4.2566

0.0744

n=50

4.1197

0.1257

0.864

1.2937

4.2546

0.0651

Table 2.6 Case 3, Cauchy distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

3.9584

0.4703

0.930

1.8622

4.1578

0.3164

n=20

3.9696

0.3201

0.876

1.5494

4.1449

0.1857

n=30

3.9642

0.2728

0.808

1.4293

4.1311

0.1437

n=40

3.9651

0.2351

0.786

1.3626

4.1258

0.1186

n=50

3.9550

0.2363

0.730

1.3166

4.1212

0.1110

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.7 Case 4, normal distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

3.8968

0.4800

0.870

1.7948

4.0754

0.3173

n=20

4.1847

0.1934

0.944

1.5037

4.3324

0.1371

n=30

4.1207

0.1627

0.898

1.3934

4.2596

0.0949

n=40

4.1114

0.1351

0.902

1.3351

4.2406

0.0737

n=50

4.0807

0.1462

0.826

1.2989

4.2185

0.0727

Table 2.8 Case 4, Cauchy distribution.

Dragstedt-Behrens Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

Mean

MSE

n=10

3.7138

0.6480

0.884

1.8735

3.8938

0.4568

n=20

4.0620

0.2609

0.914

1.5452

4.2217

0.1608

n=30

3.9849

0.2528

0.832

1.4218

4.1430

0.1314

n=40

3.9607

0.2306

0.790

1.3618

4.1102

0.1220

n=50

3.9181

0.2627

0.660

1.3209

4.0866

0.1276

Conclusion: From above simulation study by normal distribution, one can see that in the case of equally spaced dose levels, Dragstedt-Beherens method and the cubic polynomial fitting method provide similar results; in the case of unequally spaced dose levels, the cubic polynomial fitting method outperforms Dragstedt-Behrens method.

If the case of equally spaced dose levels, Dragstedt-Beherens method provides confidence intervals with empirical levels which is slightly larger than the nominal level 95%, but it performs poorly in the case of unequally spaced dose levels, the empirical confidence levels are all significantly smaller than the nominal level 95%, and more seriously, the empirical confidence levels are decreasing when the numbers of subject at each dose level increases.

When it comes to Cauchy distribution, the result is not as good as the normal distribution.

2.1.5 R-code

The simulation study is conducted by using the following R program.

###########################

# Dragstedt-Behrens Method #

###########################

rm(list=ls())

set.seed(987654)

Eqx.or.Not=1 # 0=No, 1=Yes

Eqn.or.Not=0 # 0=No, 1=Yes

# Define Dose Levels

# Equally Spaced Transformed Dose Levels

Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

(1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

Num.Dose.Level=length(Dose.Level);

# Dose Level Transformations

Tran.Dose.Level=log10(Dose.Level);

# Define Tolerance Distribution

Log.True.Mean=log10(mean(Dose.Level));

Log.True.Sd=log10(sd(Dose.Level));

Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)}; # normal distribution

#Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)};

# Cauchy distribution

# Generate True Response Rate

True.Res.Rate=Tol.Dist(Tran.Dose.Level);

# Generate Sample

Mean.LC50=Mean.LC50.new=Ep.Level=Mean.Ep.Length=MSE.LC50=MSE.LC50.new =rep(0,5);

j=1;

for(size in c(10,20,30,40,50))

{

n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1- Eqn.or.Not)*(size+round(3*runif(Num.Dose.Level,-1,1)))

Rep.in.Sim=500;

m50=m50.new=m25=m75=SEs=rep(0,Rep.in.Sim);

# Initial values of Log(LC50,LC25,LC75)

for(run in seq(Rep.in.Sim))

{

Sim.Res=rep(0,Num.Dose.Level)

for(i in seq(Num.Dose.Level))

{

Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

}

# Dragstedt-Behrens Method

T1=cumsum(Sim.Res);

T2=rev(cumsum(rev(n-Sim.Res)));

pstar=T1/(T1+T2);

plow50=max(pstar[pstar<=0.5])

pupper50=min(pstar[pstar>=0.5])

xlow50=Tran.Dose.Level[pstar==plow50]

xupper50=Tran.Dose.Level[pstar==pupper50]

m50[run]=ifelse(

pupper50==plow50, plow50,

{

w=(pupper50-0.5)/(pupper50-plow50);

w*xlow50+(1-w)*xupper50

}

)

z1=Tran.Dose.Level;

z2=z1*z1;

z3=z2*z1;

my3reg=lm(pstar~z1+z2+z3);

a0=my3reg$coefficient[1]-0.5;

a1=my3reg$coefficient[2];

a2=my3reg$coefficient[3];

a3=my3reg$coefficient[4];

myfun=function(x){a0+a1*x+a2*x^2+a3*x^3};

root=uniroot(myfun, lower=min(z1),upper=max(z1),tol = 0.0001)

m50.new[run]=root$root;

plow25=max(pstar[pstar<=0.25])

pupper25=min(pstar[pstar>=0.25])

xlow25=Tran.Dose.Level[pstar==plow25]

xupper25=Tran.Dose.Level[pstar==pupper25]

m25[run]=ifelse(

pupper25==plow25, plow25,

{

w=(pupper25-0.25)/(pupper25-plow25);

w*xlow25+(1-w)*xupper25

}

)

plow75=max(pstar[pstar<=0.75])

pupper75=min(pstar[pstar>=0.75])

xlow75=Tran.Dose.Level[pstar==plow75]

xupper75=Tran.Dose.Level[pstar==pupper75]

m75[run]=ifelse(

pupper75==plow75, plow75,

{

w=(pupper75-0.75)/(pupper75-plow75);

w*xlow75+(1-w)*xupper75

}

)

# h=(max(Tran.Dose.Level)-min(Tran.Dose.Level))/(Num.Dose.Level-1);

# SE[run]=(h/size)*sqrt(sum(Sim.Res*(size-Sim.Res))/(size))

}

h=(max(Tran.Dose.Level)-min(Tran.Dose.Level))/(Num.Dose.Level-1);

IR=m75-m25;

SE=sqrt(0.79*h*IR/round(mean(n)))

L.end=10^(m50-1.96*SE)

R.end=10^(m50+1.96*SE)

freq=sum((10^Log.True.Mean>=L.end)*(10^Log.True.Mean<=R.end))

Mean.Ep.Length[j]=mean(10^(2*1.96*SE))

Ep.Level[j]=freq/Rep.in.Sim;

Mean.LC50[j]=mean(10^m50)

Mean.LC50.new[j]=mean(10^m50.new)

MSE.LC50[j]=mean((10^m50-10^Log.True.Mean)^2);

MSE.LC50.new[j]=mean((10^m50.new-10^Log.True.Mean)^2);

j=j+1;

}

Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50,Ep.Level,Mean.Ep.Length,Mean.LC50.new,MSE.LC50.new))

dimnames(Result)=list(c("","","","",""), c("Sample.Size","Mean","MSE","Empirical.CL","Mean.Ep.Length","Mean.Cubic","MSE.Cubic"))

{

cat(" Simulation Result","\n");

cat(" Equally Spaced X:", Eqx.or.Not,","," ", "Same Number of Subjects:", Eqn.or.Not,"\n\n")

Result

}

pp=as.matrix(cbind(z1,pstar))

plot(z1,pstar,xlab="Transformed Dose Level",ylab="p*")

xx=c(pp[4,1],pp[5,1])

yy=c(pp[4,2],pp[5,2])

lines(xx,yy,col="blue",lwd=2)

xx=seq(0,1,by=0.01);

lines(xx,a0+0.5+a1*xx+a2*xx^2+a3*xx^3,col="red",lwd="2");

2.2 The Spearman-Kärber Method

The Spearman-Kärber method was first proposed by Spearman (1908) and independently reintroduced by Kärber (1931). This method is very precise and particularly easy to use, so it is still very popular in many fields, such as analysis of psychometric functions. In some experimental conditions, it can provide more power to detect differences in parameters across. There are some methods that are quite easy to calculate, such as Reed-Muench method and Dragstedt-Behrens method, but Spearman-Kärber is easier and often markedly superior.

The Spearman-Kärber method obtains the estimate of m (m=log(LC50)). Although the method does not need equal spacing of doses level and equal replications, it has two requirements:1) the doses extend over the range from 0% to 100% response; 2) the response increasing rate between successive doses should concentrate on the center of the interval.

2.2.1 Algorithm

(1). Calculate

Xi= log (ri/ni), pi=ri/n.

Let p1=0% and pk=100%.

(2). The estimation of m= log10(LD50) is

.

(3). In the case of equally spaced doses with Xi+1-Xi=d, the formula can be reduced to the simple expression:

.

(4). If the number of replications is also constant, then the m is

.

2.2.2 Confidence Interval

The variance of m is obtained by replacing the variance of pi with pi(1-pi):

,

, the 95% confidence interval of m is m±zα/2SE(m), where zα/2 is the upper 100 zα/2 percentile of the standard normal distribution. Then, the 95% confidence interval of LC50 is 10m± zα/2SE(m).

2.2.3 Simulation study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.9 Spearman-Kärber method, Case 1, normal distribution

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

2.8114

0.0385

0.920

1.276

n=20

2.8054

0.0189

0.926

1.1937

n=30

2.8046

0.0126

0.944

1.1576

n=40

2.8043

0.0087

0.946

1.1357

n=50

2.8004

0.0074

0.944

1.1203

Table 2.10 Spearman-Kärber method, Case 1, Cauchy distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

2.8331

0.076

0.428

1.1127

n=20

2.8273

0.0355

0.368

1.0566

n=30

2.8242

0.0245

0.26

1.0379

n=40

2.8227

0.018

0.252

1.0284

n=50

2.8175

0.0148

0.172

1.0226

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.11 Spearman-Kärber method, Case 2, normal distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

2.8140

0.0413

0.914

1.2818

n=20

2.8044

0.0179

0.924

1.1906

n=30

2.8046

0.0120

0.934

1.1544

n=40

2.8059

0.0088

0.952

1.1357

n=50

2.8004

0.0076

0.946

1.1215

Table 2.12 Spearman-Kärber method, Case 2, Cauchy distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

2.8329

0.0774

0.452

1.1121

n=20

2.8268

0.0355

0.346

1.0567

n=30

2.8268

0.0243

0.264

1.0379

n=40

2.8234

0.0185

0.246

1.0284

n=50

2.8176

0.0151

0.188

1.0226

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.13 Spearman-Kärber method, Case 3, normal distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

3.9410

0.4000

0.848

1.5753

n=20

3.9368

0.2952

0.758

1.3927

n=30

3.9280

0.2629

0.666

1.3132

n=40

3.9192

0.2493

0.570

1.2684

n=50

3.9290

0.2332

0.522

1.2364

Table 2.14 Spearman-Kärber method, Case 3, Cauchy distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

3.7014

0.6491

0.188

1.1667

n=20

3.6955

0.5572

0.080

1.0826

n=30

3.6865

0.5279

0.016

1.0549

n=40

3.6752

0.5235

0.010

1.0412

n=50

3.6835

0.5042

0.004

1.0328

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.15 Spearman-Kärber method, Case 4, normal distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

3.9417

0.3743

0.842

1.5533

n=20

3.9384

0.2917

0.754

1.3886

n=30

3.9302

0.2584

0.678

1.3085

n=40

3.9203

0.2484

0.586

1.2692

n=50

3.9250

0.2354

0.516

1.2363

Table 2.16 Spearman-Kärber method, Case 4, Cauchy distribution.

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Empirical Length

n=10

3.7073

0.6345

0.196

1.1677

n=20

3.6984

0.5535

0.080

1.0827

n=30

3.6857

0.5261

0.014

1.0549

n=40

3.6748

0.5241

0.012

1.0412

n=50

3.6840

0.5023

0.006

1.0328

Conclusion: The Spearman-Kärber method provides smaller MSE when the dose levels are equal. In addition, in the case of equally spaced dose levels, the confidence intervals with empirical levels are almost the nominal level 95%, but they are significantly smaller than the nominal 95% when the space of log dose levels are not equal. One can get the conclusion that the Spearman-Kärber method is more suitable for the equally spaced dose levels cases.

2.2.4 R-code

The simulation study is conducted by using the following R program.

#########################

# Spearman-Karber Method #

#########################

rm(list=ls())

set.seed(987654)

Eqx.or.Not=0 # 0=No, 1=Yes

Eqn.or.Not=0 # 0=No, 1=Yes

# Define Dose Levels

Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

(1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

Num.Dose.Level=length(Dose.Level);

# Dose Level Transformations

Tran.Dose.Level=log10(Dose.Level);

# Define Tolerance Distribution

Log.True.Mean=log10(mean(Dose.Level));

Log.True.Sd=log10(sd(Dose.Level));

Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)};# normal distribution

#Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)};

# Cauchy distribution

# Generate True Response Rate

True.Res.Rate=Tol.Dist(Tran.Dose.Level);

# Generate Sample

m=Mean.LC50=Ep.Level=Mean.Ep.length=MSE.LC50=Var=rep(0,5);

j=1;

for(size in c(10,20,30,40,50))

{

n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1-Eqn.or.Not) * ( size+round( 3*runif(Num.Dose.Level,-1,1)))

#n=rep(size,Num.Dose.Level); # Population Size at Each Dose Level

Rep.in.Sim=500;

m=rep(0,Rep.in.Sim); # Initial values of Log(LC50)

for(run in seq(Rep.in.Sim))

{

Sim.Res=rep(0,Num.Dose.Level)

for(i in seq(Num.Dose.Level))

{

Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

}

# Spearman-Karber Method

p=Sim.Res/n;

P=c(0,p,1);

K=length(P);

X=c(0,Tran.Dose.Level,1);

N=rep(size,K)

m[run]=sum((P[2:K]-P[1:(K-1)])*(X[1:(K-1)]+X[2:K]))/2

Var[run]=sum(P[2:(K-1)]*(1-P[2:(K-1)])*(X[3:(K)]-X[1:(K-2)])^2/(4*N[2:(K-1)]))

# Var of each m

}

SE=sqrt(Var/N)

L.end=10^(m-1.96*SE)

R.end=10^(m+1.96*SE)

freq=sum((10^Log.True.Mean>=L.end)*(10^Log.True.Mean<=R.end))

Mean.Ep.length[j]=mean(10^(2*1.96*SE))

Ep.Level[j]=freq/Rep.in.Sim;

Mean.LC50[j]=mean(10^m)

MSE.LC50[j]=mean((10^m-10^Log.True.Mean)^2);

j=j+1;

}

Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50,Ep.Level, Mean.Ep.length))

dimnames(Result)=list(c("","","","",""), c("Sample.Size","Mean","MSE","Empirical.CL", "Mean.Ep.length"))

{

cat(" Simulation Result","\n");

cat(" Equally Spaced X:", Eqx.or.Not,","," ", "Same Number of Subjects:", Eqn.or.Not,"\n\n")

Result

}

2.3 The Reed-Muench Method

Reed and Muench (1938) proposed “a simple method of estimating fifty per cent endpoints”. It does not need large numbers of animal tests at dilutions near the value for LD50, but use a wide range of possible variations. This method was very popular because of its numerical simplicity and no sound theoretical basis. Pittman and Lieberman (1948) showed evidence of the inferiority of the method relative to maximum likelihood. However, this method (as well as the extreme effective dose and Dragstedt-Behrens methods) is not recommended. They “ought never to be used” because their tests were not valid and less efficient then other simple methods. They should be forgotten, “except as part of statistical history” (Finney 1978). Reed-Muench method assumes that any subject responding to a given dose of an agent would respond to all higher doses; and that any subject not responding to a given dose would not respond to a lower dose.

2.3.1 Algorithm

(1). Calculate

T1(x)= Total of all values of r for doses equal to or less than x,

T2(x)= Total of all values of r for doses equal to or greater than x,

(2). If one of the doses used in the experiment has T1(xi)= T2(xi) , then let m=xi.

(3). If there is no xi satisfying T1(xi)= T2(xi), then find xi and xi+1 such that T1(xi)<T2(xi), T1(xi+1)>T2(xi+1).

(4). Using bi-linear interpolation to find out the x-coordinate of the intersection. When T1(x)=T2(x), m is the solution of

m will be estimated by this x- value.

(5). The LC50 is estimated by

.

2.3.2 Modification and Improvement

Reed-Muench method only used the information around the intersection. We expect a better estimator that will obtain if it will involve all the information in the data set. For example, we can fit quadratic curves for {xi, T1(xi)}, i=1, 2, …, k, and {xi, T2(xi)}, i=1, 2, …, k, respectively. The x-coordinated of the intersection of these two quadratic curves can be used to estimate LC50. If the estimated curves from(x, T1) is a1x2+b1x+c1, and the estimated curves from (x, T2) is a2x2+b2x+c2, then

If xi<m’<xi+1

otherwise

Figure 2.2 Plot of T versus x for Reed-Muench method and Modified Reed-Muench method.

2.3.3 Simulation study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.17 Reed- Muench method, Case 1, normal distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

2.7901

0.0368

2.8069

0.0377

n=20

2.7898

0.0176

2.8053

0.0183

n=30

2.7860

0.0118

2.8004

0.0126

n=40

2.7875

0.0078

2.8021

0.0084

n=50

2.7858

0.0070

2.8003

0.0073

Table 2.18 Reed- Muench method, Case 1, Cauchy distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

2.7400

0.0531

2.7489

0.0573

n=20

2.7408

0.0273

2.7472

0.0291

n=30

2.7355

0.0192

2.7408

0.0203

n=40

2.7366

0.0140

2.7435

0.0144

n=50

2.7346

0.0128

2.7411

0.0129

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.19 Reed- Muench method, Case 2, normal distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

2.7721

0.0379

2.7828

0.0358

n=20

2.7944

0.0159

2.8106

0.0183

n=30

2.7836

0.0108

2.8100

0.0127

n=40

2.7888

0.0074

2.8054

0.0089

n=50

2.7831

0.0074

2.8000

0.0076

Table 2.20 Reed- Muench method, Case 2, Cauchy distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

2.6712

0.0628

2.6714

0.0626

n=20

2.7618

0.0239

2.7686

0.0276

n=30

2.7416

0.0175

2.7567

0.019

n=40

2.7369

0.0128

2.7415

0.0146

n=50

2.7258

0.0142

2.7337

0.0138

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.21 Reed- Muench method, Case 3, normal distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

4.1725

0.3609

4.4655

0.3094

n=20

4.1472

0.2160

4.4592

0.1592

n=30

4.1286

0.1701

4.4514

0.1088

n=40

4.1275

0.1357

4.4510

0.0785

n=50

4.1251

0.1260

4.4492

0.0683

Table 2.22 Reed- Muench method, Case 3, Cauchy distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

4.0279

0.4144

4.321

0.318

n=20

3.9946

0.3018

4.3156

0.1629

n=30

3.9774

0.2616

4.3084

0.111

n=40

3.9702

0.2334

4.3062

0.0791

n=50

3.9628

0.2292

4.3006

0.0689

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.23 Reed- Muench method, Case 4, normal distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

3.9529

0.4272

4.2815

0.2797

n=20

4.2132

0.1811

4.5148

0.1678

n=30

4.1335

0.1582

4.4771

0.1107

n=40

4.1187

0.1323

4.4530

0.0787

n=50

4.0849

0.1467

4.4220

0.0638

Table 2.24 Reed- Muench method, Case 4, Cauchy distribution.

Reed- Muench method

Modified method

mean

MSE

mean

MSE

n=10

3.7643

0.5892

4.0705

0.385

n=20

4.086

0.2411

4.3926

0.1613

n=30

3.9969

0.2392

4.3445

0.1049

n=40

3.966

0.2289

4.3041

0.0788

n=50

3.9217

0.2626

4.2724

0.0724

Conclusion: From the above simulation study, one can see that in the case of equally spaced dose levels, Reed-Muench method and the modified method provide similar results; In the case of unequally spaced dose levels, the modified method outperforms Reed-Muench method.

The comparing result sufficiently indicates that the modified method is a good substitute for the original one in most of conditions.

2.3.4 R-code

The simulation study is conducted by using the following R program.

######################

# Reed-Muench Method #

######################

rm(list=ls())

set.seed(987654)

Eqx.or.Not=1 # 0=No, 1=Yes

Eqn.or.Not=1 # 0=No, 1=Yes

# Define Dose Levels

Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

(1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

Num.Dose.Level=length(Dose.Level);

# Dose Level Transformations

Tran.Dose.Level=log10(Dose.Level);

# Define Tolerance Distribution

Log.True.Mean=log10(mean(Dose.Level));

Log.True.Sd=log10(sd(Dose.Level));

Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)}; # normal distribution

#Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)};

# Cauchy distribution

# Generate True Response Rate

True.Res.Rate=Tol.Dist(Tran.Dose.Level);

# Generate Sample

Mean.LC50=Mean.LC50.new=Ep.Level=MSE.LC50=MSE.LC50.new=rep(0,5);

j=1;

for(size in c(10,20,30,40,50))

{

n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1-Eqn.or.Not) * ( size+round( 3*runif(Num.Dose.Level,-1,1)))

#n=rep(size,Num.Dose.Level); # Population Size at Each Dose Level

Rep.in.Sim=500;

m=m.new=rep(0,Rep.in.Sim); # Initial values of Log(LC50,LC25,LC75)

for(run in seq(Rep.in.Sim))

{

Sim.Res=rep(0,Num.Dose.Level)

for(i in seq(Num.Dose.Level))

{

Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

}

# Dragstedt-Behrens Method

X=Tran.Dose.Level;

T1=cumsum(Sim.Res);

T2=rev(cumsum(rev(n-Sim.Res)));

T1.low=max(T1[T1<=T2])

T2.low=min(T2[T1<=T2])

Xlow=Tran.Dose.Level[T1==T1.low]

T1.upper=min(T1[T1>T2])

T2.upper=max(T2[T1>T2])

Xupper=Tran.Dose.Level[T1==T1.upper]

m[run]=(Xlow+(Xupper-Xlow)*(T2.low-T1.low)/(T1.upper-T2.upper-T1.low+T2.low))

# Modified Reed-muench Method

X=Tran.Dose.Level;

T1=cumsum(Sim.Res);

T2=rev(cumsum(rev(n-Sim.Res)));

X2=X^2;

myreg1=lm(T1~X+X2);

myreg2=lm(T2~X+X2);

summary(myreg1)

c1=myreg1$coeff[1] #coefficients for reg1

b1=myreg1$coeff[2]

a1=myreg1$coeff[3]

summary(myreg2)

c2=myreg2$coeff[1] #coefficients for reg2

b2=myreg2$coeff[2]

a2=myreg2$coeff[3]

delta=sqrt((b1-b2)^2-4*(a1-a2)*(c1-c2))

root1=(b2-b1-delta)/(2*(a1-a2))

root2=(b2-b1+delta)/(2*(a1-a2))

if(root1>=Xlow & root1<=Xupper) #get the real root by regression methord

{

m.new[run]=root1

}

else

{

m.new[run]=root2

}

}

Mean.LC50[j]=mean(10^m)

Mean.LC50.new[j]=mean(10^m.new)

MSE.LC50[j]=mean((10^m-10^Log.True.Mean)^2);

MSE.LC50.new[j]=mean((10^m.new-10^Log.True.Mean)^2);

j=j+1;

}

plot(X, T1, xlim=range(X), ylim=range(T1, T2), xlab=' Trasformed Dose Level', ylab='T')

points(X, T2, xlab='Trasformed Dose Level', ylab='T')

lines(X,myreg1$fitted,col=c("red"))

lines(X,myreg2$fitted,col=c("red"))

lines(c(Xlow, Xupper),c(T1.low,T1.upper),col=c("blue"));

lines(c(Xlow, Xupper),c(T2.low,T2.upper),col=c("blue"))

Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50,Mean.LC50.new,MSE.LC50.new))

dimnames(Result)=list(c("","","",""),("Sample.Size","Mean","MSE","Mean.LC50.new","MSE.LC50.new"))

Result

2.4 The Thompson Moving Average Method

Moving average method was firstly discussed by Sheppard (1914). Thompson (1947) proposed a method utilizing moving average to estimate the LC50 of bioassay data. It has some relationship with the Spearman–Kärber method. Since the Spearman–Kärber method needs the dose interval from almost 0 to 100%, it’s wasteful if there is an existing approximation to the LC50. The Thompson method is sufficiently different because it would average all sets of successive values of p and plot average versus middle dose. Finally an estimate of the LC50 can be obtained by linear interpolation. Unlike other methods of estimating LC50, the Thompson method cannot estimate any other percentage points. The Thompson moving average method was considered a ‘basic’ one in finding LC50. For example, log10(LC50) can be obtained by plotting (pi+pi+1)/2 versus (xi+xi+1)/2 and simple linear interpolation. A three-term moving average was recommended by Topley and Wilson.

2.4.1 Algorithm

(1). Calculate

,

.

(2). If there are P*i and P*i+1following P*i <0.5< P*i+1, then the estimated value of log10(LC50) is

.

(3). If there is a successive proportions P*i =0.5, then m=xi.

(4). The LC50 is estimated by

.

2.4.2 Modification and Improvement

(1). P* is not monotonic increasing near the 0.5.

If the value of P* is not monotonic increasing near the 0.5, there will be more than one m, which is unreasonable. We drop a P* to make the rest of them to be monotonic increasing, then the Thompson moving average method can be used.

(2). Qubic-Fitting method of finding LC50.

Thompson moving average method using linear interpolation to find the estimation of LC50, which only uses the information from two points( xi, P*i)and ( xi+1, P*i+1), where P*i <0.5< P*i+1. The following is a typical plot of P*versus x.

Figure 2.3 Plot of P* versus x for moving average method

The general pattern of the plot shows a sigmoid trend. Therefore, a cubic polynomial, estimated from all data, can be used to find a better estimation of LC50.

One can easily show that P*(x) is an increasing function of the dose level .

2.4.3 Simulation study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.25 Case 1

Thompson Moving Average Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Mean

MSE

n=10

2.7936

0.0460

2.7957

0.0476

n=20

2.7934

0.0203

2.7939

0.0216

n=30

2.7862

0.0141

2.7850

0.0150

n=40

2.7894

0.0091

2.7883

0.0098

n=50

2.7883

0.0085

2.7873

0.0091

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.26 Case 2

Thompson Moving Average Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Mean

MSE

n=10

2.7952

0.0466

2.7979

0.0480

n=20

2.7920

0.0190

2.7913

0.0206

n=30

2.7853

0.0133

2.7833

0.0143

n=40

2.7906

0.0090

2.7890

0.0098

n=50

2.7876

0.0087

2.7867

0.0093

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.27 Case 3

Thompson Moving Average Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Mean

MSE

n=10

4.4014

0.7447

4.4273

0.7973

n=20

4.5839

0.7979

4.6754

0.9481

n=30

4.3549

0.3565

4.3667

0.3167

n=40

4.5744

0.2913

4.5073

0.2797

n=50

4.1539

0.2196

4.1215

0.1977

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.28 Case 4

Thompson Moving Average Method

Cubic Polynomial Fitting

Sample Size

Mean

MSE

Mean

MSE

n=10

4.4312

0.8003

4.4223

0.8800

n=20

4.6886

0.8087

4.7146

1.0268

n=30

4.3464

0.3731

4.3713

0.3420

n=40

4.5828

0.3029

4.5164

0.2944

n=50

4.1673

0.2197

4.1356

0.1985

Conclusion: The modified method doesn't have significant improvement comparing to the original one. Both methods perform better when the space of log dose level is equal.

2.4.4 R-code

The simulation study is conducted by using the following R program.

#################################

# Thompson Moving Average Method #

#################################

rm(list=ls())

set.seed(987654)

Eqx.or.Not=1 # 0=No, 1=Yes

Eqn.or.Not=0 # 0=No, 1=Yes

# Define Dose Levels

# Equally Spaced Transformed Dose Levels

Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

(1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

Num.Dose.Level=length(Dose.Level);

# Dose Level Transformations

Tran.Dose.Level=log10(Dose.Level);

# Define Tolerance Distribution

Log.True.Mean=log10(mean(Dose.Level));

Log.True.Sd=log10(sd(Dose.Level));

Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)}; # normal distribution

#Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)};

# Cauchy distribution

# Generate True Response Rate

True.Res.Rate=Tol.Dist(Tran.Dose.Level);

# Generate Sample

Mean.LC50=Mean.LC50.new=Ep.Level=MSE.LC50=MSE.LC50.new=rep(0,5);

j=1;

for(size in c(10,20,30,40,50))

{

n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1- Eqn.or.Not)*(size+round(3*runif(Num.Dose.Level,-1,1)))

Rep.in.Sim=25;

m=m.new=SE=rep(0,Rep.in.Sim); # Initial values of Log(LC50)

for(run in seq(Rep.in.Sim))

{

Sim.Res=rep(0,Num.Dose.Level)

for(i in seq(Num.Dose.Level))

{

Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

}

# Thompson Moving Average Method

X=Tran.Dose.Level;

p=Sim.Res/n

n.mean=mean(n)

k=Num.Dose.Level

P=rep(0,k-2)

for(i in seq(k-2))

{

P[i+1]=(p[i]+n[i+1]*p[i+1]/n.mean+p[i+2])/(2+n[i+1]/n.mean)

}

D=rep(0,k-2)

for(i in seq(k-2))

{

D[i+1]=(X[i]+n[i+1]*X[i+1]/n.mean+X[i+2])/(2+n[i+1]/n.mean)

}

Plow=max(P[P<=0.5])

Pupper=min(P[P>=0.5])

I=length(Plow)

Dlow=D[P==Plow]

Dupper=D[P==Pupper]

Xlow=X[P==Plow]

Xupper=X[P==Pupper]

plow=p[P==Plow]

pupper=p[P==Pupper]

if (Plow==0.5)

{

m[run]=Xlow

}

else if (I==1) #when p's are increasing

{

m[run]=Dlow+(Dupper-Dlow)*(0.5-Plow)/(Pupper-Plow)

}

else #when p's are not increasing

{

p=p(-c(pupper))

X=X(-c(Xupper))

}

# Quadratic-Fitting Method

z1=Tran.Dose.Level[2:(k-1)];

z2=z1*z1;

z3=z2*z1;

Pstar=P[2:(k-1)]

my3reg=lm(Pstar~z1+z2+z3);

a0=my3reg$coefficient[1]-0.5;

a1=my3reg$coefficient[2];

a2=my3reg$coefficient[3];

a3=my3reg$coefficient[4];

myfun=function(x){a0+a1*x+a2*x^2+a3*x^3};

root=uniroot(myfun, lower=min(z1),upper=max(z1),tol = 0.0001)

m.new[run]=root$root;

}

Mean.LC50[j]=mean(10^m)

Mean.LC50.new[j]=mean(10^m.new)

MSE.LC50[j]=mean((10^m-10^Log.True.Mean)^2);

MSE.LC50.new[j]=mean((10^m.new-10^Log.True.Mean)^2);

j=j+1;

}

Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50,Mean.LC50.new,MSE.LC50.new))

dimnames(Result)=list(c("","","","",""), c("Sample.Size","Mean","MSE","Mean.Cubic","MSE.Cubic"))

{

cat(" Simulation Result","\n");

cat(" Equally Spaced X:", Eqx.or.Not,","," ", "Same Number of Subjects:", Eqn.or.Not,"\n\n")

Result

}

pp=as.matrix(cbind(z1,Pstar))

plot(z1,Pstar,xlab="Transformed Dose Level",ylab="p")

xx=c(pp[2,1],pp[3,1],pp[4,1])

yy=c(pp[2,2],pp[3,2],pp[4,2])

lines(xx,yy,col="blue",lwd=2)

xx=seq(0,1,by=0.01);

lines(xx,a0+0.5+a1*xx+a2*xx^2+a3*xx^3,col="red",lwd="2");

2.5 The Shuster-Dietrich Method

The method was recommended by Shuster and Dietrich (1976) for estimating the dose response curves in quantal response assays. It is a general inverse regression procedure.

2.5.1 Algorithm

(1). Calculate

where i=1,2,…, k, di is the dose level.

(2). The estimated of log10LC50, m, is

.

(3). The LC50 is estimated by

.

(4) The 95% confidence interval for m is . Hence the approximate 95% confidence interval for LC50 is .

2.5.2 Simulation study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.29 Case 1

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Epirical Length

n=10

2.7315

0.0264

0.918

1.231082

n=20

2.7379

0.0157

0.910

1.167366

n=30

2.7476

0.0110

0.902

1.136606

n=40

2.7462

0.0084

0.912

1.123817

n=50

2.7431

0.0079

0.802

1.090995

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.30 Case 2

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Epirical Length

n=10

2.7161

0.0296

0.916

1.244353

n=20

2.7419

0.0148

0.912

1.164231

n=30

2.7526

0.0104

0.908

1.135364

n=40

2.7478

0.0084

0.880

1.112813

n=50

2.7440

0.0080

0.812

1.091147

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.31 Case 3

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Epirical Length

n=10

4.2048

0.4309

0.916

1.7450

n=20

4.2723

0.2583

0.972

1.6843

n=30

4.2850

0.1771

0.954

1.4784

n=40

4.2927

0.1310

0.972

1.4521

n=50

4.2923

0.1105

0.878

1.2772

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.32 Case 4

Sample Size

Mean

MSE

Empirical 95% CI

Mean of Epirical Length

n=10

4.1218

0.4972

0.910

1.7299

n=20

4.2904

0.2420

0.968

1.6163

n=30

4.2979

0.1683

0.950

1.4534

n=40

4.2936

0.1309

0.972

1.4484

n=50

4.2866

0.1113

0.880

1.2864

Conclusion: When the space of log dose level is equal, the Shuster-Dietrich method estimates LC50 with very small MSEs, but it obtains much larger MSEs in the case of unequally spaced dose level. There is no significant differences between different or same number of subject at each dose level. The Shuster-Dietrich method is more suitable for the case of equally spaced dose level.

2.5.3 R-code

###############################

# Shuster-Dietrich Method #

###############################

rm(list=ls())

set.seed(987654)

Eqx.or.Not=1 # 0=No, 1=Yes

Eqn.or.Not=0 # 0=No, 1=Yes

# Define Dose Levels

# Equally Spaced Transformed Dose Levels

Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

(1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

Num.Dose.Level=length(Dose.Level);

# Dose Level Transformations

Tran.Dose.Level=log10(Dose.Level);

# Define Tolerance Distribution

Log.True.Mean=log10(mean(Dose.Level));

Log.True.Sd=log10(sd(Dose.Level));

#Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)}; # normal distribution

Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)};

# Cauchy distribution

# Generate True Response Rate

True.Res.Rate=Tol.Dist(Tran.Dose.Level);

# Generate Sample

Mean.LC50=Ep.Level= Mean.Ep.Length =MSE.LC50=rep(0,5);

j=1;

for(size in c(10,20,30,40,50))

{

n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1- Eqn.or.Not)*(size+round(3*runif(Num.Dose.Level,-1,1)))

Rep.in.Sim=500;

m=SE=rep(0,Rep.in.Sim); # Initial values of Log(LC50)

for(run in seq(Rep.in.Sim))

{

Sim.Res=rep(0,Num.Dose.Level)

for(i in seq(Num.Dose.Level))

{

Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

}

# Shuster-Dietrich Method

N=sum(n);

Y=(sqrt(n/N))*Tran.Dose.Level;

Yhat=sum((sqrt(n/N))*Y);

p=Sim.Res/n;

Z=(sqrt(n/N))*asin(sqrt(p))*180/pi;

Zhat=sum(Z*sqrt(n/N));

Syz=sum(Y*Z)-Yhat*Zhat;

Szz=sum(Z*Z)-Zhat*Zhat;

bhat=Syz/Szz;

m[run]=Yhat-bhat*(Zhat-45);

}

Syy= sum(Y^2)-Yhat^2

Sn=sqrt(820.7*((Zhat-45)^2*Syy/Szz^2+bhat^2))

SE=Sn/sqrt(N)

L.end=10^(m-1.96*SE)

R.end=10^(m+1.96*SE)

freq=sum((10^Log.True.Mean>=L.end)*(10^Log.True.Mean<=R.end))

Ep.Level[j]=freq/Rep.in.Sim;

Mean.Ep.Length[j]=mean(10^(2*1.96*SE));

Mean.LC50[j]=mean(10^m);

MSE.LC50[j]=mean((10^m-10^Log.True.Mean)^2);

j=j+1;

}

Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50,Ep.Level, Mean.Ep.Length))

dimnames(Result)=list(c("","","","",""), c("Sample.Size","Mean","MSE","Empirical.CL "," Mean.Ep.Length "))

{

cat(" Simulation Result","\n");

cat(" Equally Spaced X:", Eqx.or.Not,","," ", "Same Number of Subjects:", Eqn.or.Not,"\n\n")

Result

}

2.6 The Shuster-Yang Method

Shuster and Yang (1975) have developed a nonparametric method to estimate the minimum dose level to induce a given response rate. They made a comparison between this method and other well-known nonparametric methods to show that the Shuster-Dietrich method is better in certain conditions.

They also introduced some practical applications of this method in particular conditions. It can estimate LC50 when the dose level is discrete. Even based on a preliminary scan of the dosage levels, this method can get exact inferences about the location of LC50. Since interpolation is applied to this procedure, there is no need to make a distribution-free estimator for LC50, but it needs to construct a distribution-free statement for the LC50 interval.

This method provides an easy algorithm for locating the interval where the estimate of LC50 can be found by liner interpolation. For this method, the only assumption is that the response rate is no decreasing in a function with respect to dose level.

2.6.1 Algorithm

(1). Calculate

,

.

(2). Find the first integer s such that .

(3). LC50 is in the interval (xs, xs+1).

(4). By linear interpolation on x and p, the estimated of LC50 is

.

2.6.2 Simulation study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.33 Case 1

Sample Size

Mean

MSE

10

2.7961

0.0773

20

2.8052

0.0434

30

2.8018

0.0283

40

2.8027

0.0192

50

2.8053

0.0170

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.34 Case 2

Sample Size

Mean

MSE

10

2.803388

0.0812

20

2.806847

0.0406

30

2.797690

0.0253

40

2.804706

0.0183

50

2.804898

0.0177

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.35 Case 3

Sample Size

Mean

MSE

10

4.4294

1.6622

20

4.5022

1.0458

30

4.4032

0.6570

40

4.4610

0.5157

50

4.4382

0.4872

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.36 Case 4

Sample Size

Mean

MSE

10

4.5263

1.7224

20

4.5546

1.0631

30

4.4123

0.6190

40

4.4656

0.5010

50

4.4496

0.4850

Conclusion: Comparing the case of equally spaced dose levels, the unequal space case get significant greater MSEs. It indicates that the Shuster-Yang Method works much better in equal space case than in unequal case.

2.6.3 R-code

The simulation study is conducted by using the following R program.

########################

# Shuster-Yang Method #

########################

rm(list=ls())

set.seed(987654)

Eqx.or.Not=0 # 0=No, 1=Yes

Eqn.or.Not=0 # 0=No, 1=Yes

# Define Dose Levels

# Equally Spaced Transformed Dose Levels

Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

(1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

Num.Dose.Level=length(Dose.Level);

# Dose Level Transformations

Tran.Dose.Level=log10(Dose.Level);

# Define Tolerance Distribution

Log.True.Mean=log10(mean(Dose.Level));

Log.True.Sd=log10(sd(Dose.Level));

Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)}; # normal distribution

#Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)};

# Cauchy distribution

# Generate True Response Rate

True.Res.Rate=Tol.Dist(Tran.Dose.Level);

# Generate Sample

Mean.LC50=Mean.LC50.new=Ep.Level=Mean.Ep.Length=MSE.LC50=MSE.LC50.new=rep(0,5);

j=1;

for(size in c(10,20,30,40,50))

{

n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1-Eqn.or.Not)*(size+round(3*runif(Num.Dose.Level,-1,1)))

Rep.in.Sim=500;

LC50=rep(0,Rep.in.Sim);

for(run in seq(Rep.in.Sim))

{

Sim.Res=rep(0,Num.Dose.Level)

for(i in seq(Num.Dose.Level))

{

Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

}

# Shuster-Yang Method

d=Dose.Level

p=Sim.Res/n;

A=Sim.Res-0.5*n

S=cumsum(A)

Ss=min(S)

s=which(S==Ss)

J=function(s)

(

if (length(s)>=2)

(return(min(s)+1))

else

(return(s+1))

)

h=J(s)

LC50[run]=d[h-1]+(d[h]-d[h-1])*(0.5-p[h-1])/(p[h]-p[h-1])

}

Mean.LC50[j]=mean(LC50,na.rm=TRUE) #remove the unavailable value of LC50 because h>length of dose levers

MSE.LC50[j]=mean((LC50-10^Log.True.Mean)^2, na.rm=TRUE); #remove the unavailable value of LC50 because h>length of dose levers

j=j+1;

}

Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50))

dimnames(Result)=list(c("","","","",""), c("Sample.Size","Mean","MSE"))

{

cat(" Simulation Result","\n");

cat(" Equally Spaced X:", Eqx.or.Not,","," ", "Same Number of Subjects:", Eqn.or.Not,"\n\n")

Result

}

2.7 Maximum Likelihood Method

Fisher recommended and analyzed Maximum-likelihood estimation in 1912. This method became very popular since then. Researchers always prefer it for estimating LC50 today because of its high accuracy.

2.7.1 Algorithm

(1) For the normal model

,

where .

(2) the likelihood function is proportional to, i= 1, 2,… ,k

(3) The log-Likelihood function is ,

where .

(4) Let . The parameters and are the solution of

where , .

2.7.2 Simulation Study

Case 1. Equally spaced dose levels, same number of subject at each dose level

Table 2.37 Case 1

Sample Size

Mean

MSE

10

2.7979

0.0346

20

2.7920

0.0168

30

2.7901

0.0110

40

2.7894

0.0074

50

2.7868

0.0064

Case 2. Equally spaced dose levels, different number of subject at each dose level

Table 2.38 Case 2

Sample Size

Mean

MSE

10

2.7985

0.0356

20

2.7912

0.0162

30

2.7901

0.0106

40

2.7909

0.0075

50

2.7870

0.0066

Case 3. Unequally spaced dose levels, same number of subject at each dose level

Table 2.39 Case 3

Sample Size

Mean

MSE

10

4.5441

1.1543

20

4.4308

0.3693

30

4.3997

0.2149

40

4.3949

0.1517

50

4.3801

0.1199

Case 4. Unequally spaced dose levels, different number of subject at each dose level

Table 2.40 Case 4

Sample Size

Mean

MSE

10

4.5754

1.4268

20

4.4271

0.3420

30

4.4040

0.2080

40

4.3981

0.1535

50

4.3752

0.1205

Conclusion: Comparing the case of equally spaced dose levels, the equal space case get similar MSEs. In addition, they give the smallest MSEs compared to other methods, indicating that the Maximum Likelihood method works much better than other methods.

2.7.3 R-code

######################

# Maximum Likelihood    #

######################

  rm(list=ls())

  set.seed(987654)

  Eqx.or.Not=0      # 0=No, 1=Yes

  Eqn.or.Not=0     # 0=No, 1=Yes

# Define Dose Levels

  Dose.Level=Eqx.or.Not*c(10^0.1, 10^0.2, 10^0.3, 10^0.4, 10^0.5, 10^0.6, 10^0.7)+

                 (1-Eqx.or.Not)*c(1,2,3,3.5,5,7,9);

 

   Num.Dose.Level=length(Dose.Level);

 

 # Dose Level Transformations

 

   Tran.Dose.Level=log10(Dose.Level);

 

 # Define Tolerance Distribution

 

   Log.True.Mean=log10(mean(Dose.Level));

   Log.True.Sd=log10(sd(Dose.Level));

   Tol.Dist=function(u){pnorm(u,Log.True.Mean,Log.True.Sd)};

 # normal distribution

#Tol.Dist=function(u){pcauchy(u, location=Log.True.Mean, scale, scale=Log.True.Sd)}; # Cauchy distribution

 # Generate True Response Rate

 

   True.Res.Rate=Tol.Dist(Tran.Dose.Level);

 

 # Generate Sample

 

   Mean.LC50=Mean.LC50.new=Ep.Level=MSE.LC50=MSE.LC50.new=rep(0,5);

   j=1;

   for(size in c(10,20,30,40,50))

     {

      

      n=(Eqn.or.Not)*rep(size,Num.Dose.Level)+(1-Eqn.or.Not) * ( size+round( 3*runif(Num.Dose.Level,-1,1)))

    

      #n=rep(size,Num.Dose.Level); # Population Size at Each Dose Level

      Rep.in.Sim=500;

      m=m.new=rep(0,Rep.in.Sim);)

       for(run in seq(Rep.in.Sim))

        {

         Sim.Res=rep(0,Num.Dose.Level)

         for(i in seq(Num.Dose.Level))

               {

                 Sim.Res[i]=rbinom(1, n[i], True.Res.Rate[i])

               }

 

          # Maximum Likelihood Method

          X=Tran.Dose.Level;

          RRes=Sim.Res/n

         

          myglm=glm(RRes~X, family=binomial(link=probit))

         

          m[run]=-myglm$coefficients[1]/myglm$coefficients[2]

                   }

        Mean.LC50[j]=mean(10^m)

        MSE.LC50[j]=mean((10^m-10^Log.True.Mean)^2);

 

        j=j+1;

        }

             

        Result=as.matrix(cbind(seq(5)*10,Mean.LC50,MSE.LC50))

        dimnames(Result)=list(c("","","",""),("Sample.Size","Mean","MSE"))

        Result

Conclusion and Discussion

Case 1 di=xi+1-xi is constant, ni=50

Case 2 di is a constant, ni is not a constant.

Case 3 di is not a constant, ni =50.

Case 4 di is not a constant, ni is not a constant

From the data in the above table, one can see the Modified Dragstedt-Behrens method and Maximum likelihood perform the best in case of equally spaced log dose levels. In addition, the Modified Dragstedt-Behrens method and Modified Reed-Muench method outperform others when the number of objects in each dose level is constant. However, in practice, we never know the exactly real model. The best data for estimation should be sufficiently checked before fitting the parametric model. Different data demand different methods.

References Or Bibliography

Jonathan J. Shuster and Mark C. K. Yang (1975). A Distribution-Free Approach to Quantal Response Assays. The Canadian Journal of Statitics, 3(1), 57-70.

L. J. Reed and H. Muench (1938). A simple method of estimating fifty per cent endpoints. The American Journal of Hygiene. 27(3), 493-497.

Gabrielle E. Kelly (2001). The median lethal dose- design and estimation. The Statistician. 50(1), 41-50.

J. T. Litchfield, JR. and F. Wilcoxon (1948). A simplified method of evaluating dose-effect experiments. Journal of Pharmacology and Experimental Therapy. 96, 99-113.

Robert D. Bruce (1985). An up-and-down procedure for acute toxicity testing. Fundamental and Applied Toxicology, 5, 151-157.

B. M. Bennety (1952). Estimation of LD50 by moving averages. The Journal of Hygiene, 50(2), 157-164.

J. A. Hoekstra. (1991). Estimation of the LC50, a review. Environment, 2(2), 139-152.

David J. Finney (1978). Statistical method in Biological assay. 3rd ed. Charles Griffin & Company LTD, London and High Wycombe.

David J. Finney (1971). Probit Analysis. 3rd ed. Cambridge University. Press, Cambridge, England.


Request Removal

If you are the original writer of this essay and no longer wish to have the essay published on the UK Essays website then please click on the link below to request removal:

Request the removal of this essay


More from UK Essays