Overview of Major Approaches

Before going into a detailed cost-benefit analysis of embryo selection, I’ll give a rough overview of the various developing approaches for genetic engineering of complex traits in humans, compare them, and briefly discuss possible timelines and outcomes. (References/analyses/code for particular claims are generally provided in the rest of the text, or in some cases, buried in my embryo editing rough notes, and omitted here for clarity.)

The past 2 decades have seen a revolution in molecular genetics: the sequencing of the human genome kicked off an exponential reduction in genetic sequencing costs which have dropped the cost of genome sequencing from millions of dollars to $25.01^$20₂₀₁₆ (SNP genotyping)–$625.32^$500₂₀₁₆ (whole genomes). This has enabled the accumulation of datasets of millions of individuals’ genomes which allow a range of genetic analyses to be conducted, ranging from SNP heritabilities to detection of recent evolution to GWASes of traits to estimation of the genetic overlap of traits.

The simple summary of the results to date is: behavioral genetics was right. Almost all human traits, simple and complex, are caused by a joint combination of environment, stochastic & randomness, and genes. These patterns can be studied by methods such as family, twin, adoption, or sibling studies, but ideally are studied directly by reading the genes of hundreds of thousands of unrelated people, which yield estimates of the effects of specific genes and predictions of phenotype values from entire genomes. Across all traits examined, genes cause ~50% of differences between people in the same environment, factors like randomness & measurement-error explain much of the rest, and whatever is left over is the effect of nurture. Evolution is true, and genes are discrete physical patterns encoded in chromosomes which can be read and edited, with simple traits such as many diseases being determined by a handful of genes, yielding complicated but discrete behavior, while complex traits are instead governed by hundreds or thousands of genes whose effects sum together and produce a normal distribution such as IQ or risk of developing a complicated disease like schizophrenia. This allows direct estimation of their genetic contribution to a phenotype, as well as that of their children. These genetic traits contribute to many observed societal patterns, such as the children of the rich also being richer and smarter and healthier, why poorer neighborhoods have sicker people, relatives of schizophrenics are less intelligent, etc; these traits are substantially heritable, and traits are also interconnected in an intricate web of correlations where one trait causes another and both are caused by the same genetic variants. For example, intelligence-related variants are uniformly inversely correlated with disease-related variants, and positively correlated with desirable traits. These results have been validated by many different approaches and the existence of widespread large heritabilities, genetic correlations, and valid PGSes are now academic consensus.

Because of this pervasive genetic influence on outcomes, genetic engineering is one of the great open questions in transhumanism: how much is possible, with what, and when?

Suggested interventions can be broken down into a few categories:

• cloning (copying)

• selection (variation with ranking)

• editing (rewriting)

• synthesis (writing)

Each of these has different potentials, costs, and advantages & disadvantages:

Overall I would summarize the state of the field as:

• cloning: is unlikely to be used at any scale for the foreseeable future despite its power, and so can be ignored (except inasmuch as it might be useful in another technology like IES or genome synthesis)

• simple single-trait embryo selection: is strictly inferior to simple multiple embryo selection, and there is no reason to use it other than the desire to save a tiny bit of statistical effort, and much reason to use it (larger and safer gains), so it need not be discussed except as a strawman.

• simple multiple-trait embryo selection: available & profitable now (first baby was born mid-2020), but is too limited in possible gains, requires a far too onerous process (IVF) for more than a small percentage of the population to use it, and is more or less trivial in consequences.

  As median embryo count in IVF hovers around 5, the total gain from selection is small, and much of the gain is wasted by losses in the IVF process (the best embryo doesn’t survive storage, the second-best fails to implant, and so on). One of the key problems is that polygenic scores are the sum of many individual small genes’ effects and form a normal distribution, which is tightly clustered around a mean. A polygenic score is attempting to predict the net effect of thousands of genes which almost all cancel out, so even accurate identification of many relevant genes still yields an apparently unimpressive predictive power. The fact that traits are normally distributed also creates difficulties for selection: the further into the tail one wants to go, the larger the sample required to reach the next step—to put it another way, if you have 10 samples, it’s easy (a 1 in 10 probability) that your next random sample will be the largest sample yet, but if you have 100 samples, now the probability of an improvement is the much harder 1 in 100, and if you have 1,000, it’s only 1 in 1,000; and worse, if you luck out and there’s an improvement, the improvement is ever tinier. After taking into account existing PGSes, previously reported IVF process losses, costs, and so on, the implication that it is moderately profitable and can increase traits perhaps 0.1SD, rising somewhat over the next decade as PGSes continue to improve, but never exceeding, say, 0.5SD.

  Embryo selection could have substantial societal impacts in the long run, especially over multiple generations, but this would both require IVF to become more common and for no other technology to supersede it (as they certainly shall). When IVF began, many pundits proclaimed it would “forever change what it means to be human” and other similar fatuosities; it did no such thing, and has since productively helped countless parents & children, and I fully expect embryo selection to go the same way. I would consider embryo selection to have been considerably overhyped (by those hyperventilating about “Gattaca being around the corner”), and, ironically, also underhyped (by those making arguments like “trait X is so polygenic, therefore embryo selection can’t work”, which is statistically illiterate, or “traits are complex interactions between genes and environment most of which we will never understand”, which is obfuscating irrelevancy and FUD).

  Embryo selection does have the advantage of being the easiest to analyze & discuss, and the most immediately relevant.

• massive multiple embryo selection: the single most binding constraint on simple embryo selection (single or multiple trait), is the number of embryos to work with, which, since paternal sperm is effectively infinite, means number of eggs.

  For selection, the key question is what is the most extreme or maximum item in the sample; a small sample will not spread wide, but a large sample will have a bigger extreme. The more lottery tickets you buy, the better the chance of getting 1 ticket which wins a lot. Whereas, the PGS, to peoples’ general surprise, doesn’t make all that much of a difference after a little while. If you have 3 embryos, even going from a noisy to a perfect predictor, it doesn’t make much of a difference, because no matter how flawless your prediction, embryo #1 (whichever it is) out of 3 just isn’t going to be all that much better than average; if you have 300 embryos, then a perfect predictor becomes more useful.

  There is no foreseeable way to safely extract more eggs from a donor: standard IVF cycle approaches appear to have largely reached their limit, and stimulating more eggs’ release in a harvesting cycle is dangerous. A different approach is required, and it seems the only option may be to make more eggs. One possibility is to not try to stimulate release of a few eggs and collect them, but instead biopsy samples of proto-eggs and then hurry them in vitro to maturity as full eggs, and get many eggs that way; biopsies might be compelling without selection at all: the painful, protracted, failure-prone, and expensive egg harvesting process to get ~5 embryos, which then might yield a failed cycle anyway, could be replaced by a single quick biopsy under anesthesia yielding hundreds of embryos effectively ensuring a successful cycle. Less invasively, laboratory results in inducing regression to stem cell states and then oogenesis have made steady progress over the past decades in primarily rat/mice but also human cells, and researchers have begun to speak of the possibility in another 5 or 10 years of enabling infertile or homosexual couples to conceive fully genetically-related children through somatic ↔︎ gametic cell conversions. This would also likely allow generating scores or hundreds of embryos by turning easier-to-acquire cells like skin cells or extracted eggs into stem cells which can replicate and then be converted into egg cells & fertilized. While it is still fighting the normal distribution with brute force, having 500 embryos works a lot better than having just 5 embryos to choose from. The downside is that one still needs to biopsy and sequence each embryo in order to compute their particular PGS; since one is still fighting the thin tail, at some point the cost of creating & testing another embryo exceeds the expected gain (probably somewhere in the hundreds of embryos).

  Unlike simple embryo selection, this could yield immediately important gains like +2SD. IVF yield ceases to be much of a problem (the second/third/fourth-best embryos are now almost exactly as good as the first-best was and they probably won’t all fail), and enough brute force has been applied to reach potentially 1-2SD in practice. If taken up by only the current IVF users and applied to intelligence alone, it would immediately lead to the next generation’s elite positions being dominated by their kids; if taken up by more and done properly on multiple traits, the advantage would be greater.

• Gamete selection/Optimal Chromosome Selection: only a theoretical possibility at the moment, as there is no direct way to sequence individual sperm/eggs or manipulate chromosome choice. GS/OCS are interesting more for the points they make about variance & order statistics & the CLT: it results in a much larger gain than one would expect simply by switching perspectives and focusing on how to select earlier in the ‘pipeline’, so to speak, where variance is greater because sets of genes haven’t yet been combined in one package & canceled each other out. If someone did something clever to allow inference on gametes’ PGSes or select individual chromosomes, then it could yield an immediate discontinuously large boost in trait-value of +2SD in conjunction with whatever embryo selection is available at that point.

