EPA/600/R-XX/XXX | October 2010 | www.epa.gov/nrmrl

Risk Assessment of the Evolution of 

Resistance in European Corn Borer(Ostrinia 

nubilalis Hubner) and Southwestern Corn 

Borer(Diatraea grandiosella) to Smartstax® corn 

Letter Report 

Sustainable Technology Division, Cincinnati, OH 45268 

National Risk Management Research Laboratory 

Office of Research and Development
EPA/600/R-XX/XXX 

October 2010

Risk Assessment of the Evolution of 

Resistance in European Corn Borer(Ostrinia 

nubilalis Hubner) and Southwestern Corn 

Borer(Diatraea grandiosella) to Smartstax® corn 

Letter Report 

October 2010 

 

Michael Caprio1, John A. Glaser2 

1

Mississippi State University, Department of Entomology and Plant Pathology, Mississippi 

State, MS 39762-9775 

2

US Environmental Protection Agency, National Risk Management Research Laboratory, 

Sustainable Technology Division, Cincinnati, OH 45268 

Sustainable Environments Branch 

Sustainable technology Division 

National Risk Management Research Laboratory 

Office of Research and Development 

Cincinnati, Ohio 45258 

U.S. Environmental Protection Agency 

Office of Research and Development 

Washington, DC 20460 

2
Notice

The U.S. Environmental Protection Agency through its Office of Research and Development 

funded and managed the research described here. It has been subjected to the Agency=s peer 

and administrative review and has been approved for publication as an EPA document. Mention 

of trade names or commercial products does not constitute endorsement or recommendation for 

use 

. 

The information provided in this paper is not intended for determining compliance with state or 

federal regulations for Plant Iincorporated Pesticide resistance management. For these purposes, 

information should be obtained from the appropriate local, state and federal regulatory agencies, 

the Federal Insecticide Fungicide and Rodenticide Act, the United States Code, the Code of 

Federal Regulations, the Federal Register and other sources. 

3
Abstract 

Genetically-Modified(GM) crops containing pesticidal traits are valued as environmental 

assets since they directly replace broad spectrum chemical pesticides having detrimental 

environmental and human health effects. As part of the risk assessment for each new GM 

crop registration requires an evaluation of the sustainable lifetime for each new GM crop 

variety. Models employing principles of population genetics have been used to simulate field 

conditions where the development of insect pest resistance to the control traits of the crop is 

encountered. The Monsanto Company’s submission of the Superstax® corn event for control 

of the European corn borer(Ostrinia nubilalis Hubner) and the Southwestern Corn Borer 

(Diatraea grandiosella) has been evaluated for the development of insect pest resistance 

using the EPA developed deterministic two locus simulation model POPGEN D2L. The design 

of research is presented with the interpreted results. 

4
Foreword 

The U. S. Environmental Protection Agency (USEPA) is charged by Congress with 

protecting the Nation's land, air, and water resources. Under a mandate of national 

environmental laws, the Agency strives to formulate and implement actions leading to a 

compatible balance between human activities and the ability of natural systems to support and 

nurture life. To meet this mandate, USEPA’s research program is providing data and technical 

support for solving environmental problems today and building a science knowledge base 

necessary to manage our ecological resources wisely, understand how pollutants affect our 

health, and prevent or reduce environmental risks in the future. 

The National Risk Management Research Laboratory (NRMRL) is the Agency's center for 

investigation of technological and management approaches for preventing and reducing risks 

from pollution that threaten human health and the environment. The focus of the Laboratory's 

research program is on methods and their cost-effectiveness for prevention and control of 

pollution to air, land, water, and subsurface resources; protection of water quality in public water 

systems; remediation of contaminated sites, sediments and ground water; prevention and control 

of indoor air pollution; and restoration of ecosystems. NRMRL collaborates with both public and 

private sector partners to foster technologies that reduce the cost of compliance and to anticipate 

emerging problems. NRMRL's research provides solutions to environmental problems by: 

developing and promoting technologies that protect and improve the environment; advancing 

scientific and engineering information to support regulatory and policy decisions; and providing 

the technical support and information transfer to ensure implementation of environmental 

regulations and strategies at the national, state, and community levels. 

This publication has been produced as part of the Laboratory's strategic long-term 

research plan. It is published and made available by USEPA's Office of Research and 

Development to assist the user community and to link researchers with their clients. 

Sally Gutierrez, Director 

National Risk Management Research Laboratory 

5
Contents 

Notice............................................................................................................................... 

Abstract............................................................................................................................ 

Foreword.......................................................................................................................... 

Contents............................................................................................................................ 

List of Figures..................................................................................................................... 

List of Tables....................................................................................................................... 

Terminology.................................................................................................................... 

CHAPTER 1 Introduction............................................................................... 

1.1 Background.................................................................................. 

1.2 Purpose....................................................................................... 

CHAPTER 2 Study design............................................................................ 

 2.1 RM Models 

 2.2 Parameters 

 2.3 Methods and materials 

 2.4 Results 

CHAPTER 3 Discussion and Conclusions 

References....................................................................................................................... 

Appendix 1 Resistance Management Models Used for this Analysis................................ 

Appendix 2 Deterministic Parameter Description.............................................................. 

Appendix 3 European CornBorer Sensitivity Analysis 

Appendix 4 Southwestern Corn Borer Sensitivity Analysis 

6
List of Figures 

Figure 1: Schematic of the deterministic frequency based model for both seed-blend and 

block refuge deployments 

7
List of Tables 

Table 1: PERT parameters for the variables in the structured and larval movement simulations 

of European corn Borer 

Table 2: PERT parameters for the variables in the structured and larval movement simulations 

of Southwestern corn Borer 

Table 3: The fitnesses and dominance values drawn in this example of fitness calculations 

Table 4a: The fitness values calculated for Cry1A.105 and Cry2Ab2 genotypes when the 

Cry1F genotype was SS 

Table 4b: The fitness values calculated for Cry1A.105 and Cry2Ab2 genotypes when the 

Cry1F genotype was RS 

Table 4c: The fitness values calculated for Cry1A.105 and Cry2Ab2 genotypes when the 

Cry1F genotype was RR 

8
Risk Assessment of the Evolution of 

Resistance in European Corn Borer(Ostrinia 

nubilalis Hubner) and Southwestern Corn 

Borer(Diatraea grandiosella) to Smartstax® corn 

Letter Report 

October 29, 2010 

 

Michael Caprio1, John A. Glaser2 

1

Misssissippi State University, Department of Entomology and Plant Pathology, Mississippi 

State, MS 39762-9775 

2

US Environmental Protection Agency, National Risk Management Research Laboratory, 

Sustainable Technology Division, Cincinnati, OH 45268 

Keywords: PIP corn variety, PIP crop registration, resistance management, simulation modeling, refuge in bag

CHAPTER 1 Introduction 

The USEPA, under its administration of the Federal Insecticide Fungicide and Rodenticide 

Act (FIFRA), requires the registration of all pesticides and pesticidal materials. The GM or 

Plant Incorporated Protectant (PIP) crops containing pesticidal traits are subject to FIFRA 

registration requirements. Concerns relating to future environmental effects of these crops 

can be investigated only using simulation models. Model application in this instance offers an 

analysis of the useful lifetime for the pesticidal traits. USEPA has clearly supported the use of 

these traits as replacements of broad spectrum pesticides that have significant environmental 

footprints. 

1.1 Background

Resistance management (RM) simulation models are designed as analytical methodologies 

