Because the GMO Panel of the EFSA scientific opinion on statistical consideration for the safety evaluation of GMO’s requires that the mean difference between the GM crop and its conventional counterpart should be within the natural variability of available commercially varieties. We proposed a new criterion to achieve this goal. In addition, equivalence hypothesis was formulated based on the linearization of proposed criterion. A statistical inference procedure was proposed to test the
equivalence hypothesis. This method rejects the null hypothesis in Equation (3.4) at the
significance level if the 100(1)% upper confidence bound constructed by the
MLS method is less than 0.
Although the simulation studies show that our proposed method can adequately control the size for moderate to large number of sites and replicates, the power increases as the variance of the commercial varieties increases. In other words, equivalence can easily be established when the variance of commercial varieties becomes large. This shortcoming of the proposed criterion can be resolved if we can restrict the variance of commercial varieties to an upper limit
0
2
g . In other words, the scaled criterion can be revised as
74
The performance of the revised criterion required further research in terms of the choice of 0
2
g , size and power.
On the other hand, one of the drawbacks of the EFSA method is that both the confidence interval for the mean difference between the GM crop and its conventional counterpart and the adjusted equivalence limits given in Equation (2.4) are random intervals. Therefore, the probability of a random interval within another random interval is difficult to evaluate. In addition, the adjusted equivalence interval becomes narrow as the number of sites or the number of replicates increases. It turns out, as shown by the simulation results, that the empirical size decreases as the number of site or the number of replicate increases. This is a counter-intuitive phenomenon which requires further research.
The EFSA statistical considerations also suggest that at least 6 commercial varieties should be compared concurrently in the same field trial using a randomized completely block design. The EFSA statistical consideration also provide an example of design for multiple crops with 4 commercial varieties. Since the block size increases as the numbers of GM crops, their conventional counterparts, and commercial varieties increase, then heterogeneity within each block also increases. Therefore, we suggest that
75
the experiment should be restricted to one GM crop, its conventional counterpart, and 6 commercial varieties for a total of 8 plots per block. For the situation where multiple GM crops with more than 6 commercial varieties are compared in a field trial at a site, a balanced incomplete block design should be used.
As mentioned before, our research is confined only to one variable. However for assessment of equivalence between the GM crop and its counterparts consists of many variables such as agronomic, phenotypes and compositional characteristics, usually exceeding 25. If the requirement for approved of a GM crop is that all characteristics must meet equivalence criterion at the 5% significance level, then its corresponding overall type Ⅰ error rate is 3.0 10 33 which is extremely low assuming that all endpoints are independent. This is a problem of multiple correlated must-win comparisons (Julious and Mclntyre, 2012). Further research in this area for evaluation of equivalence for multiple endpoints between the GM crop and its counterpart is required.
76
References
Altman D. and Bland J.M. (1995) Absence of evidence is not evidence of absence.
British Medical Journal, 311:485
Altman D. and Bland J.M. (2004) Confidence intervals illuminate absence of evidence.
British Medical Journal, 328:1016-1017
Chow S. C., Liu J. P. (2010) Design and Analysis of Bioavailability and Bioequivalence Studies, 3th Ed., CRC/Chapman and Hall, Taylor and France, New York.
Chow S. C., Shao J., Wang H. (2003) Individual bioequivalence testing under 2x3 designs. Statistics in Medicine 2002; 21:629-648
European Food Safety Authority (EFSA) (2010). EFSA GMO Panel opinion on Statistical considerations for the safety evaluation of GMOs. European Food Safety Authority, Parma, Italy
European Food Safety Authority (EFSA) (2001) Directive 2001/18/EC. European Food Safety Authority, Parma, Italy
European Food Safety Authority (EFSA) (2010) EFSA GMO Panel opinion on Statistical considerations for the Safety Evaluation of GMOs. Presentation of the full statistical code used for the example. European Food Safety Authority, Parma, Italy
European Agency for the Evaluation of Medicinal Products (EMEA). 2001. Note for guidance on the investigation of bioavailability and bioequivalence. Document CPMP/EWP/QWP/1401/98, European Agency for the Evaluation of Medicinal Products, Committee for Proprietary Medicinal Products, London, UK.
Food and Drug Administration (FDA). (2001). Guidance for Industry - Statistical approaches to establishing bioequivalence. U.S. Department of Health and Human Services, Center for Drug Evaluation and Research. Food and Drug Administration, Rockville, Maryland, 2001.
Hyslop T., Hsuan F., Holder D.J. (2000) Asmall sample confidence interval approach to assess individual bioequivalence. Statistics in Medicine; 19:2885-2897
Hsu, J.C., Hwang, J.T.G., Liu, H.K., and Ruberg, S.J. (1994) Confidence intervals
77
associated with tests for bioequivalence. Biometrika, 81, 103-114.
James, C. (2012) Global Status of Commercialized Biotech/GM Crops: 2012. ISAAA Brief No. 44. ISAAA: Ithaca, NY
Julious, S. A., McIntyre, N. E.. (2012) Sample sizes for trials involving multiple correlated must-win comparisons. Pharmaceutical Statistics, 11 (2):177-85
Kenward M.G. and Roger J.H. 1997. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53: 983-997.
Lee, A.F.S. and Fineberg, N.S. 1991. A fitted test for the Behrens-Fisher problem.
Communications in Statistics-Theory and Methods, 20, 653-666.
Lee Y.H, Shao J. and Chow S.C. (2004). Modified large-sample confidence intervals for linear combinations of variance components: extension, theory, and application.
Journal of the American Statistical Association, 99, 467-478
Patterson S.D., Jones B. (2012) Observations on scaled average bioequivalence.
Pharmaceutical Statistics. 11,1-7
Schuirmann D.J. (1987) A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability. Journal of Pharmacokinetics and Biopharmaceutics, 15: 657-680.
Tang F.W. (2009) A study on application of Hotelling T-squared to statistical evaluation of genetically modified products. Master thesis, National Taiwan University
Tsai P.Y. (2009) A study on statistical evaluation of in vitro bioequivalence with consideration of variance components. Master thesis, National Taiwan University World Health Organization (WHO) (2005) World Health Organization Draft Revision
on Multisource (Generic) Pharmaceutical Products: Guidelines on Registration Requirements to Establish Interchangeability, Geneva, Switzerland.
World Health Organization (WHO) (2002) 20 questions on genetically modified foods, Geneva, Switzerland.
78
Appendix
SAS Codes for Simulation Program
/*Dataset simulation*/
data aaa;
do site=1 to 4;
site_value=rand("normal",0,sqrt(3));
do rep=1 to 4;
rep_value=rand("normal",0,sqrt(1));
do k=1 to 8;
if k=1 then treat= -0.073;
else if k=2 then treat=0.15;
else treat=-0.077;
output;
end;
end;
end;run;
79
data bbb;
do k=1 to 8;
if k=1 then genotype=0;
else if k=2 then genotype=0;
else genotype=rand("normal",0,sqrt(2));
do rep=1 to 4;
output;
end;
end;
proc sort data=aaa; by rep k;
proc sort data=bbb; by rep k;
data ccc; merge aaa bbb; by rep k;
e=rand("normal",0,sqrt(0.6));
m=20;
y=site_value+rep_value+treat+genotype+e+m;
proc sort data=ccc; by site rep;
run;
proc sql;
80
create table ln as
select site, site_value, rep, rep_value, k, treat, genotype, e, m, log(y) as y from ccc;
run;
/*End of dataset simulation*/
/*Calculation for mean of conventional counterparts*/
proc sql;
create table ccon as select mean(y) as y1mean from ln
where k eq 1;
run;
/* Calculation for mean of GMOs */
proc sql;
create table gmo as select mean(y) as y2mean
81
from ln where k eq 2;
run;
/*Calculation for variance of between GMO and its conventional counterpart*/
proc sql;
create table ccon_y as select y as y1
from ln where k eq 1;
run;
proc sql;
create table gmo_y as select y as y2
from ln where k eq 2;
run;
data y12;
82
merge ccon_y gmo_y;
run;
proc sql;
create table y12var as select var(y1-y2) as var_y12 from y12;
run;
/* Logarithmic transformation of dataset*/
data ss;
set ln;
if k eq 1 then genotypegroup="comp";
else if k eq 2 then genotypegroup="gmo";
else genotypegroup="ref";
run;
/*Definition of variance components*/
data sss;
set ss;
83
if (k eq 1) or (k eq 2) then gmo_genotype=0;
if k eq 3 then gmo_genotype=1;
if k eq 4 then gmo_genotype=2;
if k eq 5 then gmo_genotype=3;
if k eq 6 then gmo_genotype=4;
if k eq 7 then gmo_genotype=5;
else if k eq 8 then gmo_genotype=6;
run;
/*Analysis of varieties for fixed effects*/
proc varcomp method=type1 data=sss;
class site rep genotypegroup gmo_genotype;
model y=genotypegroup site rep(site) gmo_genotype / fixed=1;
ods output Estimates=est_varcomp ANOVA=ddd;
run;
data comvarcomp;
set est_varcomp(where=(varcomp="Var(gmo_genotype)"));
comvarcomp=Estimate;
84
keep comvarcomp;
run;
/*Capture for required data for MLS method*/
data eee (keep= source MS);
set ddd ; run;
proc transpose data=eee
out=eee_transposed name=MS prefix=MS;
id source;
run;
data fff (keep= source DF);
set ddd;
run;
proc transpose data=fff
out=fff_transposed name=DF prefix=DF;
id source;
run;
85
data ggg;
merge fff_transposed eee_transposed ccon gmo comvarcomp y12var;
run;
/*Construction of 95% upper bound by MLS method*/
data hhh;
set ggg;
n=4;
p=4;
t=3;
q=6;
s=(t-1)+q;
diff=abs(y2mean-y1mean);
A=diff**2;
tinv=tinv(0.95,2*(n*p-1));
/*var_diff=2*(MSsite+(p-1)*MSrep_site_+p*(s-1)*MSError)/(n*p*p*s);*/
ucl=(diff+tinv*sqrt((var_y12)/(n*p)))**2;
C=(ucl-A)**2;
86
theta=0.2;
lamda1=-theta*MSgmo_genotype/(n*p);
lamda2=theta*MSError/(n*p);
B=-theta*comvarcomp;
chi95=cinv(0.95,q-1);
chi5=cinv(0.05,DFError);
D=(lamda1**2)*((q-1)/chi95-1)**2+(lamda2**2)*(DFError/chi5-1)**2;
CL=A+B+sqrt(C+D);
proc print data=hhh;run;
---EFSA guidance of statistical consideration--- data s;
set ln;
if k eq 1 then genotypegroup="comp";
if k eq 2 then genotypegroup="gmo";
if k ge 3 then genotypegroup="ref";
if genotypegroup eq "ref" then indref=1;
else indref=0;
/*Mixed model analysis for differences*/
87
proc mixed data=s CL=WALD;
class site site_value rep rep_value k treat genotype e m genotypegroup y;
model y=genotypegroup/s covb outp=out ddfm=kenwardroger;
random site_value site_value*rep_value indref*genotype;
estimate "gmo_comp" genotypegroup -1 1 0/CL alpha=0.1;
lsmeans genotypegroup;
ods output lsmeans=lsmeans_g estimates=estdif covparms=covparms;
run;
/*Geometric means*/
data lsmeans_g2;
set lsmeans_g;
gm=exp(estimate);
keep y genotypegroup gm;
run;
/*Assessment of differences*/
data estdif1;
set estdif;
ratiod=exp(estimate);lowd=exp(lower);uppd=exp(upper);
88
dgc=estimate;sed=stderr;dfd=df;
lsddif=sed*tinv(0.95,dfd);
keep y dgc ratiod lowd uppd sed dfd lsddif;
run;
/*Mixed model for setting equivalence limits*/
proc mixed data=s CL=WALD;
class site_value rep_value genotype genotypegroup;
model y=genotypegroup /s covb outp=out ddfm=kenwardroger;
random site_value site_value*rep_value genotype;
estimate "gmo_ref" genotypegroup 0 1 -1/CL alpha=0.05;
ods output estimates=esteq covparms=covparms;
run;
/*Keep estimated variance components for genotype*/
data varg;
set covparms(where=(covparm="genotype"));
varg=estimate;
89
keep varg;
run;
/*95% equivalence limits ratio GMO/ref*/
data esteq1;
merge esteq varg;
ratioe=exp(estimate);lowe=exp(lower);uppe=exp(upper);
dgr=estimate;see=stderr;dfe=df;
lsdeq=see*tinv(0.975,dfe);
keep varg dgr ratioe lowe uppe see dfe lsdeq;
run;
/*Mixed model analysis for testing non-equivalence*/
proc mixed data=s CL=WALD;
class site_value rep_value genotype genotypegroup;
model y=genotypegroup/s covb outp=out ddfm=kenwardroger;
random site_value site_value*rep_value indref*genotype;
estimate "gmo_ref" genotypegroup 0 1 -1/CL alpha=0.1;
90
ods output estimates=estteq;
run;
/*Assessment of equivalences(adjusted scale)*/
data estteq1;
merge estdif1 esteq1 estteq;
lsdteq=stderr*tinv(0.95,df);
fac=lsddif/lsdteq;
lower=-dgr-lsdeq;upper=-dgr+lsdeq;
adjlow=dgc+fac*lower;adjupp=dgc+fac*upper;
adjlow=exp(adjlow);adjupp=exp(adjupp);
keep dgc stderr df lsdteq fac adjlow lowd uppd adjupp varg;
run;
/*Classification of results*/
data classification;
set estteq1;
tol=1e-6;
91
noequivtest=(varg<=tol);
sigdif=((lowd>1) or (uppd<1));
equiv=((lowd>=adjlow) and (uppd<=adjupp));
warning=((exp(dgc)<adjlow) or (lowd>adjupp));
nonequiv=((uppd<adjlow) or (lowd>adjupp));
if ((not sigdif) and equiv) then type=1;
if ( sigdif and equiv) then type=2;
if ((not sigdif) and (not equiv) and (not noequivtest) and (not warning)) then type=3;
if ( sigdif and (not equiv) and (not noequivtest) and (not warning)) then type=4;
if ((not sigdif) and (not equiv) and (not noequivtest) and warning) then type=5;
if ( sigdif and (not noequivtest) and warning and nonequiv) then type=6;
if ( noequivtest) then type=7;
keep noequivtest sigdif equiv warning nonequiv type;
proc sort data=classification;
by type;
proc print data=classification;
run;
---End of EFSA Guidance of Statistical Consideration---
92
---Size/Power Simulation--- data aaa;
do n=10 to 100;
do p=10 to 100;
do q=6,12,18,24;
t=3;
s=(t-1)+q;
dfe=(n*p*s-1)-(n-1)-(n*p-n)-(t-1)-(q-1);
delta=0; /*let delta=0.223 for size simulation*/
sigma_site=2;
sigma_rep=1;
sigma_genotype=0.2; /*let sigma_genotype=0.248645 for size simulation*/
sigma_e=0.6;
theta=0.2;
t=quantile("t",0.95,2*(n*p-1));
q1=((q-1)/quantile("chisq",0.95,(q-1))-1)**2;
q2=(dfe/quantile("chisq",0.05,dfe)-1)**2;
do j=1 to 10000;
93
diff=rand("normal",delta,sqrt(2*(sigma_site+sigma_rep+sigma_e)/(n*p)));
MSgenotype=rand("chisq",q-1)*(sigma_e+n*p*sigma_genotype)/(q-1);
MSe=rand("chisq",dfe)*sigma_e/dfe;
A=diff**2;
u=(abs(diff)+t*sqrt(2*(sigma_site+sigma_rep+sigma_e)/(n*p)))**2;
C=(u-A)**2;
B=-theta*(MSgenotype-MSe)/(n*p);
lamda1=-theta*MSgenotype/(n*p);
lamda2=theta*MSe/(n*p);
D=lamda1**2*q1+lamda2**2*q2;
CL=A+B+sqrt(C+D);
if CL<0 then type=1; else type=0;
output;
end;end;end;end;
proc means data=aaa mean; by p;
output out=result mean(type)=;
run;
---End of Power Simulation---
94
---Response Surface--- goptions reset=global gunit=pct border
ftext=swissb htitle=6 htext=3;
data one;
input site replicate power ; cards;