• Iterated Embryo Selection: If IES were to happen, it would allow for almost arbitrarily large increases in trait-values across the board in a short period of time, perhaps a year. While IES has major disadvantages (extremely costly to produce the first optimized embryos, depending on how many generations of selection are involved; selection has some inherent speed limits trading off between accidentally losing possibly useful variants & getting as large a gain each generation as possible; embryos are unlikely to resemble the original donors at all without an additional generation ‘backcrossed’ with the original donor cells, undoing most of the work), the extreme increases may justify use of IES and create demand from parents. This could then start a tsunami. Depending on how far IES is pushed, the first release of IES-optimized embryos may become one of the most important events in human history.

  IES is still distant and depends on a large number of wet lab breakthroughs and finetuned human-cell protocols. Coaxing scores or hundred of cells through all the stages of development and fertilization, for multiple generations, is no easy task. When will IES be possible? The relevant literature is highly technical and only an expert can make sense of it, and one should have hands-on expertise to even try to make forecasts. There are no clear cost curves or laws governing progress in stem cell/gamete research which can be used to extrapolate. Perhaps no one will ever put all the money and consistent research effort into developing it into something which could be used clinically. Just because something is theoretically possible and has lots of lab prototypes doesn’t mean that the transition will happen. (Look at human cloning; everyone assumed it’d happen long ago, but as far as anyone knows, it never has.) On the other hand, perhaps someone will.

  IES is one of the scariest possibilities on the list, and the hardest to evaluate; it seems clear, at least, that it will certainly not happen in the next decade, but after that…? IES has been badly under-discussed to date.

• Gene Editing: the development of CRISPR has led to more hype than embryo selection itself. However, the current family of CRISPR techniques & previous alternatives & future improvements, can be largely dismissed on statistical grounds alone. Even if we hypothesized some super-CRISPR which could make a handful of arbitrary SNP edits with zero risk of mutation or other forms of harm, it would not be especially useful and would struggle to be competitive with embryo selection, let alone IES/OCS/genome synthesis. The unfixable root cause is the polygenicity of the most important polygenic traits (which is a blessing for selection or synthesis approaches, as it creates a vast reservoir of potential improvements, but a curse for editing), and to a lesser extent, the asymmetry of effect sizes (harmful variants are more harmful than beneficial ones are beneficial).

  The benefit of gene editing a SNP is the number of edits times the SNP effect of each edit times the probability the effect is causal. Probability it’s causal? Can’t we assume that the top hits from large GWASes these days have a posterior probability ~100% of having a non-zero effect? No. This is because of a technical detail which is largely irrelevant to selection processes but is vitally important to editing: the hits identified in a PGS are not necessarily the exact causal base-pair(s). Often they are, but more often they are not. They are instead proxies for a neighboring causal variant which happens to usually be inherited with it, as genomes are inherited in a chunky fashion, in big blocks, and do not split & recombine at every single base-pair. This is no problem for selection—it predicts great and is cheaper & easier to find a correlated SNP than the true causal variant. But it is fatal to editing: if you edit a proxy, it’ll do nothing (or maybe it’ll do the opposite).

  How fatal is this? Attempts at “fine-mapping” or using large datasets to distinguish which of a few SNPs is the real culprit or seeing how PGSes’ performance shrinks when going from the original GWAS population to a deeply genetically different population like Subsaharan Africans who have totally different proxy patterns (if there is non-zero prediction power, it must be thanks to the causal hits, which act the same way in both populations), we can estimate that the causal probability may be as low as 10%. Combine this with the few edits safely permitted, perhaps 5, the small effect size of each genetic variant, like 0.2 IQ points for intelligence, and the effect becomes dismal. A tenth of a point? Not much. Even if we had all causal variants, the small average effect size, combined with few possible edits, is no good. Fix the causal variant problem, and it’s still only 5 edits at 0.2 points each. Nor is IQ at all unique in this respect—it’s somewhat unusually polygenic, but a cleaner trait like height still implies small gains such as half an inch.

  What about rare variants? The problem with rare variants is that they are rare, and also not of especially large beneficial effect. Being rare makes them hard to find in the first place, and the lack of benefit (as compared to a baseline human without said variant) means that they are not useful for editing. We might find many variants which damage a trait by a large amount, say, increasing abdominal fat mass by a kilogram or lowering IQ by a dozen points, but of course, we don’t want to edit those in! (They also aren’t that important for any embryo selection method, because they are rare, not usually present, and thus there is usually no selection to be done.) We could hope to find some variant which increases IQ by several points—but none have been found, if they were at all common they would’ve been found a long time ago, and indirect methods like DeFries-Fulker regression suggest that there are few or no such rare variants. Nor is measuring other traits a panacea: if there were some variant which increased IQ by a medium amount by increasing a specific trait like working memory which has not been studied in large GWASes or DeFries-Fulker regressions to date, then such a WM-boosting variant should’ve been detected through its mediated effect, and to the extent that it has no effect on hard endpoints like IQ or education or income, it then must be questioned how useful it is in the first place. The situation may be somewhat better with other traits (there’s still hope for finding large beneficial effects¹, and in the other direction, disease traits tend to have more rare variants of larger effects which might be worth fixing in relatively many individual cases, like BRCA or APOE) but I just don’t see any realistic way to reach gains like +1SD on anything with gene editing methods in the foreseeable future using existing variants.

  What about non-existing variants ie brand-new variants based on extrapolation from human genetic history or animal models? These hypothetical mutations/edits could have large effects even if we have failed to find any in the wild. But the track record of animal models in predicting complex human systems such as the brain is not good at all, and such large novel mutations would have zero safety record, and how would you prove any were safe without dozens of live births and long-term followup—which would never be permitted? Given the poor prior probability of both safety & efficacy, such mutations would simply remain untried indefinitely.

  It is difficult to see how to remedy this in any useful way. The causal probability will creep up as datasets expand & cross-racial GWASes become more common, but that doesn’t resolve the issue after we increase the gain by a factor of 10. The limit is still the edit count: the unique edit limit of ~5 is not enough to work with. Can this be combined usefully with IES to do edits per generation? Likely but you still need IES first! Can the edit limit be lifted? …Maybe. Genetic editing follows no predictable improvement curve, or learning curve, and doesn’t benefit directly from any exponentials. It is hard to forecast what improvements may happen. 2019 saw a breakthrough from a repeated-edit SOTA of ~60 edits in a cell to ~2,600 (Smith et al 2019), which no one forecast, but it’s unclear when if ever that would transfer to useful per-SNP edits; but nevertheless, the possibility of mass editing cannot be ruled out.

  So, CRISPR-style editing may be revolutionary in rare genetic diseases, agriculture, & research, but as far as we are concerned, it has been grossly overhyped: there is a chance it will live up to the most extreme claims, but not a large one.

• Genome synthesis: the simple answer to gene editing’s failure is to observe that if you have to make possibly thousands of edits to fix up a genome to the level you want it, why not go out and make your own genome? (with blackjack and hookers…) That is the audacious proposal of genome synthesis. It sounds crazy, since genome synthesis has historically been mostly used to make short segments for research, or perhaps the odd pandemic virus, but unnoticed by most, the cost per base-pair has been crashing for decades, allowing the creation of entire yeast genomes and leading to the recent HGP-Write proposal from George Church & others to invest in genome synthesis research with the aim of inventing methods which can create custom genomes at reasonable prices. Such an ability would be staggeringly useful: custom organisms designed to produce arbitrary substances, genomes with the fundamental encoding all swapped around rendering them immune to all viruses ever, organisms with a single giant genome or with all mutations replaced with the modal gene, among other crazy things. One could also, incidentally, use cheap genome synthesis for bulk storage of data in a dense, durable, room-temperature format (explaining both Microsoft & IARPA’s interest in funding genome synthesis research).

  Of course, if you can synthesize an entire genome—a single chromosome would be almost as good to some extent—you can take a baseline genome and make as many ‘edits’ as you please. Set all the top variants for all the relevant traits to the estimated best setting. The possible gains are greater than IES (since you are not limited by the initial gene pool of starting variants nor by the selection process itself), and one can increase traits by hundreds of SDs (whatever that means).

  

  

  Genome synthesis, unlike IES, has historically proceeded on a smooth cost-curve, has many possible implementations, and has many research groups & startups involved due to its commercial applications. A large-scale “HGP-Write” has been proposed to scale genome synthesis up to yeast sized organisms and eventually human-sized genomes. The cost curve suggests that around 2035, whole human genomes reach well-resourced research project ranges of $10-30m; some individuals in genome synthesis tell me they are optimistic that new methods can greatly accelerate the cost-curve. (Unlike IES, genome synthesis is not committed to a particular workflow, but can use any method which yields, in the end, the desired genome; all of these methods can be researched in parallel & only one needs to work, representing a major advantage.) Genome synthesis has many challenges before one could realistically implant an embryo, such as ensuring all the relevant structural features like methylation are correct (which may not have been necessary for earlier more primitive/robust organisms like yeast), and so on, but whatever the challenges for genome synthesis, the ones for IES appear greater. It is entirely possible that IES will develop too slowly and will be obsoleted by genome synthesis in 10-20 years. The consequences of genome synthesis would be, if anything, larger than IES because the synthesis technology will be distributed in bulk, will probably continue decreasing in cost due to the commercial applications regardless of human use, and don’t require rare specialized wet lab expertise but like genome sequencing, will almost certainly become highly automated & ‘push button’.

  If IES has been under-discussed and is underrated, genome synthesis has not been discussed at all & vastly more underrated.

To sum up the timeline: CRISPR & cloning are already available but will remain unimportant indefinitely for various fundamental reasons; multiple embryo selection is useful now but will always be minor; massive multiple embryo selection is some ways off but increasingly inevitable and the gains are large enough on both individual & societal levels to result in a shock; IES will come sometime after massive multiple embryo selection but it’s impossible to say when, although the consequences are potentially global; genome synthesis is a similar level of seriousness, but is much more predictable and can be looked for, very loosely, 2030^–₁₀2040 (and possibly sooner).

FAQ: Frequently Asked Questions

Readers already familiar with the idea of embryo selection may have some common misconceptions which would be good to address up front:

1. IVF Costs: IVF is expensive, somewhat dangerous, and may have worse health outcomes than natural childbirth

   I agree, but we can consider the case where these issues are irrelevant. It is unclear what the long-run effects of IVF on children may be, other than the harm probably isn’t too great; the literature on IVF suggests that the harms are probably very small and smaller than, for example, paternal age effects, but it’s hard to be sure given that IVF usage is hardly exogenous and good comparison groups for even just correlational analysis are hard to come by. (Natural-born children are clearly not comparable, but neither are natural-born siblings of IVF children—why was their mother able to have one child naturally but needed IVF for the next?) I would not recommend anyone do IVF solely to benefit from embryo selection (as opposed to doing PGD to avoid passing a horrible genetic disease like Huntington’s, where it is impossible for the hypothetical harms of IVF to outweigh the very real harm of that genetic disease). Here I consider the case where parents are already doing IVF, for whatever reason, and so the potential harms are a “sunk cost”: they will happen regardless of the choice to do embryo selection, and can be ignored. This restricts any results to that small subset (~1% of parents in the USA as of 2016), of course, but that subset is the most relevant one at present, is going to grow over time, and could still have important societal effects.

   An interesting question would be, at what point does embryo selection become so compelling that would-be parents with a family history of disease (such as schizophrenia) would want to do it? (Because of the nonlinear nature of liability-threshold polygenic traits and relatively rare diseases like schizophrenia, someone with a family history benefits far more than someone with an average risk; see the truncation selection/multiple-trait selection on why this implies that selection against diseases is not as useful as it seems.) What about would-be parents with no particular history? How good does embryo selection need to be for would-be parents who could conceive naturally to be willing to undergo the cost (~$10k even at the cheapest fertility clinics) and health risks (for both mother & child) to benefit from embryo selection? I don’t know, but I suspect “simple embryo selection” is too weak and it will require “massive embryo selection” (see the overview for definitions & comparisons).

2. PGSes Don’t Work: GWASes merely produce false positives and can’t do anything useful for embryo selection because they are false positives/population structure/publication bias/etc…

   Some readers overgeneralize the debacle of the candidate-gene literature, which is almost 100% false-positive garbage, to GWASes; but GWASes were designed in response to the failure of candidate-genes by much more stringent thresholds & large datasets & more population structure correction, and have performed well as datasets reached necessary sizes. Their PGSes predict out-of-sample increasingly large amounts of variance, the PGSes have high genetic correlations between cohorts/countries/times/measurement methods, and they work within-family between siblings, who by definition have identical ancestries/family backgrounds/SES/etc but have randomized inheritance from their parents. For a more detailed discussion, see the section, “Why Trust GWASes?”. (While GWASes are indeed highly flawed, those flaws typically work in the direction of inefficiency/reducing their predictive power, not inflating them.)

3. The Prediction Is Noncausal: GWASes may be predictive but this is irrelevant because the SNPs in a PGS are merely non-causal variants which proxy for causal variants

   Background: in a GWAS, the measured SNPs may cause the outcome or they may merely be located on a genome nearby a genetic variant which has the causal effect; because genomes are inherited in a ‘chunky’ fashion, a measured SNP may almost always be found alongside the causal genetic variant within a particular population. (Over a long enough timeframe, as organisms reproduce, that part of the genome will be broken up, but this may take centuries or millennia.) Such a SNP is in “linkage disequilibrium” or just LD. Such a scenario is quite common, and may in fact be the case for the overwhelming majority of SNPs in human GWASes. This is both a blessing and a curse for GWASes: it means that easy cheaply-measured SNPs can probe harder-to-find genetic variants, but it also means that the SNPs are not causal themselves. So for example, if one took a list of SNPs from a GWAS, and used CRISPR to edit them, most of the edits would do nothing. This is a serious concern for genetic engineering approaches—just because you have a successful GWAS doesn’t mean you know what to edit!

   But is this a problem for embryo selection? No. Because you are not engaged in any editing or causal manipulation. You are passively observing and predicting what is the best embryo in a sample. This does not disturb the LD patterns or break any correlations, and the predictions remain valid. Selection doesn’t care what the causal variants are, it cares only that, whatever they are or wherever they are on the genome, the chosen embryo has more of them than the not-chosen embryos. Any proxy will do, as long as it predicts well. In the long run, changes in LD will gradually reduce the PGS’s predictive power as the SNPs become better/worse proxies, but this is unimportant since there will be many GWASes in between now and then, and one would be upgrading PGSes for other reasons (like their steadily improving predictive power regardless of LD patterns).

4. PGSes Predict Too Little: Embryo selection can’t be useful with PGSes predicting only X% [where X% > state-of-the-art] of individual variance

   The mistake here is confusing a statistical measure of error with the goal. Any default summary statistic like R² or RMSE is merely a crutch with tenuous connections to optimal decisions. In embryo selection, the goal is to choose better embryos than average to implant rather than implant random embryos, to get a gain which pays for the costs involved. A PGS only needs to be accurate enough to select a better embryo out of a (typically small) batch. It doesn’t need to be able to predict future, say, IQ, within a point. Estimating the precise future trait value of an embryo may be quite difficult, but it’s much easier to predict which of two embryos will have a higher trait value. (It’s the difference between predicting the winner of a soccer game and predicting the exact final score; the latter does let one do the former, but the former is what one needs and is much easier.) Once your PGS is good enough to pick the best or near-best embryo, even a far better PGS makes little difference—after all, one can’t do any better than picking the best embryo out of a batch. And due to diminishing returns/tail effects, the larger the batch, the smaller the difference between the best and the 4th-best etc, reducing the regret. (In a batch of 2, there’s little difference between a poor and a perfect predictor; and in a batch of 1, there’s none.)

   To decide whether a PGS of X% is adequate cannot be done in a vacuum; the necessary performance will depend critically on the value of a trait, the cost of embryo selection, the losses in the IVF pipeline, and most importantly of all, the number of embryos in each batch. (The final gain depends the most on the embryo count—a fact lost on most people discussing this topic.) As embryo selection is cheap at the margin, and ranking is easier than regression, this can be done with surprisingly poor PGSes, and the bar of profitability is easy to meet, and for embryo selection, has been met for some years now (see the rest of this page for an analysis of the specific case of IQ).

   • The genome-wide statistically-significant hits explain <X% of individual variance:

     Statistical-significance thresholds are essentially arbitrary. There is no need to fetishize them: they do not correspond to any posterior probability of a hit being “real”, introduce many serious difficulties of interpretation due to power (if a GWAS has a hit on an SNP with an estimated effect size of X, and a second GWAS also estimates it at X but due to a slightly higher standard error, it is no longer “statistically-significant”, what does that mean, exactly?) and even if they did, the number of false positives has little relationship to the predictive power, much less selection gain of a PGS, much less the final profit of embryo selection. The relevant question is what are the best predictions which can be made? For human complex traits, the most accurate predictions typically use a PGS based on most of or all measured variants. Anything less is less.

5. Unintended Consequences: Selection on traits, especially intelligence, will backfire horribly

   It is hypothetically possible for selection on one trait, which happens to be inversely correlated on a genetic level, with another important trait, to backfire by increasing the first trait but then doing much more damage by decreasing the second trait. This occurs occasionally in long-term or intense breeding programs, and has been demonstrated by very carefully-designed group-selection experiments such as the famous chicken-crate experiment.

   However, for humans, such genetic correlations are highly unlikely a priori as we can simply observe broad patterns like the global correlations of SES/wealth/intelligence/health with all desirable outcomes (“Cheverud’s conjecture”), and countless human genetic correlations have already been calculated by various methods and are now routinely reported in GWASes, and invariably diseases positively correlate with diseases and good things correlate with other good things. Whatever harmful backfire effects there may be are far outweighed by the beneficial backfire effects, so selection on a single trait, especially intelligence, is not going to incur these speculative hypothetical harms.

   If there are any such harms, they can be reduced or eliminated by simply taking into account multiple traits while selecting, and doing multi-trait selection. This is easy to do with the present availability of PGSes on hundreds of traits—given that all the hard work is in the genotyping step, why would one ignore all traits but one and throw away all that data? In fact, even if there were no possibility of backfire effects, embryo selection would be done with multi-trait selection anyway, simply because it is so easy and the benefits are so compelling: using multiple traits allows for much greater overall gains because two embryos similar or identical on one trait may differ a great deal on another trait, and when traits are genetically correlated, they can serve as proxies for each other, producing effective boosts in predictive power. For all these reasons, most breeding programs use multi-trait selection. For more details and an example of the benefits in embryo selection, see the multiple-selection section.

Embryo Selection Cost-Effectiveness

In vitro fertilization (IVF) is a medical procedure for infertile women in which eggs are extracted, fertilized with sperm, allowed to develop into an embryo, and the embryo injected into their womb to induce pregnancy. The choice of embryo to implant is usually arbitrary, with some simple screening for gross abnormalities like missing chromosomes or other cellular defects, which would either be fatal to the embryo’s development (so useless & wasteful to implant) or cause birth defects like Down syndrome (so much preferable to implant a healthier embryo).

However, various tests can be run on embryos, including genome sequencing after extracting a few cells from the embryo, which is called: Preimplantation genetic profiling / preimplantation genetic diagnosis / preimplantation genetic screening (PGD; review)—when genetic information is measured and used to choose which embryo to implant. PGD has historically been used primarily to detect and select against a few rare recessive genetic diseases with single-gene causes like the fatal Huntington’s disease: if both parents are carriers, an embryo without the recessive can be chosen, or at least, an embryo which is heterozygous and won’t develop the disease. This is useful for those unlucky enough to have a family history or be known carriers, but while initially controversial, is now merely an obscure & useful part of fertility medicine.