to provide realistic assessment of the risk of resistance evolution given the allele frequencies 

used as initial conditions. New concepts of PIP crop deployment such as new trait structures 

9
and refuge-in-the-bag require a detailed evaluation to ensure the claims associated with each 

concept offer the sustainable protection of the crop. Simulation models offer a means to 

evaluate the threat of crop loss through the development of resistance in the near future. 

1.2 Purpose

This research effort was designed to evaluate the relative merit of resistance management 

concepts as submitted by Monsanto Company relative to the Superstax® corn registration. 

This research specifically analyzed the relative risks related to the control of western corn 

rootworm in seed mixture scenarios 

CHAPTER 2 Study Design 

To understand the relative risks related to the evolution of resistance in western corn 

rootworm for the Superstax® corn variety as submitted by Monsanto Company, a two locus 

frequency-based deterministic model was developed and evaluated to simulate the 

development of insect pest resistance.. 

2.1 RM Models 

POPGEN-D2L 

The deterministic model POPGEN-D2L is a modification of a previously developed simple 

frequency-based model that incorporates noncompliance, non-random mating in compliant 

fields and larval movement. This dual locus simulation model designed to evaluate the 

population dynamics of crop pest insects relating to the development of resistance to 

pesticidal controls installed in the crop. The model includes a switch to evaluate a blockrefuge design and a seed-mixture design. Biological parameters were primarily derived from 

10
that model with the exception of mortality rates, block size, dispersal rates and spatial aspects 

of refuge placement and movement.

2.2 Materials and Methods 

 The same deterministic risk assessment model was used for these simulations as was 

used for the western corn rootworm simulations (though all three loci were used). An excerpt 

from previous research is also included at the end of this report that more fully presents 

details of the deterministic model (Caprio &Glaser 2010a&b). Parameters were adjusted for 

the unique biology of these insects and the impacts of the three toxins on each species (see 

Tables 1 and 2). Some specific adjustments included an increase in the effective refuge size 

because these insects are more mobile and a more even deposition of eggs across the field 

would be expected (see Caprio & Glaser 2010a&b). The base larval and asymmetrical 

movement rates were both increased. The proportion of consistent non-compliance was 

maintained at the same level as the corn rootworm simulations, reflecting a concern over 

areas of concentrated non-compliance, but the amount of mixing between these fields and 

other fields was increased to reflect the greater overwintering dispersal of the two corn borer 

species. 

 The greatest difficulty in simulating the time to resistance based on the submitted data 

is partitioning the non-multiplicative interactions between the loci. If the mortality caused by 

each toxin were independent of the mortality caused by the other toxins, one would expect 

the overall mortality of the product to be the product of the survivorship to each toxin. If two 

toxins each kill 95% of the population, then one would expect the second toxin to kill 95% of 

the 5% of the population that survived the first toxin and the joint survivorship should be 

0.25%. These simulations can estimate this non-multiplicative interaction by comparing the 

11
overall survivorship drawn from our PERT distribution with product of the individual toxin 

mortalities based on the draws from their respective PERT distributions. In a two locus 

simulation, there is only one interaction, and all of the interaction can be attributed to that pair 

of loci. In three locus simulations as used with Superstax®, there are three interactions 

(Cry1A.105<->Cry2Ab2, Cry1A.105<->Cry1F, Cry2Ab2<->Cry1F), and it is impossible to 

accurately distribute the interactions without additional data on the two locus interactions. In 

the absence of data the mortality effects of pairs of toxins, the attempt was made to distribute 

the interactions evenly between the three pairs of loci. Note that we assumed if an insect is 

resistant to two loci, then the mortality rates for the remaining locus are the same as if that 

toxin were used alone. The mortality of W[SS][RR][RR] is therefor estimated as W[SS]. The 

overall survivorship is given by W[SS][SS][SS]. The fitnesses that are difficult to determine based 

on the submitted data are: W[SS][SS][RR], W[RR][SS][SS] and W[SS][RR][SS]. Our approach was to 

estimate the individual interactions by taking the cube root of the total non-multiplicative 