However, with ever cheaper SNP arrays and the advent of large GWASes in the 2010s, large amounts of subtler genetic information becomes available, and one could check for abnormalities and also start making useful predictions about adult phenotypes: one could choose embryos with higher/lower probability of traits with many known genetic hits such as height or intelligence or alcoholism or schizophrenia—thus, in effect, creating “designer babies” with proven technology no more exotic than IVF and 23andMe. Since such a practice is different in so many ways from traditional PGD, I’ll call it “embryo selection”.

Embryo selection has already begun to be used by the most sophisticated cattle breeding programs (Mullaart & Wells2018) as an adjunct to their highly successful genomic selection & embryo transfer programs, and use in humans poses no notable challenges. The first human baby whose IVF selection process involved polygenic scores², not just standard monogenic PGD selection, was born in mid-2020.

What traits might one want to select on? For example, increases in height have long been linked to increased career success & life satisfaction with estimates like +$800 per inch per year income, which combined with polygenic scores predicting a decent fraction of variance, could be valuable³ But height, or hair color, or other traits are in general zero-sum traits, often easily modified (eg. hair dye or contact lenses), and far less important to life outcomes than personality or intelligence, which profoundly influence an enormous range of outcomes ranging from academic success to income to longevity to violence to happiness to altruism (and so increases in which are far from “frivolous”, as some commenters have labeled them); since the personality GWASes have had difficulties (probably due to non-additivity of the relevant genes connected to predicted frequency-dependent selection, see Penke et al 2007/Penke & Jokela2016), that leaves intelligence as the most important case.

Discussions of this possibility have often led to both overheated prophecies of “genius babies” or “super-babies”, and to dismissive scoffing that such methods are either impossible or of trivial value; unfortunately, specific numbers and calculations backing up either view tend to be lacking, even in cases where the effect can be predicted easily from behavioral genetics and shown to be not as large as laymen might expect & consistent with the results (for example, the “genius sperm bank”⁴).

In “Embryo Selection for Cognitive Enhancement: Curiosity or Game-changer?”, Shulman & Bostrom2014 consider the potential of embryo selection for greater intelligence in a little detail, ultimately concluding that in the most applicable current scenario of minimal uptake (restricted largely to those forced into IVF use) and gains of a few IQ points, embryo selection is more of “curiosity” than “game-changer” as it will be “Socially negligible over one generation. Effects of social controversy more important than direct impacts.” Some things are left out of their analysis which I’m interested in:

1. they give the upper bound on the IQ gain that can be expected from a given level of selection & then-current imprecise GCTA heritability estimates, but not the gain that could be expected with updated figures: is it a large or small fraction of that maximum? And they give a general description of what societal effects might be expected from combinations of IQ gains and prevalence, but can we say something more rigorously about that?

2. their level of selection may bear little resemblance to what can be practically obtained given the realities of IVF and high embryo attrition rates (selecting from 1 in 10 embryos may yield x IQ points, but how many real embryos would we need to implement that, since if we extract 10 embryos, 3 might be abnormal, the best candidate might fail to implant, the second-best might result in a miscarriage, etc?)

3. there is no attempt to estimate costs nor whether embryo selection right now is worth the costs, or how much better our selection ability would need to be to make it worthwhile. Are the advantages compelling enough that ordinary parents, who are already using IVF and could use embryo selection at minimal marginal cost, would pay for it and take the practice out of the lab? Under what assumptions could embryo selection be so valuable as to motivate parents without fertility problems into using IVF solely to benefit from embryo selection?

4. if it is not worthwhile because the genetic information is too weakly predictive of adult phenotype, how much additional data would it take to make the predictions good enough to make selection worthwhile?

5. What are the prospects for embryo editing instead of selection, in theory and right now?

I start with Shulman & Bostrom2014’s basic framework, replicate it, and extend it to include realistic parameters for practical obstacles & inefficiencies, full cost-benefits, and extensions & possible improvements to the naive univariate embryo selection approach, among other things. (A subsequent 2019 analysis, Karavani et al 2019 (code/supplement), while concluding that the glass is half-empty, reaches similar results within its self-imposed analytical limits; their 2021 followup paper recapitulates my liability-threshold results about traits like schizophrenia, etc. (For further glass-half-empty takes, see Turley et al 2021.) Similarly, Treff et al 2019/Treff et al 2020/Tellier et al 2021. These largely recapitulate the expected results from the many sibling PGS comparison studies discussed later, such as Lello et al 2020.)

gcta <- read.csv(stdin(), header=TRUE)
Study, N, HSNP, SE, Age.mean, Twin, Country
Davies et al 2011, 2139, 0.17, 0.2, 70.82, FALSE, Scotland
Davies et al 2011, 6063, 0.99, 0.22, 65, FALSE, England
Plomin et al 2013, 3154, 0.35, 0.117, 12, TRUE, England
Benyamin et al 2014, 3376, 0.40, 0.21, 14, TRUE, USA
Benyamin et al 2014, 5517, 0.46, 0.06, 9, FALSE, England
Rietveld et al 2013, 7959, 0.224, 0.042, 57.47, FALSE, international
Marioni et al 2014, 6609, 0.29, 0.05, 57, FALSE, Scotland
Kirkpatrick et al 2014, 3322, 0.35, 0.11, 14.63, FALSE, USA
Toro et al 2014, 1765, 0.54, 0.25, 14.5, FALSE, international
Davies et al 2015, 6617, 0.29, 0.05, 57.2, FALSE, USA
Davies et al 2015, 5976, 0.28, 0.07, 70, FALSE, USA
Davies et al 2016, 30801, 0.31, 0.018, 56.91, FALSE, England
Robinson et al 2015, 3689, 0.36, 0.108, 13.7, FALSE, USA

library(metafor)
## Model as continuous normal variable; heritabilities are ratios 0-1,
## but metafor doesn't support heritability ratios, or correlations with
## standard errors rather than <em>n</em>s (which grossly overstates precision)
## so, as is common and safe when the estimates are not near 0/1, we treat it
## as a standardized mean difference
rem <- rma(measure="SMD", yi=HSNP, sei=SE, data=gcta); rem
# ...estimate       se     zval     pval    ci.lb    ci.ub
#  0.3207   0.0253  12.6586   <.0001   0.2711   0.3704
remAge <- rma(yi=HSNP, sei=SE, mods = Age.mean, data=gcta); remAge
# Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#
# tau^2 (estimated amount of residual heterogeneity):     0.0001 (SE = 0.0010)
# tau (square root of estimated tau^2 value):             0.0100
# I^2 (residual heterogeneity / unaccounted variability): 2.64%
# H^2 (unaccounted variability / sampling variability):   1.03
# R^2 (amount of heterogeneity accounted for):            96.04%
#
# Test for Residual Heterogeneity:
# QE(df = 11) = 15.6885, p-val = 0.1531
#
# Test of Moderators (coefficient(s) 2):
# QM(df = 1) = 6.6593, p-val = 0.0099
#
# Model Results:
#
#          estimate      se     zval    pval    ci.lb    ci.ub
# intrcpt    0.4393  0.0523   8.3953  <.0001   0.3368   0.5419
# mods      −0.0025  0.0010  -2.5806  0.0099  −0.0044  −0.0006
remAgeT <- rma(yi=HSNP, sei=SE, mods = ~ Age.mean + Twin, data=gcta); remAgeT
# intrcpt      0.4505  0.0571   7.8929  <.0001   0.3387   0.5624
# Age.mean    −0.0027  0.0010  -2.5757  0.0100  −0.0047  −0.0006
# Twin TRUE   −0.0552  0.1119  −0.4939  0.6214  −0.2745   0.1640
gcta <- gcta[order(gcta$Age.mean),] # sort by age, young to old
forest(rma(yi=HSNP, sei=SE, data=gcta), slab=gcta$Study)
## so estimated heritability at 30yo:
0.4505 + 30*-0.0027
# [1] 0.3695
## Take a look at the possible existence of a quadratic trend as suggested
## by conventional IQ heritability results:
remAgeTQ <- rma(yi=HSNP, sei=SE, mods = ~ I(Age.mean^2) + Twin, data=gcta); remAgeTQ
# Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#
# tau^2 (estimated amount of residual heterogeneity):     0.0000 (SE = 0.0009)
# tau (square root of estimated tau^2 value):             0.0053
# I^2 (residual heterogeneity / unaccounted variability): 0.83%
# H^2 (unaccounted variability / sampling variability):   1.01
# R^2 (amount of heterogeneity accounted for):            98.87%
#
# Test for Residual Heterogeneity:
# QE(df = 10) = 16.1588, p-val = 0.0952
#
# Test of Moderators (coefficient(s) 2,3):
# QM(df = 2) = 6.2797, p-val = 0.0433
#
# Model Results:
#
#                estimate      se     zval    pval    ci.lb    ci.ub
# intrcpt          0.4150  0.0457   9.0879  <.0001   0.3255   0.5045
# I(Age.mean^2)   −0.0000  0.0000  -2.4524  0.0142  −0.0001  −0.0000
# Twin TRUE       −0.0476  0.1112  −0.4285  0.6683  −0.2656   0.1703
## does fit better but enough?

exactMax <- Vectorize(function (n, mean=0, sd=1) {
if (n>2000) { ## avoid lmomco bugs at higher _n_, where the approximations are near-exact anyway
    chen1999 <- function(n,mean=0,sd=1){ mean + qnorm(0.5264^(1/n), sd=sd) }
    chen1999(n,mean=mean,sd=sd) } else {
    if(n>200) { library(lmomco)
        exactMax_unmemoized <- function(n, mean=0, sd=1) {
            expect.max.ostat(n, para=vec2par(c(mean, sd), type="nor"), cdf=cdfnor, pdf=pdfnor) }
        exactMax_unmemoized(n,mean=mean,sd=sd) } else {

 lookup <- c(0,0,0.5641895835,0.8462843753,1.0293753730,1.1629644736,1.2672063606,1.3521783756,1.4236003060,
  1.4850131622,1.5387527308,1.5864363519,1.6292276399,1.6679901770,1.7033815541,1.7359134449,1.7659913931,
  1.7939419809,1.8200318790,1.8444815116,1.8674750598,1.8891679149,1.9096923217,1.9291617116,1.9476740742,
  1.9653146098,1.9821578398,1.9982693020,2.0137069241,2.0285221460,2.0427608442,2.0564640976,2.0696688279,
  2.0824083360,2.0947127558,2.1066094396,2.1181232867,2.1292770254,2.1400914552,2.1505856577,2.1607771781,
  2.1706821847,2.1803156075,2.1896912604,2.1988219487,2.2077195639,2.2163951679,2.2248590675,2.2331208808,
  2.2411895970,2.2490736293,2.2567808626,2.2643186963,2.2716940833,2.2789135645,2.2859833005,2.2929091006,
  2.2996964480,2.3063505243,2.3128762306,2.3192782072,2.3255608518,2.3317283357,2.3377846191,2.3437334651,
  2.3495784520,2.3553229856,2.3609703096,2.3665235160,2.3719855541,2.3773592389,2.3826472594,2.3878521858,
  2.3929764763,2.3980224835,2.4029924601,2.4078885649,2.4127128675,2.4174673530,2.4221539270,2.4267744193,
  2.4313305880,2.4358241231,2.4402566500,2.4446297329,2.4489448774,2.4532035335,2.4574070986,2.4615569196,
  2.4656542955,2.4697004768,2.4736966781,2.4776440650,2.4815437655,2.4853968699,2.4892044318,2.4929674704,
  2.4966869713,2.5003638885,2.5039991455,2.5075936364,2.5111482275,2.5146637581,2.5181410417,2.5215808672,
  2.5249839996,2.5283511812,2.5316831323,2.5349805521,2.5382441196,2.5414744943,2.5446723168,2.5478382097,
  2.5509727783,2.5540766110,2.5571502801,2.5601943423,2.5632093392,2.5661957981,2.5691542321,2.5720851410,
  2.5749890115,2.5778663175,2.5807175211,2.5835430725,2.5863434103,2.5891189625,2.5918701463,2.5945973686,
  2.5973010263,2.5999815069,2.6026391883,2.6052744395,2.6078876209,2.6104790841,2.6130491728,2.6155982225,
  2.6181265612,2.6206345093,2.6231223799,2.6255904791,2.6280391062,2.6304685538,2.6328791081,2.6352710490,
  2.6376446504,2.6400001801,2.6423379005,2.6446580681,2.6469609341,2.6492467445,2.6515157401,2.6537681566,
  2.6560042252,2.6582241720,2.6604282187,2.6626165826,2.6647894763,2.6669471086,2.6690896839,2.6712174028,
  2.6733304616,2.6754290533,2.6775133667,2.6795835873,2.6816398969,2.6836824739,2.6857114935,2.6877271274,
  2.6897295441,2.6917189092,2.6936953850,2.6956591311,2.6976103040,2.6995490574,2.7014755424,2.7033899072,
  2.7052922974,2.7071828562,2.7090617242,2.7109290393,2.7127849375,2.7146295520,2.7164630139,2.7182854522,
  2.7200969934,2.7218977622,2.7236878809,2.7254674700,2.7272366478,2.7289955308,2.7307442335,2.7324828686,
  2.7342115470,2.7359303775,2.7376394676,2.7393389228,2.7410288469,2.7427093423,2.7443805094,2.7460424475)

 return(mean + sd*lookup[n+1]) }}})

exactMax(4) * sqrt(0.05)
# [1] 0.230175331
exactMax(2) * sqrt(0.10)
# [1] 0.178412412
0.230175331 / 0.178412412
# [1] 1.29013071

## show the locations of expected maxima/minima, demonstrating diminishing returns/thin tails:
x <- seq(-3, 3, length=1000)
y <- dnorm(x, mean=0, sd=1)
extremes <- unlist(Map(exactMax, 1:100))
plot(x, y, type="l", lwd=2,
    xlab="SDs", ylab="Normal density", main="Expected maximum/minimums for Gaussian samples of size n=1-100")
abline(v=c(extremes, extremes*-1), col=rep(c("black","gray"), 200))

exactMax(10, mean=0, sd=1)
# [1] 1.53875273
exactMax(10, mean=1, sd=1)
# [1] 2.53875273
exactMax(10, mean=0, sd=2)
# [1] 3.07750546

library(ggplot2)
library(gridExtra)

compareDistributions <- function(n=10, varianceMultiplier=2) {
    baselineMax <- exactMax(n, mean=0, sd=1)
    increasedVarianceMax <- exactMax(n, mean=0, sd=varianceMultiplier)
    baselineAdjusted <- increasedVarianceMax - baselineMax

    width <- increasedVarianceMax*1.2
    x1 <- seq(-width, width, length=1000)
    y1 <- dnorm(x1, mean=baselineAdjusted, sd=1)

    x2 <- seq(-width, width, length=1000)
    y2 <- dnorm(x2, mean=0, sd=varianceMultiplier)

    df <- data.frame(X=c(x1, x2), Y=c(y1, y2), Distribution=c(rep("baseline", 1000), rep("variable", 1000)))

    return(qplot(X, Y, color=Distribution, data=df) +
        geom_vline(xintercept=increasedVarianceMax, color="blue") +
        ggtitle(paste0("Variance Increase: ", varianceMultiplier, "x (Difference: +",
            round(digits=2, baselineAdjusted), "SD)")) +
        geom_text(aes(x=increasedVarianceMax*1.01, label=paste0("expected maximum (n=", n, ")"),
            y=0.3), colour="blue", angle=270))
}
p0 <- compareDistributions(varianceMultiplier=1.25) +
 ggtitle("Mean increase required to have equal expected maximum as a more\
 variable distribution\nVariance increase: 1.25x (Difference: +0.38SD)")
p1 <- compareDistributions(varianceMultiplier=1.50)
p2 <- compareDistributions(varianceMultiplier=1.75)
p3 <- compareDistributions(varianceMultiplier=2.00)
p4 <- compareDistributions(varianceMultiplier=3.00)
p5 <- compareDistributions(varianceMultiplier=4.00)
p6 <- compareDistributions(varianceMultiplier=5.00)
grid.arrange(p0, p1, p2, p3, p4, p5, p6, ncol=1)

## select 1 out of N embryos (default: siblings, who are half-related)
embryoSelection <- function(n, variance=1/3, relatedness=1/2) {
    exactMax(n, mean=0, sd=sqrt(variance*relatedness)); }
embryoSelection(n=10) * 15
# [1] 9.422897577
embryoSelection(n=10, variance=0.444) * 15
# [1] 10.87518323
embryoSelection(n=5, variance=0.444) * 15
# [1] 8.219287927

embryoSelection(n=10, variance=0.8) * 15
# [1] 14.59789016

library(MASS)
library(ggplot2)
plotSelection <- function(n, variance, relatedness=1/2) {
    r = sqrt(variance*relatedness)

    data = mvrnorm(n=n, mu=c(0, 0), Sigma=matrix(c(1, r, r, 1), nrow=2), empirical=TRUE)
    df <- data.frame(Trait=data[,1], PGS=data[,2], Selected=max(data[,2]) == data[,2])

    trueMax <- max(df$Trait)
    selected <- df[df$Selected,]$Trait
    regret <- trueMax - selected

    return(qplot(PGS, Trait, color=Selected, size=I(9), data=df) +
        coord_cartesian(ylim = c(-2.5,2.5), xlim=c(-2.5,2.5), expand=FALSE) +
        geom_hline(yintercept=0, color="red") +
        labs(title=paste0("Selection hypothetical (higher=better): with n=", n, " samples & PGS variance=", round(variance,digits=2),
            ". Performance: true max: ", round(trueMax, digits=2), "; selected: ", round(selected, digits=2),
            "; regret: ", round(regret, digits=2)))
     )
}
library(animation)
saveGIF(
    for (i in 1:100) {
      n   <- max(3, round(rnorm(1, mean=6, sd=3)))
      pgs <- runif(1, min=0, max=0.5)
      p <- plotSelection(n, pgs)
      print(p)
      },
    interval=0.8, ani.width = 1000, ani.height=800,
    movie.name = "embryo-selection.gif")

## consider first column as true latent genetic scores, & the second column as noisy measurements correlated _r_:
library(MASS)
generateCorrelatedNormals <- function(n, r) {
    mvrnorm(n=n, mu=c(0, 0), Sigma=matrix(c(1, r, r, 1), nrow=2)) }