interaction (INI, measured as Overallsurvivorship/(W[SS][RR][RR] * W[RR][SS][RR] * W[RR][RR][SS]). The 

fitness of W[SS][SS][RR] was then W[SS][RR][RR] * W[RR][SS][RR] * INI. 

 nonAdditiveInteraction = OverallSurvivorship/(fit1[0] * fit2[0] *fit3[0]); 

 nonAdditiveInteraction2L = Math.exp(Math.log(nonAdditiveInteraction)/3); 

 Survivorships[0][0][0] = OverallSurvivorship; 

 Survivorships[2][0][0] = fit1[0] * fit2[0] * nonAdditiveInteraction2L; 

In this code fit1[0] is W[SS][RR][RR] and fit2[0] is W[RR][SS][RR]. There were several additional 

constraints on the new fitness. The fitness of the double susceptible genotype had to be less 

than or equal to either of the component single susceptible genotypes. This ensures that the 

12
fitness of a genotype does not decrease as additional resistance alleles are added. 

 if (Survivorships[2][0][0] > Math.min(fit2[0],fit3[0])) 

 { Survivorships[2][0][0] = Math.min(fit2[0],fit3[0]); } 

Finally, the new fitness of the double susceptible genotype could not be less than the overall 

fitness. Note that the overall survivorship was already constrained to be at least as low as the 

lowest single toxin fitness. 

 if(Survivorships[2][0][0] < OverallSurvivorship) { 

 Survivorships[2][0][0] = OverallSurvivorship; } 

This process was repeated for all three double susceptible genotypes. The remaining 

fitnesses involving heterozygous genotypes could then be calculated based on dominance 

values and fitness values already estimated. The relevant code is (hl1, hl2 and hl3 are the 

dominance values drawn for dominance of the Cry1A.105, Cry2Ab2 and Cry1F resistance 

loci): 

 Survivorships[0][2][2] = fit1[0]; //This is the equivalent of W[SS][RR][RR] 

 Survivorships[1][2][2] = (1.0 - fit1[0]) * hl1 + fit2[0]; 

 Survivorships[2][2][2] = 1.0; 

 Survivorships[2][0][2] = fit2[0]; 

 Survivorships[2][1][2] = (1.0 - fit2[0]) * hl2 + fit2[0]; 

 Survivorships[2][2][0] = fit3[0]; 

13
 Survivorships[2][2][1] = (1.0 - fit3[0]) * hl3 + fit3[0]; 

 Survivorships[0][1][2] = (Survivorships[0][2][2] - Survivorships[0][0][2]) * hl2 + 

Survivorships[0][0][2]; 

 Survivorships[1][0][2]= (Survivorships[2][0][2] - Survivorships[0][0][2]) * hl1 + 

Survivorships[0][0][2]; 

 Survivorships[2][1][0] = (Survivorships[2][2][0] - Survivorships[2][0][0]) * hl2 + 

Survivorships[2][0][0]; 

 Survivorships[2][0][1] = (Survivorships[2][0][2] - Survivorships[2][0][0]) * hl3 + 

Survivorships[2][0][0]; 

 Survivorships[0][1][0] = (Survivorships[0][2][0] - OverallSurvivorship) * hl2 + 

OverallSurvivorship; 

 Survivorships[1][2][0] = (Survivorships[2][2][0] - Survivorships[0][2][0]) * hl2 + 

Survivorships[0][2][0]; 

 Survivorships[1][0][0] = (Survivorships[2][0][0] - OverallSurvivorship) * hl1 + 

OverallSurvivorship; 

 Survivorships[2][2][1] = (1.0 - Survivorships[2][2][0]) * hl3 + Survivorships[2][2][0]; 

 Survivorships[2][1][1] = (Survivorships[2][1][2] - Survivorships[2][1][0]) * hl3 + 

Survivorships[2][1][0]; 

 Survivorships[1][0][1] = (Survivorships[1][0][2] - Survivorships[1][0][0]) * hl3 + 

Survivorships[1][0][0]; 

 Survivorships[0][0][1] = (Survivorships[0][0][2] - OverallSurvivorship) * hl3 + 

OverallSurvivorship; 

14
 Survivorships[0][2][1] = (Survivorships[0][2][2] - Survivorships[0][2][0]) * hl3 + 

Survivorships[0][2][0]; 

 Survivorships[1][2][1] = (Survivorships[1][2][2] - Survivorships[1][2][0]) * hl3 + 

Survivorships[1][2][0]; 

 Survivorships[0][1][0] = (Survivorships[0][2][0] - OverallSurvivorship) * hl2 + 

OverallSurvivorship; 

 Survivorships[0][1][1] = (Survivorships[0][1][2] - Survivorships[0][1][0]) * hl3 + 

Survivorships[0][1][0]; 

 Survivorships[1][1][0] = (Survivorships[2][1][0] - Survivorships[0][1][0])*(hl1+hl2)/2 + 

Survivorships[0][1][0]; 

 Survivorships[1][1][1] = (Survivorships[2][1][1] - Survivorships[0][1][1])*(hl1+hl2)/2 + 

Survivorships[0][1][1]; 

 Survivorships[1][1][2] = (Survivorships[2][1][2] - Survivorships[0][1][2])*(hl1+hl2)/2 + 

Survivorships[0][1][2]; 

As noted in the CRW modeling report, there can be multiple ways to estimate double 

heterozygotes (Caprio & Glaser, 2010c). For example, the fitness W[RS][RS][RR] could be 

calculated as (W[RR][RS][RR] – W[SS][RS][RR]) * hl1 + W[SS][RS][RR], or it could be calculated as 

(W[RS][RR][RR] – W[RS][SS][RR]) * hl2 + W[RS][SS][RR]. Another alternative is to use a mean 

dominance for double heterozygotes. Other options certainly exist. For example, the 

minimum, mean or maximum of the two calculated heterozygotes fitnesses could be used 

depending on the degree of conservatism one desires. The options became clearer as the 

modeling research progressed, and we were not entirely consistent in their use. Testing the 

15
two options suggested the model was insensitive to each option, and given the short time 

frame we chose to proceed with the calculations as shown. 

All simulations were extended for at least 1000 generations. Assuming that each gene should 

ideally last at least 15 years, the target time horizon for a three gene product would be 45 

years. Conservatively assuming that there are 3 generations per year for both species, this 

gives a target time goal of 135 generations. For both species, at least 1000 simulations were 

conducted for both the structured refuge and seed-blend options. 

2.3 Results 

 The European corn borer simulations did not show any cases of resistance prior to the 

time goal of 135 generations. The seed-blend option was less durable, but the shortest time 

taken for resistance to evolve in any single simulation was 158 generations. The shortest 

time for the structured refuge option was longer than 1000 generations (see Appendix 3). The 

risk of resistance was somewhat greater with Southwestern corn borer than with European 

corn borer. The seed-blend option was less durable and 4.5% of the simulations resulted in 

resistance evolving in less than 135 generations (ca. 45 years). For the structured refuge 

option, the shortest time to resistance was 118 generations and 0.2% of the simulations 

resulted in the evolution of resistance in less than 135 generations (see Appendix 4). 

Chapter 3 - Discussion 

 The results presented here suggest that there is a low risk of resistance evolving within 

45 years to SmartStax corn in either species (<5%). These may be considered somewhat 

conservative estimates as potential fitness costs were not included in the simulations. Given 

16
the long time frames observed for the evolution of resistance, it is likely that fitness costs 

could have a significant impact in further delaying resistance. The genetics and interactions 

between fitness costs to different toxins are, however, unknown and would add considerably 

to the uncertainty of the results. Given the time frames already observed, the inclusion of 

fitness costs seemed superfluous. 

 The inclusion of three loci significantly increased the complexity of these simulations. 

The non-multiplicative interactions between loci were impossible to accurately partition given 

the information submitted. To the degree possible, inclusion of data on two locus interactions 

in future submissions would help to more accurately evaluate the data. The solutions to the 

problems of three locus simulations presented here are merely a first attempt and will 

hopefully provide a basis for further improvement in methods and data needs to understand 

how resistance evolves in these complex systems. 

References 

Caprio, M. A. and D. Sumerford. 2007. Evaluating transgenic plants for suitability in pest and 

resistance management programs, 2nd edition. L. Lacey & H. Kaya [eds]. Field Manual of 

Techniques for the Application and Evaluation of Entomopathogens. 

Caprio, M.A. and J.A. Glaser 2010a, The Use of Simulation Models to Evaluate Pest 

Resistance Development to Refuge in the Bag Concepts related to Pioneer Submission: 

Letter Report, January 2010, US EPA National Risk Management Research Laboratory, 

Cincinnati, OH. 

Caprio, M.A. and J.A. Glaser 2010b, The Use of Simulation Models to Evaluate Pest 

Resistance Development to Refuge in the Bag Concepts related to Pioneer Submission: 

Letter Report Addendum, April 2010, US EPA National Risk Management Research 

Laboratory, Cincinnati, OH. 

Caprio, M.A. and J.A. Glaser 2010c, Risk Assessment of the Evolution of Resistance in 

Western corn rootworm (Diabrotica virgifera virgifera LeConte) to Smartstax® corn, Letter 

Report, October 2010, US EPA National Risk Management Research Laboratory, Cincinnati, 

OH. 

17
Hibbard, BE, MI Higdon, DP Duran, YM Schweikert, and MR Ellersieck. 2004. Role of egg 

density on establishment and plant-to-plant m7(3): 871-882 

Hibbard, BE, TT Vaughn, IO Oyediran, TL Clark, MR Ellersieck. 2005. Effect of Cry3Bb1 

expressing transgenic corn on plant-to-plant movement by western corn rootworm larvae 

(Coleoptera Chrysomeldiae). J. Econ. Entomol. 98(4): 1126-1138 

Vose, 2000. Risk Analysis: A quantitative guide. Wiley, New York. 

18
 

Variable Min Mode Max Comments

Cry1A.105 0.00001 0.001 0.005

Cry2Ab2 0.00001 0.001 0.01

Cry1F 0.001 0.005 0.01

Overall Survivorship 0.00001 0.001 0.005This will always be as low as lowest of the three loci

Dominance 0.01 0.05 0.2Drawn independently for both loci

Initial frequency 0.001 0.005 0.02

These are variables specific to spatial refuges

Refuge size 0.025 0.04 0.06Always set to 0.05 if seed-mix

Non-compliance 0 0.3 0.6 Always set to 0 for seed-mix

Non-compliance mixing 0.8 0.9 1

Non-random mating 0.2 0.5 0.7Always set to 0 for seed-mix

These are v ariables specific to seed-mixes

Base larval movement 0.3 0.5 0.7Always set to 0 for spatial refuges

Days of larval selection 5 10 15

Asymmetrical larval movement 0.8 0.9 0.95Always set to 0 for spatial refuges

Table 15: PERT parameters for the variables in the structured and larval movement simulations of 

European corn Borer

Variable Min Mode Max Comments

Cry1A.105 0.001 0.05 0.08

Cry2Ab2 0.001 0.05 0.1

Cry1F 0.001 0.005 0.01

Overall Survivorship 0.00001 0.005 0.008This will always be as low as lowest of the three loci

Dominance 0.01 0.05 0.2Drawn independently for both loci

Initial frequency 0.0005 0.001 0.002

These are variables specific to spatial refuges

Refuge size 0.025 0.04 0.06Always set to 0.05 if seed-mix

Non-compliance 0 0.3 0.6 Always set to 0 for seed-mix

Non-compliance mixing 0.8 0.9 1

Non-random mating 0.05 0.3 0.5Always set to 0 for seed-mix

These are v ariables specific to seed-mixes

Base larval movement 0.3 0.5 0.7Always set to 0 for spatial refuges

Days of larval selection 5 10 15

Asymmetrical larval movement 0.8 0.9 0.95Always set to 0 for spatial refuges

Table 26: PERT parameters for the variables in the structured and larval movement simulations of 

Southwestern corn Borer

19
 Survivorship of homozygous 

susceptibles 

Dominance 

Cry1A.105 0.00275 0.11500 

Cry2Ab2 0.00190 0.09976 

Cry1F 0.00502 0.03600 

Overallsurvivorship 0.00140 N/A 

Table 3: The fitnesses and dominance values drawn in this example of fitness calculations 

 Cry1A.105 

 SS RS RR 

Cry2Ab2 SS 0.00140 0.00140 0.00140 

 RS 0.00140 0.00144 0.00176 

RR 0.00140 0.00176 0.00502 

Table 4a: The fitness values calculated for Cry1A.105 and Cry2Ab2 genotypes when the Cry1F 

genotype was SS 

 Cry1A.105 

 SS RS RR 

Cry2Ab2 SS 0.00140 0.00140 0.00142 

 RS 0.00140 0.00182 0.00535 

RR 0.00145 0.00593 0.04084 

Table 4b: The fitness values calculated for Cry1A.105 and Cry2Ab2 genotypes when the Cry1F 

genotype was RS 

 Cry1A.105 

 SS RS RR 

Cry2Ab2 SS 0.00140 0.00146 0.00190 

 RS 0.00153 0.01228 0.10148 

RR 0.00275 0.11764 1.0 

Table 4c: The fitness values calculated for Cry1A.105 and Cry2Ab2 genotypes when the Cry1F 

genotype was RR 

20
Figure 12: Schematic of the deterministic frequency based model for both seed-blend and 

block refuge deployments. 

21
Appendix 1 - Resistance Management Model Used for this Analysis 

The deterministic model is a relatively simple frequency-based model that incorporates noncompliance, non-random mating in compliant fields and larval movement. It includes a switch 

to vary between a block-refuge design and a seed-mixture design. Appendix 2 includes 

descriptions of the different parameters. The primary difference between the Bt population in 

the non-compliant population and the Bt-population in the non-random mating population is 

that no refuge was planted for the non-compliant fields. Non-compliant populations reduce 

both non-random mating and refuge size. (Caprio & Glaser, 2010a&b). 

While this may overestimate movement of refuges by growers, this was incorporated into the 

model by defining the refuge size not as the physical area devoted to the refuge, but rather 

the proportion of the eggs in a field that start out underneath the new refuge. This is the 

effective refuge size to distinguish it from physical refuge size (the proportion of an area 

planted to a refuge). These small refuge sizes are independent of non-compliance and result 

from growers moving their refuge to areas of the field were few eggs were deposited the 

previous year. To estimate this important parameter, we used the stochastic model to estimate 

egg densities at the start of a year across our four different field sizes. In each case, we 

estimated a best and worst scenario. For the best case scenario, the grower moved their 

refuge to the opposite side of the field, next to where the refuge would have been planted in 

the neighboring field the previous year (these simulations were run with all growers planting 

refuges to the same sides of their fields and swapping them every year). The worst case 

scenario was based on planting the refuge in the area of the field where the egg density was 

the lowest. Because the egg density distribution across fields has never been measured, it is 

unlikely that growers could accomplish this, but it does represent and theoretical lower bound 

for effective refuge size. The worst case effective refuge sizes for the four different field sizes 

22
were (from smallest to largest field size): 15%, 11%, 3% and 0.055%. Few insects were noted 

to be moving into the center of large fields as reflected by the effective refuge size of 0.055% 

for 144 ha fields. Other estimates for dispersal are higher and may depend on the region of 

the country. To test the effect of higher dispersal rates we doubled the vector of daily 

movement (to 29m/day) and repeated the large field size simulation In this case the density 

of eggs in the center of the field was fifty times greater under the worst case scenario and the 

effective refuge size was 1%. Most fields in the corn belt, however, are smaller, in the range of 

30-50 ha, and we decided to use a lower limit of 1% for this parameter. 

Asymmetrical movement was included in the model, though the amount movement of 

susceptible larvae from transgenic corn is uncertain (Hibbard et al 2004, Hibbard et al 2005). 

Asymmetrical movement was defined as the probability that larvae that have survived 

exposure on a Bt-plant will move. When larvae move, they move onto plant genotypes 

according to the ratio of those genotypes in the field. In this case, a larva leaving a transgenic 

plant had a 95% probability of moving onto another transgenic plant. As with other rates in the 

model, this movement rate is converted into a daily rate (depending on the days of larval 

movement parameter). 

23
Appendix 2 - Deterministic Parameter Description 

Wss – The fitness of homozygous susceptible genotypes on the transgenic plants. Without 

better knowledge of field tests it's difficult to develop a range, but a specific factor over and 

under was used to deve;op our estimate. One question to be asked might be what is the 

lowest/highest performance under field conditions. Instead of 98.75% mortality, could it be 

95% (as an example)? 

Dominance – The relative fitness of heterozygotes. The traditional population genetics 

equation was employed: 

 Wrs = (Wrr – Wss) * h + Wss 

where h is the dominance and Wrr = 1.0. Our distribution is derived from an estimate of 

h=0.05, but that seems low compared to many estimates based on LC50s. LC50s, however, 

are based on a log scale, however this research has used a linear scale. In general, for 

products with lower effectiveness, dominance tends to be higher ( Caprio et al 2007). 

IRF – The initial resistance allele frequency. My impression from dealing with CRW folks is 

that this value has tended to be higher in empirical estimates than we have seen for other 

insects. 

Block refuge model parameterization: 

Values cited above are shared with the larval movement model and changes in would effect 

both the block and larval movement models to a similar degree. 

24
Refuge proportion – Initially set this to 20%. However, it can be viewed as the proportion of 

eggs in a compliant field that end up being in a refuge (perhaps after being moved), and the 

non-random oviposition patterns can be addressed using this parameter. If more individuals 

survive in refuges, and they tend to lay eggs near that area, then the number of eggs in the 

soil to start the new growing season will be greatest near the refuge from the previous year. If 

the refuge is moved, it stands to reason that the new refuge, even if it covers 20% of the area, 

will receive less than 20% of the eggs oviposited in the field. Based on this modeling 

research, this fact and the likelihood that refuges are moved are identified as key 

assumptions of previous modeling efforts. Unfortunately, there is little empirical data on the 

non-random oviposition patterns of CRW. 

Non-compliance – This is the proportion of fields that are planted without a refuge. 

Individuals emerging from this area mate as a separate population. The most likely value is a 

literature value (CSPI). This parameter should be based on a field by field basis and 

incorporate some estimate of amount of non-compliance. A person planting a 3% refuge in a 

field is not the same as someone not planting any refuge. Alternatively, it shouldn't measure if 

a field is compliant (as yes./no answer), but how compliant a non-compliant field is. General 

comments have been encountered which indicated that 30% non-compliance was too low. It 

is not obvious how non-compliance was estimated. This parameter may be important 

because it determines how much Bt-corn is planted without any refuge. 

Non-compliance mixing – This determines how much movement there is between compliant 

and non-compliant fields. If non-compliant growers are accidentally non-compliant (i.e. they 

never got around to planting the refuge because of weather), then they will generally be 

25
compliant in the next year, and this parameter should be close to 1.0. On the other hand, if 

non-compliant growers are growers who purposefully don't plant refuges and continue this 

practice year-by-year, then this parameter should be low. How low would depend on how 

much movement naturally occurs between non-compliant and compliant fields. Inevitably this 

leads to discussion of the clustering of non-compliant fields. 

Non-random mating – This parameter addresses the proportion of the population in the Btportion of compliant fields that is far enough away from the refuges that they do not mate at 

random. These differ from non-compliant fields in that refuge is planted, but it was too far 

away. The non-random oviposition that results from this non-random mating (and limited 

dispersal capability) is covered by the reductions in refuge size discussed above. 

Larval movement model parameters 

The refuge proportion is fixed at 0.05 in the larval movement model. 

Base larval movement – This is the amount of movement that might occur in a non-transgenic 

field with populations low enough so that significant density-dependent dispersal is not a 

factor. In the stochastic models, population decline is found to be a greater problem than 

density-dependent factors such as mortality or dispersal. The model does include 

asymmetrical movement, with the movement of individuals on Bt plants increased by an 

amount related to their daily mortality rates. Susceptibles therefore move more than 

heterozygotes, and resistant homozygotes would only move at the base larval movement 

rate. Of course, an individual leaving a transgenic plant only has a 5% change of moving 

onto a non-Bt plant. Various estimates of larval dispersal rates have been encountered but 

26
mostly on the low side. Generally, this research assumes that populations will be low in seedmix fields and there will be little to no density-dependent larval dispersal. This could be an 

area of concern since the selection of alternate values is relatively unimportant compared to 

other parameters and a precise estimate would be a misapplication of resources. 

Days of larval movement – Movement and larval development are assumed to occur over this 

period of days. Movement occurs once per day. Daily survivorship is constant and is 

determined as the rate that raised to the days power results in the overall fitness listed for the 

genotype in the original table. In a previous case, that was 0.0125 for susceptible larvae. If 

the days of larval movement is 10, then the daily survivorship rate is 0.645. This means that 

on day 1, 36.5% of the SS larvae on Bt die (assume conservatively that mortality occurs 

before dispersal) , and of the remaining 64.5%, 36.5% will disperse (see above). Of these 

dispersants, only 5% will move onto non-Bt plants, increasing the refuge size from 5% to ca. 

6%. If the base larval movement is low (remember it is spread over days as well), then this 

process can lead to an increase in the refuge population and actually decrease the rate of 

resistance evolution. This feature was inserted in the model to see what effect it might have. 

The value was varied between 5, 10 and 15 and it had very little effect on the results. 

27
Appendix 3 – European Corn Borer Sensitivity Analysis 

Structured Refuge 

> numSummary(ECBStruct2[,"time"], statistics=c("mean", "sd", 

"quantiles"), 

+ quantiles=c(0,.05,.1,.25,.5)) 

 mean sd 0% 5% 10% 25% 50% n 

 48653.23 39796.83 1125 3055.7 4259.8 9187.5 37657 1099 

Call: 

lm(formula = log(time) ~ dominance + IRF + l2dom + l2wss + 

Non.compliance + 

 non.compliance.mixing + non.random.mating + Overall + Refuge + 

 w3dominance + w3ss + wss, data = ECBStruct2) 

Residuals: 

 Min 1Q Median 3Q Max 

-2.42872 -0.58097 0.09332 0.55130 2.42397 

Coefficients: 

 Estimate Std. Error t value Pr(>|t|) 

(Intercept) 12.04053 0.67677 17.791 < 2e-16 *** 

dominance -4.65073 0.71820 -6.476 1.43e-10 *** 

IRF -106.55591 8.46078 -12.594 < 2e-16 *** 

l2dom -7.28007 0.73243 -9.940 < 2e-16 *** 

l2wss -436.50924 17.99658 -24.255 < 2e-16 *** 

Non.compliance -0.69809 0.22527 -3.099 0.00199 ** 

non.compliance.mixing 0.04642 0.69378 0.067 0.94666 

28
non.random.mating 0.31598 0.28462 1.110 0.26717 

Overall 2290.92934 83.97219 27.282 < 2e-16 *** 

Refuge 21.30089 3.90634 5.453 6.13e-08 *** 

w3dominance -0.82344 0.75083 -1.097 0.27301 

w3ss -252.49117 13.89056 -18.177 < 2e-16 *** 

wss -362.84562 31.62867 -11.472 < 2e-16 *** 

--- 

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.7681 on 1086 degrees of freedom 

Multiple R-squared: 0.6285, Adjusted R-squared: 0.6244 

F-statistic: 153.1 on 12 and 1086 DF, p-value: < 2.2e-16 

Stepwise regression 

> stepwise(RegModel.1, direction='backward/forward', criterion='BIC') 

Direction: backward/forward 

Criterion: BIC 

Start: AIC=-501.85 

log(time) ~ dominance + IRF + l2dom + l2wss + Non.compliance + 

 non.compliance.mixing + non.random.mating + Overall + Refuge + 

 w3dominance + w3ss + wss 

29
... 

Step: AIC=-520.22 

log(time) ~ dominance + IRF + l2dom + l2wss + Non.compliance + 

 Overall + Refuge + w3ss + wss 

 Df Sum of Sq RSS AIC 

<none> 642.32 -520.22 

- Non.compliance 1 6.11 648.43 -516.82 

+ non.random.mating 1 0.82 641.49 -514.63 

+ w3dominance 1 0.81 641.51 -514.60 

+ non.compliance.mixing 1 0.01 642.31 -513.23 

- Refuge 1 17.66 659.98 -497.41 

- dominance 1 25.13 667.44 -485.05 

- l2dom 1 58.46 700.77 -431.50 

- wss 1 79.09 721.41 -399.60 

- IRF 1 93.99 736.30 -377.15 

- w3ss 1 195.15 837.47 -235.66 

- l2wss 1 347.86 990.17 -51.58 

- Overall 1 442.41 1084.73 48.66 

Call: 

lm(formula = log(time) ~ dominance + IRF + l2dom + l2wss + 

Non.compliance + Overall + Refuge + w3ss + wss, data = 

ECBStruct2) 

30
Coefficients: 

 (Intercept) dominance IRF l2dom 

l2wss 

 12.186 -4.684 -106.721 -7.289 

-436.892 

Non.compliance Overall Refuge w3ss 

wss 

 -0.723 2296.728 21.362 -252.181 

-365.391 

Standardized regression coefficients 

> lm.beta(RegModel.1) 

 dominance IRF l2dom 

 -0.120300126 -0.234139159 -0.184230039 

 l2wss Non.compliance non.compliance.mixing 

 -0.467551121 -0.057654482 0.001239275 

 non.random.mating Overall Refuge 

 0.020643543 0.571578932 0.101041557 

 

 w3dominance w3ss wss 

 -0.020441152 -0.338343358 -0.232784395 

Bayesian Model Analysis 

31
Call: 

iBMA.bicreg.matrix(x = as.matrix(ECBStruct2[, 1:12]), Y = 

log(ECBStruct2$time)) 

completed 1 iteration 

number of iterations in last call = 100 

variables selected = 1 2 3 4 5 6 7 8 11 

statistics for variables examined so far 

 variable status max Probne0 n times in model 

1 w3ss selected 100 1 

2 Overall selected 100 1 

3 l2wss selected 100 1 

4 IRF selected 100 1 

5 l2dom selected 100 1 

6 dominance selected 100 1 

7 Refuge selected 100 1 

8 Non.compliance selected 86.1 1 

9 w3dominance dropped 4.6 1 

10 non.random.mating dropped 4.7 1 

11 wss selected 100 1 

12 non.compliance.mixing dropped 0 1 

 Results of running BMA on final set of selected variables: 

32
Call: 

bicreg(x = currentX, y = cs$Y, wt = cs$wt, maxCol = cs$maxNvar + 

1) 

 2 models were selected 

 Best 2 models (cumulative posterior probability = 1 ): 

 p!=0 EV SD model 1 model 2 

Intercept 100.0 12.1517 0.2214 12.186 11.962 

w3ss 100.0 -252.3820 13.8802 -252.181 -253.497 

Overall 100.0 2295.0778 83.9948 2296.728 2285.936 

l2wss 100.0 -436.8656 18.0020 -436.892 -436.719 

IRF 100.0 -106.7926 8.4613 -106.721 -107.188 

l2dom 100.0 -7.2734 0.7334 -7.289 -7.188 

dominance 100.0 -4.6883 0.7181 -4.684 -4.714 

Refuge 100.0 21.3890 3.9066 21.362 21.541 

Non.compliance 84.7 -0.6125 0.3323 -0.723 . 

wss 100.0 -364.6580 31.6159 -365.391 -360.598 

 

nVar 9 8 

r2 0.628 0.624 

BIC -1022.529 -1019.104 

post prob 0.847 0.153 

33
Seed-blend 

> numSummary(ECBLM3CBLM[,"time"], statistics=c("mean", "sd", 

"quantiles"), 

+ quantiles=c(0,.05,.1,.25,.5)) 

 mean sd 0% 5% 10% 25% 50% n 

 22687.02 30770.43 158 451.6 638 1303.5 7752 1000 

Call: 

lm(formula = log(time) ~ AsymmMove + BaseLarvalMovement + dominance + 

 IRF + Overall + SelectionDays + w2dominance + w2ss + w3dominance 

+ 

 w3ss + wss, data = ECBLM3CBLM) 

Residuals: 

 Min 1Q Median 3Q Max 

-5.5379 -0.8099 0.1515 0.8006 3.0982 

Coefficients: 

 Estimate Std. Error t value Pr(>|t|) 

(Intercept) 12.85794 1.23390 10.421 < 2e-16 *** 

AsymmMove 0.10445 1.32983 0.079 0.937408 

BaseLarvalMovement -1.65905 0.52120 -3.183 0.001502 ** 

dominance -3.69039 1.05205 -3.508 0.000472 *** 

IRF -137.61374 12.70775 -10.829 < 2e-16 *** 

34
Overall 4142.36059 116.90224 35.434 < 2e-16 *** 

SelectionDays -0.04638 0.01979 -2.344 0.019294 * 

w2dominance -10.82444 1.06005 -10.211 < 2e-16 *** 

w2ss -457.19541 27.30907 -16.742 < 2e-16 *** 

w3dominance -1.99457 1.04612 -1.907 0.056857 . 

w3ss -308.19289 19.43282 -15.859 < 2e-16 *** 

wss -474.23244 40.24046 -11.785 < 2e-16 *** 

--- 

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.059 on 988 degrees of freedom 

Multiple R-squared: 0.6505, Adjusted R-squared: 0.6466 

F-statistic: 167.2 on 11 and 988 DF, p-value: < 2.2e-16 

Stepwise regression: 

> stepwise(RegModel.2, direction='backward/forward', criterion='BIC') 

Direction: backward/forward 

Criterion: BIC 

Start: AIC=185.72 

log(time) ~ AsymmMove + BaseLarvalMovement + dominance + IRF + 

35
 Overall + SelectionDays + w2dominance + w2ss + w3dominance + 

 w3ss + wss 

... 

Step: AIC=173.63 

log(time) ~ BaseLarvalMovement + dominance + IRF + Overall + 

 w2dominance + w2ss + w3ss + wss 

 Df Sum of Sq RSS AIC 

<none> 1117.9 173.63 

+ SelectionDays 1 5.52 1112.4 175.58 

- BaseLarvalMovement 1 10.36 1128.3 175.94 

+ w3dominance 1 3.44 1114.5 177.45 

- dominance 1 13.50 1131.4 178.72 

+ AsymmMove 1 0.00 1117.9 180.53 

- w2dominance 1 115.59 1233.5 265.11 

- IRF 1 131.72 1249.6 278.11 

- wss 1 154.28 1272.2 296.00 

- w3ss 1 283.43 1401.3 392.69 

- w2ss 1 320.52 1438.4 418.81 

- Overall 1 1412.28 2530.2 983.56 

Call: 

lm(formula = log(time) ~ BaseLarvalMovement + dominance + IRF + 

Overall + w2dominance + w2ss + w3ss + wss, data = ECBLM3CBLM) 

36
Coefficients: 

 (Intercept) BaseLarvalMovement dominance 

IRF 

 12.306 -1.579 -3.645 

-137.652 

 Overall w2dominance w2ss 

w3ss 

 4143.697 -10.758 -460.466 

-308.774 

 wss 

 -471.161 

Standardized regression coefficients: 

> lm.beta(RegModel.2) 

 AsymmMove BaseLarvalMovement dominance 

IRF 

 0.001483205 -0.060187208 -0.066242028 -

0.204426528 

 Overall SelectionDays w2dominance 

w2ss 

 0.746663470 -0.044292524 -0.192556319 -

0.334520193 

 w3dominance w3ss wss 

 -0.036173282 -0.298782467 -0.236638839 

Bayesian Model Analysis: 

37
Call: 

iBMA.bicreg.matrix(x = as.matrix(ECBLM3CBLM[, 1:11]), Y = 

log(ECBLM3CBLM$time)) 

completed 1 iteration 

number of iterations in last call = 100 

variables selected = 1 2 3 4 5 6 7 8 9 10 

statistics for variables examined so far 

 variable status max Probne0 n times in model 

1 Overall selected 100 1 

2 w3ss selected 100 1 

3 IRF selected 100 1 

4 w2dominance selected 100 1 

5 w2ss selected 100 1 

6 dominance selected 96.3 1 

7 SelectionDays selected 26.6 1 

8 BaseLarvalMovement selected 79.9 1 

9 wss selected 100 1 

10 w3dominance selected 10.5 1 

11 AsymmMove dropped 0 1 

 Results of running BMA on final set of selected variables: 

38
Call: 

bicreg(x = currentX, y = cs$Y, wt = cs$wt, maxCol = cs$maxNvar + 

1) 

 7 models were selected 

 Best 5 models (cumulative posterior probability = 0.9276 ): 

 p!=0 EV SD model 1 model 2 

Intercept 100.0 12.26600 0.50899 12.30648 

12.74934 

Overall 100.0 4144.30464 117.20047 4143.69684 

4140.57938 

w3ss 100.0 -308.15724 19.50628 -308.77437 -

307.98185 

IRF 100.0 -137.45800 12.74766 -137.65239 -

137.50792 

w2dominance 100.0 -10.75447 1.06398 -10.75756 -

10.80279 

w2ss 100.0 -459.87983 27.36956 -460.46560 -

460.75032 

dominance 96.3 -3.47898 1.24143 -3.64453 -

3.59952 

SelectionDays 26.6 -0.01172 0.02197 . -

0.04378 

BaseLarvalMovement 79.9 -1.26788 0.78852 -1.57893 -

1.58982 

wss 100.0 -472.20715 40.34886 -471.16082 -

474.89833 

39
w3dominance 10.5 -0.19824 0.67056 . . 

 

nVar 8 9 

r2 0.648 

0.649 

BIC -987.47098 -

985.51176 

post prob 0.477 

0.179 

 model 3 model 4 model 5 

Intercept 11.48683 12.46686 11.91835 

Overall 4147.57545 4145.13265 4144.52495 

w3ss -306.74836 -309.03770 -305.95236 

IRF -137.05360 -137.78116 -136.90694 

w2dominance -10.69724 -10.77550 -10.74147 

w2ss -458.31118 -457.23042 -458.57754 

dominance -3.50910 -3.72975 -3.46375 

SelectionDays . . -0.04321 

BaseLarvalMovement . -1.64049 . 

wss -471.46061 -470.24516 -475.15152 

w3dominance . -1.82763 . 

 

nVar 7 9 8 

r2 0.644 0.649 0.646 

BIC -985.14465 -983.63185 -983.02685 

post prob 0.149 0.070 0.052 

40
Appendix 4 – Southwestern Corn Borer Sensitivity Analysis 

Structured refuge 

> numSummary(SWCBStruct2[,"time"], statistics=c("mean", "sd", 

"quantiles"), 

+ quantiles=c(0,.05,.1,.25,.5)) 

 mean sd 0% 5% 10% 25% 50% n 

 8944.133 26217.33 118 207.8 262.9 389 677 1000 

Call: 

lm(formula = log(time) ~ dominance + IRF + l2dom + l2wss + 

Non.compliance + 

 non.compliance.mixing + non.random.mating + Overall + Refuge + 

 w3dominance + w3ss + wss, data = SWCBStruct2) 

Residuals: 

 Min 1Q Median 3Q Max 

-1.5722 -0.6733 -0.2550 0.3588 3.5677 

Coefficients: 

 Estimate Std. Error t value Pr(>|t|) 

(Intercept) 11.7292 0.9155 12.812 < 2e-16 *** 

dominance 0.4475 0.9895 0.452 0.65121 

IRF -299.8843 128.9619 -2.325 0.02025 * 

l2dom -1.1349 0.9816 -1.156 0.24788 

l2wss -35.8689 1.3597 -26.380 < 2e-16 *** 

Non.compliance -1.5360 0.4908 -3.130 0.00180 ** 

41
non.compliance.mixing -3.0610 0.9301 -3.291 0.00103 ** 

non.random.mating 0.2468 0.4251 0.581 0.56168 

Overall 187.2322 18.7328 9.995 < 2e-16 *** 

Refuge 22.6041 5.3726 4.207 2.82e-05 *** 

w3dominance -8.7198 0.9968 -8.748 < 2e-16 *** 

w3ss 111.9971 18.7953 5.959 3.53e-09 *** 

wss -39.2748 1.5177 -25.878 < 2e-16 *** 

--- 

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.9758 on 987 degrees of freedom 

Multiple R-squared: 0.6329, Adjusted R-squared: 0.6285 

F-statistic: 141.8 on 12 and 987 DF, p-value: < 2.2e-16 

Stepwise regression 

> stepwise(RegModel.4, direction='backward/forward', criterion='BIC') 

Direction: backward/forward 

Criterion: BIC 

Start: AIC=27.75 

log(time) ~ dominance + IRF + l2dom + l2wss + Non.compliance + 

 non.compliance.mixing + non.random.mating + Overall + Refuge + 

 w3dominance + w3ss + wss 

... 

42
Step: AIC=7.38 

log(time) ~ l2wss + Non.compliance + non.compliance.mixing + 

 Overall + Refuge + w3dominance + w3ss + wss 

 Df Sum of Sq RSS AIC 

<none> 946.68 7.38 

+ IRF 1 5.11 941.57 8.88 

- Non.compliance 1 9.38 956.06 10.33 

- non.compliance.mixing 1 10.31 956.99 11.30 

+ l2dom 1 1.23 945.45 12.99 

+ non.random.mating 1 0.28 946.40 13.99 

+ dominance 1 0.15 946.53 14.13 

- Refuge 1 17.42 964.10 18.70 

- w3ss 1 33.82 980.51 35.58 

- w3dominance 1 73.28 1019.96 75.03 

- Overall 1 98.33 1045.01 99.29 

- wss 1 632.35 1579.03 512.07 

- l2wss 1 660.31 1607.00 529.63 

Call: 

lm(formula = log(time) ~ l2wss + Non.compliance + 

non.compliance.mixing + Overall + Refuge + w3dominance + w3ss + 

wss, data = SWCBStruct2) 

43
Coefficients: 

 (Intercept) l2wss Non.compliance 

 11.407 -35.632 -1.540 

non.compliance.mixing Overall Refuge 

 -3.057 189.786 22.918 

 w3dominance w3ss wss 

 -8.738 111.878 -38.956 

Standardized regression coefficients 

> lm.beta(RegModel.4) 

 dominance IRF l2dom 

 0.008746106 -0.045226973 -0.022345949 

 l2wss Non.compliance non.compliance.mixing 

 -0.512162511 -0.060586663 -0.063649039 

 non.random.mating Overall Refuge 

 0.011247447 0.201656003 0.081445809 

 

 w3dominance w3ss wss 

 -0.169883903 0.119719428 -0.501846586 

Bayesian Model Analysis 

Call: 

iBMA.bicreg.matrix(x = as.matrix(SWCBStruct2[, 1:12]), Y = 

44
log(SWCBStruct2$time)) 

completed 1 iteration 

number of iterations in last call = 100 

variables selected = 1 2 3 4 5 6 7 8 12 

statistics for variables examined so far 

 variable status max Probne0 n times in model 

1 wss selected 100 1 

2 l2wss selected 100 1 

3 Overall selected 100 1 

4 w3ss selected 100 1 

5 w3dominance selected 100 1 

6 Refuge selected 100 1 

7 non.compliance.mixing selected 90.4 1 

8 Non.compliance selected 84.1 1 

9 l2dom dropped 2.9 1 

10 non.random.mating dropped 0 1 

11 dominance dropped 0 1 

12 IRF selected 30.4 1 

 Results of running BMA on final set of selected variables: 

Call: 

45
bicreg(x = currentX, y = cs$Y, wt = cs$wt, maxCol = cs$maxNvar + 

1) 

 6 models were selected 

 Best 5 models (cumulative posterior probability = 0.9712 ): 

 p!=0 EV SD model 1 model 2 

model 3 

Intercept 100.0 11.208 1.1935 11.407 11.789 

11.136 

wss 100.0 -39.061 1.5229 -38.956 -39.243 

-39.077 

l2wss 100.0 -35.676 1.3609 -35.632 -35.826 

-35.492 

Overall 100.0 189.387 18.7835 189.786 187.077 

191.507 

w3ss 100.0 111.649 18.8479 111.878 112.972 

108.950 

w3dominance 100.0 -8.709 0.9999 -8.738 -8.761 

-8.616 

Refuge 100.0 22.612 5.3836 22.918 22.294 

22.340 

non.compliance.mixing 90.1 -2.766 1.2733 -3.057 -3.089 

-3.073 

Non.compliance 83.6 -1.292 0.7272 -1.540 -1.553 

. 

IRF 31.3 -93.059 155.6262 . -298.641 

. 

 

46
nVar 8 9 

7 

r2 0.630 0.632 

0.627 

BIC -939.693 -938.182 

-936.751 

post prob 0.502 0.236 

0.115 

 model 4 model 5 

Intercept 8.648 11.509 

wss -38.915 -39.360 

l2wss -35.692 -35.682 

Overall 191.424 188.857 

w3ss 110.980 110.001 

w3dominance -8.576 -8.638 

Refuge 22.754 21.721 

non.compliance.mixing . -3.104 

Non.compliance -1.548 . 

IRF . -293.772 

 

nVar 7 8 

r2 0.626 0.629 

BIC -935.760 -935.025 

post prob 0.070 0.049 

47
Seed-blend 

> numSummary(SWCBLM3CBLM[,"time"], statistics=c("mean", "sd", 

"quantiles"), 

+ quantiles=c(0,.05,.1,.25,.5)) 

 mean sd 0% 5% 10% 25% 50% n 

 8171.048 25354.55 73 138.85 158.9 214.75 346 1000 

p<135 generations = 4.5% 

Call: 

lm(formula = log(time) ~ AsymmMove + BaseLarvalMovement + dominance + 

 IRF + Overall + SelectionDays + w2dominance + w2ss + w3dominance 

+ 

 w3ss + wss, data = SWCBLM3CBLM) 

Residuals: 

 Min 1Q Median 3Q Max 

-2.0351 -0.7045 -0.2690 0.3404 3.7590 

Coefficients: 

 Estimate Std. Error t value Pr(>|t|) 

(Intercept) 9.69785 1.31742 7.361 3.83e-13 *** 

AsymmMove -0.87515 1.41242 -0.620 0.535657 

BaseLarvalMovement -1.02658 0.50214 -2.044 0.041176 * 

dominance -3.39043 1.07559 -3.152 0.001670 ** 

48
IRF -20.01400 142.62874 -0.140 0.888434 

Overall 359.61183 20.70324 17.370 < 2e-16 *** 

SelectionDays 0.02491 0.02060 1.209 0.226800 

w2dominance 0.90340 1.05706 0.855 0.392960 

w2ss -32.76668 1.43105 -22.897 < 2e-16 *** 

w3dominance -9.37930 1.04545 -8.972 < 2e-16 *** 

w3ss 74.75041 20.93917 3.570 0.000374 *** 

wss -39.69144 1.70834 -23.234 < 2e-16 *** 

--- 

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.073 on 988 degrees of freedom 

Multiple R-squared: 0.6186, Adjusted R-squared: 0.6144 

F-statistic: 145.7 on 11 and 988 DF, p-value: < 2.2e-16 

Stepwise regression: 

stepwise(RegModel.1, direction='backward/forward', criterion='BIC') 

Direction: backward/forward 

Criterion: BIC 

Start: AIC=211.72 

log(time) ~ AsymmMove + BaseLarvalMovement + dominance + IRF + 

49
 Overall + SelectionDays + w2dominance + w2ss + w3dominance + 

 w3ss + wss 

... 

Step: AIC=183.99 

log(time) ~ dominance + Overall + w2ss + w3dominance + w3ss + 

 wss 

 Df Sum of Sq RSS AIC 

<none> 1145.3 183.99 

+ BaseLarvalMovement 1 4.77 1140.5 186.73 

- dominance 1 11.46 1156.7 187.04 

+ SelectionDays 1 1.67 1143.6 189.44 

+ w2dominance 1 0.79 1144.5 190.21 

- w3ss 1 15.25 1160.5 190.31 

+ AsymmMove 1 0.48 1144.8 190.47 

+ IRF 1 0.03 1145.2 190.87 

- w3dominance 1 97.10 1242.4 258.46 

- Overall 1 346.01 1491.3 441.07 

- w2ss 1 608.00 1753.3 602.93 

- wss 1 624.28 1769.5 612.17 

Call: 

lm(formula = log(time) ~ dominance + Overall + w2ss + w3dominance + 

w3ss + wss, data = SWCBLM3CBLM) 

50
Coefficients: 

(Intercept) dominance Overall w2ss w3dominance 

w3ss 

 8.710 -3.391 358.116 -32.853 -9.551 

76.011 

 wss 

 -39.669 

Standardized regression coefficients 

> lm.beta(RegModel.1) 

 AsymmMove BaseLarvalMovement dominance 

IRF 

 -0.012216973 -0.040440872 -0.062078028 -

0.002765404 

 Overall SelectionDays w2dominance 

w2ss 

 0.351883583 0.023826750 0.016814413 -

0.451905952 

 w3dominance w3ss wss 

 -0.177519991 0.072209652 -0.459232630 

Bayesian Model Analysis 

Call: 

iBMA.bicreg.matrix(x = as.matrix(SWCBLM3CBLM[, 1:11]), Y = 

log(SWCBLM3CBLM$time)) 

51
completed 1 iteration 

number of iterations in last call = 100 

variables selected = 1 2 3 4 5 6 7 

statistics for variables examined so far 

 variable status max Probne0 n times in model 

1 w2ss selected 100 1 

2 wss selected 100 1 

3 Overall selected 100 1 

4 w3dominance selected 100 1 

5 w3ss selected 100 1 

6 dominance selected 83 1 

7 BaseLarvalMovement selected 19.3 1 

8 w2dominance dropped 0 1 

9 IRF dropped 0 1 

10 AsymmMove dropped 0 1 

11 SelectionDays dropped 4.1 1 

 Results of running BMA on final set of selected variables: 

Call: 

bicreg(x = currentX, y = cs$Y, wt = cs$wt, maxCol = cs$maxNvar + 

1) 

52
 4 models were selected 

 Best 4 models (cumulative posterior probability = 1 ): 

 p!=0 EV SD model 1 model 2 

model 3 

Intercept 100.0 8.7703 0.2905 8.710 9.202 

8.491 

w2ss 100.0 -32.8635 1.4329 -32.853 -32.780 -

32.995 

wss 100.0 -39.6477 1.7068 -39.669 -39.546 -

39.689 

Overall 100.0 358.1576 20.7353 358.116 360.542 

355.602 

w3dominance 100.0 -9.5344 1.0452 -9.551 -9.375 -

9.657 

w3ss 100.0 75.7921 20.9239 76.011 74.450 

76.545 

dominance 82.3 -2.7928 1.6221 -3.391 -3.405 

. 

BaseLarvalMovement 20.1 -0.2052 0.4667 . -1.021 

. 

 

nVar 6 7 

5 

r2 0.616 0.618 

0.612 

BIC -915.692 -912.960 -

912.649 

post prob 0.656 0.167 

53
0.143 

 model 4 

Intercept 8.977 

w2ss -32.923 

wss -39.568 

Overall 357.995 

w3dominance -9.484 

w3ss 75.002 

dominance . 

BaseLarvalMovement -1.012 

 

nVar 6 

r2 0.614 

BIC -909.772 

post prob 0.034 

54