## consider plausible scenarios for IQ-related non-massive simple embryo selection, so 2-50 embryos;
## and PGSes must max out by 80%:
scenarios <- expand.grid(Embryo.n=2:50, PGS.variance=seq(0.01, 0.50, by=0.02), Rank.mean=NA, P.max=NA, P.min=NA, P.below.mean=NA, P.minus.two=NA, Regret.SD=NA)
for (i in 1:nrow(scenarios)) {
 n = scenarios[i,]$Embryo.n
 r = sqrt(scenarios[i,]$PGS.variance * 0.5) # relatedness deflation for the ES context
 iters = 500000
 sampleStatistics <- function(n,r) {
     sim <- generateCorrelatedNormals(n, r=r)
     # max1_l  <- max(sim[,1])
     # max2_m  <- max(sim[,2])
     max1i_l <- which.max(sim[,1])
     max2i_m <- which.max(sim[,2])
     gain <- sim[,1][max2i_m]
     rank <- which(sim[,2][max1i_l] == sort(sim[,2]))

     ## P(max): if the max of the noisy measurements is a different index than the max or min of the true latents,
     ## then embryo selection fails to select the best/maximum or selects the worst.
     ## If n=1, trivially P.max/P.min=1 & Regret=0; if r=0, P.max/P.min = 1/n;
     ## if r=1, P.max=1 & P.min=0; r=0-1 can be estimated by simulation:
     P.max <- max2i_m == max1i_l
     ## P(min): if our noisy measurement led us to select the worst point rather than best:
     P.min <- which.min(sim[,1]) == max2i_m
     ## P(<avg): whether we managed to at least boost above mean of 0
     P.below.mean <- gain < 0
     ## P(IQ(70)): whether the point falls below -2SDs
     P.minus.two <- gain <= -2

     ## Regret is the difference between the true latent's maximum, and the true score
     ## for the index with the maximum of the noisy measurements, which if a different index,
     ## means a loss and thus non-zero regret.
     ## If r=0, regret = max/the n_k order statistic; r=1, regret=0; in between, simulation:
     Regret.SD <- max(gain,0)
     return(c(P.max, P.min, P.below.mean, P.minus.two, Regret.SD, rank)) }
 sampleAverages <- colMeans(t(replicate(iters, sampleStatistics(n,r))))
 # print(c(n,r,sampleAverages))
 scenarios[i,]$P.max        <- sampleAverages[1]
 scenarios[i,]$P.min        <- sampleAverages[2]
 scenarios[i,]$P.below.mean <- sampleAverages[3]
 scenarios[i,]$P.minus.two  <- sampleAverages[4]
 scenarios[i,]$Regret.SD    <- sampleAverages[5]
 scenarios[i,]$Rank.mean    <- sampleAverages[6]
 }

library(ggplot2); library(gridExtra)
p0 <- qplot(Embryo.n, Rank.mean, color=as.ordered(PGS.variance), data=scenarios) +
    theme(legend.title=element_blank()) + geom_abline(slope=1, intercept=0) +
    ggtitle("Expected true rank after selecting best out of N embryos based on PGS score (idealized, excluding IVF losses)")
p1 <- qplot(Embryo.n, P.max, color=as.ordered(PGS.variance), data=scenarios) +
    coord_cartesian(ylim = c(0,0.84)) + theme(legend.title=element_blank()) +
    ggtitle("Probability of selecting best out of N embryos as function of PGS score (*)")
p2 <- qplot(Embryo.n, P.min, color=as.ordered(PGS.variance), data=scenarios) +
    coord_cartesian(ylim = c(0,0.48)) + theme(legend.title=element_blank()) +
    ggtitle("Probability of mistakenly selecting worst out of N embryos (*)")
p3 <- qplot(Embryo.n, P.below.mean, color=as.ordered(PGS.variance), data=scenarios) +
    coord_cartesian(ylim = c(0,0.5)) + theme(legend.title=element_blank()) +
    ggtitle("Probability of mistakenly selecting below-average out of N embryos (*)")
p4 <- qplot(Embryo.n, P.minus.two, color=as.ordered(PGS.variance), data=scenarios) +
    coord_cartesian(ylim = c(0,0.02)) + theme(legend.title=element_blank()) +
    ggtitle("Probability of selecting below -2SDs out of N embryos (*)")
p5 <- qplot(Embryo.n, Regret.SD, color=as.ordered(PGS.variance), data=scenarios) +
    theme(legend.title=element_blank()) +
    ggtitle("Loss from non-omniscient selection from N embryos in SDs (*)")
grid.arrange(p0, p1, p2, p3, p4, p5, ncol=1)

selzam2016 <- 0.035
embryoSelection(n=10, variance=selzam2016) * 15
# [1] 3.053367791

# https://www.genome.gov/pages/der/seqcost2015_4.xlsx
genome <- c(9408739,9047003,8927342,7147571,3063820,1352982,752080,342502,232735,154714,108065,70333,46774,
            31512,31125,29092,20963,16712,10497,7743,7666,5901,5985,6618,5671,5826,5550,5096,4008,4920,4905,
            5731,3970,4211,1363,1245,1000)
l <- lm(log(I(tail(genome, 3))) ~ I(1:3)); l
# Coefficients:
# (Intercept)       I(1:3)
#  7.3937180   −0.1548441
exp(sapply(1:10, function(t) { 7.3937180 + -0.1548441*t } )) / 6
# [1] 232.08749421 198.79424215 170.27695028 145.85050092 124.92805739 107.00696553  91.65667754
# [8]  78.50840827  67.24627528  57.59970987

## simulate a single IVF cycle (which may not yield any live birth, in which case there is no gain returnable):
simulateIVF <- function (eggMean, eggSD, polygenicScoreVariance, normalityP=0.5, vitrificationP, liveBirth) {
  eggsExtracted <- max(0, round(rnorm(n=1, mean=eggMean, sd=eggSD)))

  normal        <- rbinom(1, eggsExtracted, prob=normalityP)

  scores        <- rnorm(n=normal, mean=0, sd=sqrt(polygenicScoreVariance*0.5))

  survived      <- Filter(function(x){rbinom(1, 1, prob=vitrificationP)}, scores)

  selection <- sort(survived, decreasing=TRUE)

  if (length(selection)>0) {
   for (embryo in 1:length(selection)) {
    if (rbinom(1, 1, prob=liveBirth) == 1) {
      live <- selection[embryo]
      return(live)
     }
    }
   }
 }
simulateIVFs <- function(eggMean, eggSD, polygenicScoreVariance, normalityP, vitrificationP, liveBirth, iters=100000) {
  return(unlist(replicate(iters, simulateIVF(eggMean, eggSD, polygenicScoreVariance, normalityP, vitrificationP, liveBirth)))); }

simulateTan <- function() { return(simulateIVFs(3, 4.6, selzam2016, 0.5, 0.96, 0.24)); }
iqTan <- mean(simulateTan()) * 15; iqTan
# [1] 0.3808377013

simulateHodesWertz <- function() { return(simulateIVFs(8.2, 4.6, selzam2016, 0.35, 0.96, 0.40)) }
iqHW <- mean(simulateHodesWertz()) * 15; iqHW
# [1] 0.684226242

dnorm((130-100)/15) * 320000000
# [1] 17277109.28
dnorm((160-100)/15) * 320000000
# [1] 42825.67224

size <- function(mean, cutoff, populationSize, useFraction=1) { if(cutoff>mean) { dnorm(cutoff-mean) * populationSize * useFraction } else
                                                                 { (1 - dnorm(cutoff-mean)) * populationSize * useFraction }}
size(0, (60/15), 67996)
# [1] 9.099920031

eliteGain <- function(ivfMean, ivfGain, ivfFraction, generation, cutoff, ivfPop, genMean, genPop) {

              ivfers      <- size(ivfMean,                      cutoff, ivfPop, 1)
              selected    <- size(ivfMean+(ivfGain*generation), cutoff, ivfPop, ivfFraction)
              nonSelected <- size(ivfMean,                      cutoff, ivfPop, 1-ivfFraction)
              gain        <- (selected+nonSelected) - ivfers

              population <- size(genMean, cutoff, genPop)
              multiplier <- gain / population
              return(multiplier) }
eliteGain(0, (9.36/15), 1, 1, (60/15), 67996, 0, 3932181)
# [1] 0.1554096565

eliteGain(0, (9.36/15), 1, 2, (60/15), 67996, 0, 3932181)
# [1] 1.151238772
eliteGain(0, (9.36/15), 1, 5, (60/15), 67996, 0, 3932181)
# [1] 34.98100356

eliteGain(0, (0.13808892057/15), 1, 1, (60/15), 67996, 0, 3932181)
# [1] 0.0006478714323
eliteGain((15/15), (0.13808892057/15), 1, 1, (60/15), 67996, 0, 3932181)
# [1] 0.01601047464

scenarios <- expand.grid(c(0.025, 0.1, 0.9), c(4/15, 12/15, 19/15, 100/15), c(1,2), c(30/15, 60/15))
colnames(scenarios) <- c("Adoption.fraction", "IQ.gain", "Generation", "Eliteness")
scenarios$Gain.fraction <- round(do.call(mapply, c(function(adoptionRate, gain, generation, selectiveness) {
                                  eliteGain(0, gain, adoptionRate, generation, selectiveness, 3932181, 0, 3932181) }, unname(scenarios[,1:4]))),
                           digits=2)

simulateIVFCB <- function (eggMean, eggSD, polygenicScoreVariance, normalityP=0.5, vitrificationP, liveBirth, fixedCost, embryoCost, traitValue) {
  eggsExtracted <- max(0, round(rnorm(n=1, mean=eggMean, sd=eggSD)))

  normal        <- rbinom(1, eggsExtracted, prob=normalityP)

  totalCost     <- fixedCost + normal * embryoCost
  scores        <- rnorm(n=normal, mean=0, sd=sqrt(polygenicScoreVariance*0.5))

  survived      <- Filter(function(x){rbinom(1, 1, prob=vitrificationP)}, scores)

  selection <- sort(survived, decreasing=FALSE)
  live <- 0
  gain <- 0

  if (length(selection)>0) {
   for (embryo in 1:length(selection)) {
    if (rbinom(1, 1, prob=liveBirth) == 1) {
      live <- selection[embryo]
      }
   }
  gain <- max(0, live - mean(selection))
  }
  return(data.frame(Trait.SD=gain, Cost=totalCost, Net=(traitValue*gain - totalCost)))  }
library(plyr)
simulateIVFCBs <- function(eggMean, eggSD, polygenicScoreVariance, normalityP, vitrificationP, liveBirth, fixedCost, embryoCost, traitValue, iters=20000) {
  ldply(replicate(simplify=FALSE, iters, simulateIVFCB(eggMean, eggSD, polygenicScoreVariance, normalityP, vitrificationP, liveBirth, fixedCost, embryoCost, traitValue))) }

iqLow <- 3270*15; iqHigh <- 16151*15
## Tan:
summary(simulateIVFCBs(3, 4.6, selzam2016, 0.5, 0.96, 0.24, 1500, 200, iqLow))
#    Trait.SD               Cost              Net
# Min.   :0.00000000   Min.   :1500.00   Min.   :-3900.0000
# 1st Qu.:0.00000000   1st Qu.:1500.00   1st Qu.:-1700.0000
# Median :0.00000000   Median :1700.00   Median :-1500.0000
# Mean   :0.02854686   Mean   :1873.05   Mean   : -472.8266
# 3rd Qu.:0.03149430   3rd Qu.:2100.00   3rd Qu.: -579.1553
# Max.   :0.42872383   Max.   :4300.00   Max.   :19076.2182
summary(simulateIVFCBs(3, 4.6, selzam2016, 0.5, 0.96, 0.24, 1500, 200, iqHigh))
#    Trait.SD               Cost              Net
# Min.   :0.00000000   Min.   :1500.00   Min.   : -4100.000
# 1st Qu.:0.00000000   1st Qu.:1500.00   1st Qu.: -1700.000
# Median :0.00000000   Median :1700.00   Median : -1500.000
# Mean   :0.02847819   Mean   :1873.08   Mean   :  5026.188
# 3rd Qu.:0.03005473   3rd Qu.:2100.00   3rd Qu.:  5143.879
# Max.   :0.48532430   Max.   :4100.00   Max.   :115677.092

## Hodes-Wertz:
summary(simulateIVFCBs(8.2, 4.6, selzam2016, 0.35, 0.96, 0.40, 1500, 200, iqLow))
#    Trait.SD                Cost              Net
# Min.   :0.000000000   Min.   :1500.00   Min.   :-4100.0000
# 1st Qu.:0.000000000   1st Qu.:1700.00   1st Qu.:-1900.0000
# Median :0.007840085   Median :2100.00   Median :-1500.0000
# Mean   :0.051678465   Mean   :2079.25   Mean   :  455.5787
# 3rd Qu.:0.090090594   3rd Qu.:2300.00   3rd Qu.: 2168.2666
# Max.   :0.463198015   Max.   :4100.00   Max.   :21019.8626
summary(simulateIVFCBs(8.2, 4.6, selzam2016, 0.35, 0.96, 0.40, 1500, 200, iqHigh))
#    Trait.SD                Cost              Net
# Min.   :0.000000000   Min.   :1500.00   Min.   : -3700.0000
# 1st Qu.:0.000000000   1st Qu.:1700.00   1st Qu.: -1700.0000
# Median :0.006228574   Median :2100.00   Median :  -650.2792
# Mean   :0.050884913   Mean   :2083.41   Mean   : 10244.2234
# 3rd Qu.:0.088152844   3rd Qu.:2300.00   3rd Qu.: 19048.4272
# Max.   :0.486235107   Max.   :4100.00   Max.   :114497.7483
## USA, youngest:
summary(simulateIVFCBs(9, 4.6, selzam2016, 0.3, 0.90, 10.8/100, 1500, 200, iqLow))
#    Trait.SD               Cost              Net
# Min.   :0.00000000   Min.   :1500.00   Min.   :-3900.0000
# 1st Qu.:0.00000000   1st Qu.:1700.00   1st Qu.:-2045.5047
# Median :0.00000000   Median :1900.00   Median :-1500.0000
# Mean   :0.03360950   Mean   :2037.22   Mean   : -388.6739
# 3rd Qu.:0.05023528   3rd Qu.:2300.00   3rd Qu.:  287.3619
# Max.   :0.52294123   Max.   :3900.00   Max.   :23950.2672
summary(simulateIVFCBs(9, 4.6, selzam2016, 0.3, 0.90, 10.8/100, 1500, 200, iqHigh))
#    Trait.SD               Cost              Net
# Min.   :0.00000000   Min.   :1500.00   Min.   : -3900.000
# 1st Qu.:0.00000000   1st Qu.:1700.00   1st Qu.: -1900.000
# Median :0.00000000   Median :1900.00   Median : -1500.000
# Mean   :0.03389909   Mean   :2044.75   Mean   :  6167.812
# 3rd Qu.:0.05115755   3rd Qu.:2300.00   3rd Qu.: 10224.781
# Max.   :0.45364794   Max.   :4100.00   Max.   :108203.019

summary(simulateIVFCBs(9, 4.6, 0.33, 0.3, 0.90, 10.8/100, 1500, 200, iqLow))
#     Trait.SD              Cost              Net
#  Min.   :0.0000000   Min.   :1500.00   Min.   :-3700.000
#  1st Qu.:0.0000000   1st Qu.:1700.00   1st Qu.:-2100.000
#  Median :0.0000000   Median :1900.00   Median :-1500.000
#  Mean   :0.1037614   Mean   :2042.24   Mean   : 3047.259
#  3rd Qu.:0.1562492   3rd Qu.:2300.00   3rd Qu.: 5516.869
#  Max.   :1.4293926   Max.   :3900.00   Max.   :68411.709
summary(simulateIVFCBs(9, 4.6, 0.33, 0.3, 0.90, 10.8/100, 1500, 200, iqHigh))
#     Trait.SD              Cost             Net
#  Min.   :0.0000000   Min.   :1500.0   Min.   : -4100.00
#  1st Qu.:0.0000000   1st Qu.:1700.0   1st Qu.: -1900.00
#  Median :0.0000000   Median :1900.0   Median : -1500.00
#  Mean   :0.1030492   Mean   :2037.6   Mean   : 22927.61
#  3rd Qu.:0.1530295   3rd Qu.:2300.0   3rd Qu.: 34652.62
#  Max.   :1.3798166   Max.   :4100.0   Max.   :331981.26
ivfBirths <- 67996; discount <- 0.05
current <- 6230; perfect <- 23650
(ivfBirths * perfect)/(log(1+discount)) - (ivfBirths * current)/(log(1+discount))
# [1] 24277235795

gwasSizeToVariance <- function(N, h2) { ((N / 67865) * h2^2) / ((N/67865) * h2 + 1) }
sampleIncrease <- 1000
original     <- gwasSizeToVariance(17000, 0.33)
originalplus <- gwasSizeToVariance(17000+sampleIncrease, 0.33)
originalGain     <- mean(simulateIVFCBs(9, 4.6, original, 0.3, 0.90, 10.8/100, 1500, 200, iqHigh)$Net)
originalplusGain <- mean(simulateIVFCBs(9, 4.6, originalplus, 0.3, 0.90, 10.8/100, 1500, 200, iqHigh)$Net)
originalGain; originalplusGain
((((originalplusGain - originalGain) * ivfBirths) / log(1+discount)) / sampleIncrease) / 6
# [1] 71716.90116

gainByEggcount <- sapply(1:300, function(egg) { mean(simulateIVFCBs(egg, 4.6, selzam2016, 0.3, 0.90, 10.8/100, 1500, 200, iqHigh)$Net) })
max(gainByEggcount); which.max(gainByEggcount)
# [1] 26657.1117
# [1] 281
plot(1:300, gainByEggcount, xlab="Average number of eggs available", ylab="Profit")
summary(simulateIVFCBs(which.max(gainByEggcount), 4.6, selzam2016, 0.3, 0.90, 10.8/100, 1500, 200, iqHigh))
#     Trait.SD              Cost              Net
#  Min.   :0.0000000   Min.   :12300.0   Min.   :-21900.00
#  1st Qu.:0.1284192   1st Qu.:17300.0   1st Qu.: 12711.92
#  Median :0.1817688   Median :18300.0   Median : 25630.74
#  Mean   :0.1845060   Mean   :18369.1   Mean   : 26330.25
#  3rd Qu.:0.2372748   3rd Qu.:19500.0   3rd Qu.: 39162.75
#  Max.   :0.5661427   Max.   :25300.0   Max.   :117856.55
max(gainByEggcount) / which.max(gainByEggcount)
# [1] 94.86516619

allegrini2018 <- sqrt(0.11*0.5)
iqLow         <- 3270*15
testCost      <- 200

exactMax(5)
# [1] 1.162964474
exactMax(5) * iqLow
# [1] 57043.40743
(exactMax(5) - exactMax(3))
# [1] 0.3166800983
(exactMax(5) - exactMax(3)) * iqLow
# [1] 15533.15882

round(sapply(seq(2, 300, by=10), function(n) { (exactMax(n) - exactMax(n-1)) * iqLow }))
#  [1] 27673  2099  1007   648   473   370   303   255   220   194   172   155   141   129   119
#        110   103    96    90    85    81    76    73    69    66    63    61
# [28]    58    56    54

round(sapply(2:26, function(n) { (exactMax(n, sd=allegrini2018) - exactMax(n-1, sd=allegrini2018)) * iqLow }))
# [1] 6490 3245 2106 1537 1199  977  822  706  618  549  492  446  407  374  346  322  300  281  265  250  236  224  213  203  194
exactMax(24, sd=allegrini2018)
# [1] 0.4567700586
exactMax(25, sd=allegrini2018)
# [1] 0.4609071309
0.4609071309 - 0.4567700586
# [1] −0.0041370723

stoppingRule <- function(predictorSD, utilityCost, utilityGain) {
 n <- 1
 while(((exactMax(n+1, sd=predictorSD) - exactMax(n, sd=predictorSD)) * utilityGain) > utilityCost) { n <- n+1 }
 return(c(n, exactMax(n), exactMax(n, sd=predictorSD))) }

round(digits=2, stoppingRule(allegrini2018, testCost, iqLow))
# [1] 25.00 1.97 0.46
round(digits=2, stoppingRule(allegrini2018, 100, iqLow))
# [1] 45.00  2.21  0.52

marginalGain <- (exactMax(48, sd=allegrini2018) - exactMax(24, sd=allegrini2018)) * iqLow
marginalCost <- (48-24) * testCost
marginalGain; marginalCost
[1] 3283.564451
[1] 4800
marginalGain - marginalCost
[1] -1516.435549

library(memoise)
library(parallel) # warning, Windows users
library(plyr)

## Memoize the Monte Carlo evaluation to save time - it's almost exact w/100k & simpler:
expectedPastThreshold <- memoise(function(maximum, predictorSD) {
    mean({ x <- rnorm(100000, sd=predictorSD); ifelse(x>maximum, x-maximum, 0) }) })

optimalSearch <- function(maxN, predictorSD, utilityCost, utilityBenefit) {

    samples <- rnorm(maxN, sd=predictorSD)

    i <- 1; maximum <- samples[1]; cost <- utilityCost; profit <- 0; gain <- max(maximum,0);
    while (i < maxN) {

        marginalGain <- expectedPastThreshold(maximum, predictorSD)

        if (marginalGain*utilityBenefit > utilityCost) {
          i <- i+1
          cost <- cost+utilityCost
          nth <- samples[i]
          maximum <- max(maximum, nth); } else { break; }
    }

    gain <- maximum * utilityBenefit; profit <- gain-cost;
    searchAllProfit <- max(samples)*utilityBenefit - maxN*utilityCost

    return(c(i, maximum, cost, gain, profit, searchAllProfit, searchAllProfit - (gain-cost)))
    }

optimalSearch(100, allegrini2018, testCost, iqLow)
# [1]    48     0  9600 22475 12875  9462 -3413

## Parallelize simulations:
optimalSearchs <- function(a,b,c,d, iters=10000) { df <- ldply(mclapply(1:iters, function(x) { optimalSearch(a,b,c,d); }));
  colnames(df) <- c("N", "Maximum.SD", "Cost.total", "Gain.total", "Profit", "Nonadaptive.profit", "Nonadaptivity.regret"); return(df) }

summary(digits=2, optimalSearchs(5,   allegrini2018, testCost, iqLow))
#       N         Maximum.SD      Cost.total     Gain.total         Profit       Nonadaptive.profit Nonadaptivity.regret
# Min.   :1.0   Min.   :-0.27   Min.   : 200   Min.   :-13039   Min.   :-14039   Min.   :-14039     Min.   : -800
# 1st Qu.:5.0   1st Qu.: 0.16   1st Qu.:1000   1st Qu.:  7978   1st Qu.:  6978   1st Qu.:  6978     1st Qu.:    0
# Median :5.0   Median : 0.26   Median :1000   Median : 12902   Median : 11902   Median : 11902     Median :    0
# Mean   :4.6   Mean   : 0.27   Mean   : 921   Mean   : 13267   Mean   : 12346   Mean   : 12306     Mean   :  -40
# 3rd Qu.:5.0   3rd Qu.: 0.37   3rd Qu.:1000   3rd Qu.: 18199   3rd Qu.: 17199   3rd Qu.: 17199     3rd Qu.:    0
# Max.   :5.0   Max.   : 1.05   Max.   :1000   Max.   : 51405   Max.   : 51205   Max.   : 50405     Max.   :14789
summary(digits=2, optimalSearchs(10,  allegrini2018, testCost, iqLow))
#       N          Maximum.SD      Cost.total     Gain.total        Profit      Nonadaptive.profit Nonadaptivity.regret
# Min.   : 1.0   Min.   :-0.06   Min.   : 200   Min.   :-2934   Min.   :-4934   Min.   :-4934      Min.   :-1800
# 1st Qu.: 7.0   1st Qu.: 0.27   1st Qu.:1400   1st Qu.:13047   1st Qu.:11047   1st Qu.:11047      1st Qu.: -400
# Median :10.0   Median : 0.35   Median :2000   Median :17275   Median :15275   Median :15275      Median :    0
# Mean   : 8.2   Mean   : 0.36   Mean   :1649   Mean   :17594   Mean   :15945   Mean   :15754      Mean   : -190
# 3rd Qu.:10.0   3rd Qu.: 0.44   3rd Qu.:2000   3rd Qu.:21718   3rd Qu.:20742   3rd Qu.:20109      3rd Qu.:    0
# Max.   :10.0   Max.   : 0.97   Max.   :2000   Max.   :47618   Max.   :46218   Max.   :45618      Max.   :20883
summary(digits=2, optimalSearchs(24,  allegrini2018, testCost, iqLow))
#       N        Maximum.SD     Cost.total     Gain.total        Profit      Nonadaptive.profit Nonadaptivity.regret
# Min.   : 1   Min.   :0.12   Min.   : 200   Min.   : 5719   Min.   :  919   Min.   :  919      Min.   :-4600
# 1st Qu.: 7   1st Qu.:0.37   1st Qu.:1400   1st Qu.:18238   1st Qu.:13438   1st Qu.:13438      1st Qu.:-2800
# Median :16   Median :0.43   Median :3200   Median :21201   Median :19223   Median :17145      Median : -600
# Mean   :15   Mean   :0.44   Mean   :3032   Mean   :21689   Mean   :18656   Mean   :17648      Mean   :-1008
# 3rd Qu.:24   3rd Qu.:0.50   3rd Qu.:4800   3rd Qu.:24527   3rd Qu.:22636   3rd Qu.:21217      3rd Qu.:    0
# Max.   :24   Max.   :1.13   Max.   :4800   Max.   :55507   Max.   :52107   Max.   :50707      Max.   :25705
summary(digits=2, optimalSearchs(100, allegrini2018, testCost, iqLow))
#       N         Maximum.SD     Cost.total      Gain.total        Profit      Nonadaptive.profit Nonadaptivity.regret
# Min.   :  1   Min.   :0.31   Min.   :  200   Min.   :15218   Min.   :-4782   Min.   :-4782      Min.   :-19800
# 1st Qu.:  7   1st Qu.:0.43   1st Qu.: 1400   1st Qu.:21223   1st Qu.:16696   1st Qu.: 5342      1st Qu.:-15507
# Median : 16   Median :0.47   Median : 3200   Median :23239   Median :19919   Median : 8266      Median :-11772
# Mean   : 23   Mean   :0.50   Mean   : 4654   Mean   :24398   Mean   :19744   Mean   : 8762      Mean   :-10983
# 3rd Qu.: 33   3rd Qu.:0.54   3rd Qu.: 6600   3rd Qu.:26504   3rd Qu.:23076   3rd Qu.:11651      3rd Qu.: -7293
# Max.   :100   Max.   :1.10   Max.   :20000   Max.   :53952   Max.   :52352   Max.   :33952      Max.   : 18226
summary(digits=2, optimalSearchs(500, allegrini2018, testCost, iqLow))
#       N         Maximum.SD     Cost.total      Gain.total        Profit       Nonadaptive.profit Nonadaptivity.regret
# Min.   :  1   Min.   :0.40   Min.   :  200   Min.   :19607   Min.   :-25265   Min.   :-76428     Min.   :-99800
# 1st Qu.:  7   1st Qu.:0.43   1st Qu.: 1400   1st Qu.:21289   1st Qu.: 16559   1st Qu.:-67982     1st Qu.:-89569
# Median : 17   Median :0.48   Median : 3400   Median :23349   Median : 19779   Median :-65471     Median :-85154
# Mean   : 24   Mean   :0.50   Mean   : 4772   Mean   :24498   Mean   : 19726   Mean   :-64955     Mean   :-84681
# 3rd Qu.: 33   3rd Qu.:0.54   3rd Qu.: 6600   3rd Qu.:26500   3rd Qu.: 23232   3rd Qu.:-62591     3rd Qu.:-80393
# Max.   :234   Max.   :1.09   Max.   :46800   Max.   :53390   Max.   : 50453   Max.   :-44268     Max.   :-37431