Critiques

• criticism of teeth-removal experiment in rats https://www.lesswrong.com/r/discussion/lw/kfb/open_thread_30_june_2014_6_july_2014/b1u3

• criticism of small Noopept self-experiment http://www.bluelight.org/vb/threads/689936-My-Paper-quot-Noopept-amp-The-Placebo-Effect-quot?p=11910708&viewfull=1#post11910708

• why Soylent is not a good idea https://www.lesswrong.com/lw/hht/link_soylent_crowdfunding/90y7

• misinterpretation of fluoridation meta-analysis and ignorance of VoI http://theness.com/neurologicablog/index.php/anti-fluoride-propaganda-as-news/#comment-76400

• https://www.lesswrong.com/lw/1lt/case_study_melatonin/8mgf

• Fulltext: https://dl.dropboxusercontent.com/u/280585369/2014-dubal.pdf is this possible? http://nextbigfuture.com/2014/05/kl-vs-gene-makes-up-six-iq-points-of.html/Nootropics/comments/25233r/boost_your_iq_by_6_points/chddd7f

• tACS causes lucid dreaming: https://www.reddit.com/r/LucidDreaming/comments/27y7n6/no_brain_stimulation_will_not_get_you_lucid/ck6isgo

• Herbalife growth patterns: https://www.reddit.com/r/business/comments/24aoo2/what_unsustainable_growth_looks_like_herbalife/ch5hwtv

• Plausible correlate of Fairtrade: https://www.reddit.com/r/Economics/comments/26jb2d/surprise_fairtrade_doesnt_benefit_the_poor/chrx9s4

• slave whippings vs cotton production https://www.lesswrong.com/r/discussion/lw/kwc/open_thread_sept_17_2014/bajv

• whether a study on mental illness & violence shows schizophrenics are not more likely to murder but rather be murdered: https://www.reddit.com/r/psychology/comments/2fwjs8/people_with_mental_illness_are_more_likely_to_be/ckdq50k / http://www.nationalelfservice.net/publication-types/observational-study/people-with-mental-illness-are-more-likely-to-be-victims-of-homicide-than-perpetrators-of-homicide/#comment-95507 (see also http://slatestarscratchpad.tumblr.com/post/120950150581/psycholar-giraffepoliceforce-museicetc https://www.reddit.com/r/slatestarcodex/comments/744rqn/violence_is_not_a_product_of_mental_illness/dnwb1kj/ )

• Fortune analysis of higher female CEO returns https://www.lesswrong.com/r/discussion/lw/l3b/contrarian_lw_views_and_their_economic/bftw

• algae/IQ: https://www.lesswrong.com/r/discussion/lw/l9v/open_thread_nov_17_nov_23_2014/bm7o

• synaesthesia/IQ: https://www.reddit.com/r/psychology/comments/2mryte/surprising_iq_boost_12_in_average_by_a_training/cm760v8

• misinterpretation: https://slatestarcodex.com/2014/12/08/links-1214-come-ye-to-bethlinkhem/

• underpowered/multiple-correction jobs program: https://slatestarcodex.com/2014/12/08/links-1214-come-ye-to-bethlinkhem/

• claimed fall in digit span backwards minuscule and non-statistically-significant, no evidence of heterogeneity beyond variability due to sample size https://drjamesthompson.blogspot.com/2015/04/digit-span-bombshell.html?showComment=1428096775425#c4097303932864318518

• Claimed randomized experiment of whether sushi tastes worse after freezing is not actually a randomized experiment https://www.reddit.com/r/science/comments/324xmf/randomized_doubleblind_study_shows_the_quality_of/cq8dmsb

• sexual openness result undermined by ceiling effect http://mindhacks.com/2015/04/28/when-society-isnt-judging-womens-sex-drive-rivals-mens/

• music study claiming WM interaction: possible ceiling effect? see FB PM

• attempt to measure effect of Nazi anti-schizophrenia eugenics program failed to use breeder’s equation to estimate possible size of effect, which is too small to detect with available data and hence attempt is foredoomed: https://rpubs.com/JLLJ/GSCZ

• claim high IQ types almost 100% failure rates due to inappropriate model assumption of normal distribution + narrow standard deviation: https://polymatharchives.blogspot.com/2015/01/the-inappropriately-excluded.html?showComment=1441741719623

• implausible claims about success rate of facial recognition applied to St Petersburg population: https://news.ycombinator.com/item?id=11491264 (see also “Facial recognition systems stumble when confronted with million-face database”)

• human Toxoplasma gondii study is not well-powered as authors claim due to incorrect power analysis, and results are evidence for harm: https://blogs.discovermagazine.com/neuroskeptic/2016/02/20/myth-mind-altering-parasite-toxoplasma-gondii/

• attempt at attributing Bitcoin price increases to technology improvements: https://www.reddit.com/r/Bitcoin/comments/5tbt8f/buzz_factor_or_innovation_potential_what_explains/ddlmzrz/

• analysis of designer drug/research chemical activity on Wikipedia is driven almost entirely by editing patterns of just 2 Wikipedia editors particularly interested in the topic: http://calib.ro/chemical-wiki/explorations/2016-09-12-emcdda-watchlist-and-wikipedia-timeline#comment-3277669328

• failure to use mediation SEM, difference-in-statistical-significance-is-not-a-significant-difference: https://www.reddit.com/r/slatestarcodex/comments/6qwb0q/critical_thinking_skills_are_more_important_than/dl51ubw/

• Neanderthal ancestry percentage & autism: https://www.reddit.com/r/slatestarcodex/comments/74fevz/findings_suggest_that_high_levels_of_neanderthal/dny3sh9/

• Anime image classification project likely undone by using non-i.i.d. images

> Suppliers usually age the mugwort and grind it up to a fluff;

How on earth is this supposed to help AD? How does burning a plant

> Ten SRs met our inclusion criteria, which related to the

That review also mentions, incidentally, that

> Many of the primary moxibustion trials originate from China

“Someone Should Do Something”: Wishlist of Miscellaneous Project Ideas

Split out to list of ideas.

Estimating Censored Test Scores

An acquaintance asks the following question: he is applying for an university course which requires a certain minimum score on a test for admittance, and wonders about his chances and a possible trend of increasing minimum scores over time. (He hasn’t received his test results yet.) The university doesn’t provide a distribution of admittee scores, but it does provide the minimum scores for 2005^–₈2013_11ya, unless all applicants were admitted because they all scored above an unknown cutoff—in which case it provides no minimum score. This leads to the dataset:

2005,NA
2006,410
2007,NA
2008,NA
2009,398
2010,407
2011,417
2012,NA
2013,NA

A quick eyeball tells us that we can’t conclude much: only 4 actual datapoints, with 5 hidden from us. We can’t hope to conclude anything about time trends, other than there doesn’t seem to be much of one: the last score, 417, is not much higher than 410, and the last two scores are low enough to be hidden. We might be able to estimate a mean, though.

We can’t simply average the 4 scores and conclude the mean minimum is 410 because of those NAs: a number of scores have been ‘censored’ because they were too low, and while we don’t know what they were, we do know they were <398 (the smallest score) and so a bunch of <398s will pull down the uncensored mean of 410.

On approach is to treat it as a Tobit model and estimate using something like the censReg library (overview).

But if we try a quick call to censReg, we are confounded: a Tobit model expects you to provide the cutoff below which the observations were censored, but that is something we don’t know. All we know is that it must be below 398, we weren’t told it was exactly 395, 394, etc. Fortunately, this is a solved problem. For example: “The Tobit model with a non-zero threshold”, Carson & Sun 2007_17ya tells us:

In this paper, we consider estimating the unknown censoring threshold by the minimum of the uncensored y_i’s. We show that the estimator γ’ of γ is superconsistent and asymptotically exponentially distributed. Carson (1988_36ya, 1989_35ya) also suggests estimating the unknown censoring threshold by the minimum of the uncensored y_i’s. In a recent paper, Zuehlke (2003_21ya) rediscovers these unpublished results and demonstrates via simulations that the asymptotic distribution of the maximum likelihood estimator does not seem to be affected by the estimation of the censoring threshold.

That seems to be almost too simple and easy, but it makes sense and reminds me a little of the German tank problem: the minimum might not be that accurate a guess (it’s unlikely you just happened to draw a sample right on the censoring threshold) and it definitely can’t be wrong in the sense of being too low. (A Bayesian method might be able to do better with a prior like an exponential.)

With that settled, the analysis is straightforward: load the data, figure out the minimum score, set the NAs to 0, regress, and extract the model estimates for each year:

scores <- data.frame(Year=2005:2013,
                     MinimumScore=c(NA,410,NA,NA,398,407,417,NA,NA));
censorThreshold <- min(scores$MinimumScore, na.rm=T)
scores[is.na(scores)] <- 0

library(censReg)
## 'censorThreshold-1' because censReg seems to treat threshold as < and not <=
summary(censReg(MinimumScore ~ Year, left=censorThreshold-1, data=scores))
# Warning message:
# In censReg(MinimumScore ~ Year, left = censorThreshold - 1, data = scores) :
#   at least one value of the endogenous variable is smaller than the left limit
#
# Call:
# censReg(formula = MinimumScore ~ Year, left = censorThreshold -
#     1, data = scores)
#
# Observations:
#          Total  Left-censored     Uncensored Right-censored
#              9              5              4              0
#
# Coefficients:
#              Estimate Std. error t value Pr(> t)
# (Intercept) -139.9711        Inf       0       1
# Year           0.2666        Inf       0       1
# logSigma       2.6020        Inf       0       1
#
# Newton-Raphson maximisation, 37 iterations
# Return code 1: gradient close to zero
# Log-likelihood: -19.35 on 3 Df
-139.9711 + (0.2666 * scores$Year)
# [1] 394.6 394.8 395.1 395.4 395.6 395.9 396.2 396.4 396.7

With so little data the results aren’t very reliable, but there is one observation we can make.

The fact that half the dataset is censored tells us that the uncensored mean may be a huge overestimate (since we’re only looking at the ‘top half’ of the underlying data), and indeed it is. The original mean of the uncensored scores was 410; however, the estimate including the censored data is much lower, 397 (13 less)!

This demonstrates the danger of ignoring systematic biases in your data.

So, trying to calculate a mean or time effect is not helpful. What might be better is to instead exploit the censoring directly: if the censoring happened because everyone got in, then if you showed up in a censored year, you have 100% chance of getting in; while in a non-censored year you have an unknown but <100% chance of getting in; so the probability of a censored year sets a lower bound on one’s chances, and this is easy to calculate as a simple binomial problem—5 out of 9 years were censored years, so:

binom.test(c(5,4))
#
#   Exact binomial test
#
# data:  c(5, 4)
# number of successes = 5, number of trials = 9, p-value = 1
# alternative hypothesis: true probability of success is not equal to 0.5
# 95% confidence interval:
#  0.212 0.863
# sample estimates:
# probability of success
#                 0.5556

So we can tell him that he may have a >55% chance of getting in.

The Traveling Gerontologist Problem

A quick probability exercise: Wikipedia mentions Finland has 566 centenarians as of 2010_14ya.

That’s few enough you could imagine visiting them all to research them and their longevity, in a sort of traveling salesman problem but with gerontologists instead. Except, because of the exponential increase in mortality, centenarians have high annual mortality rates; it depends on the exact age but you could call it >30% (eg. Finnish 99yos in 2012 had a death toll of 326.54/1000). So you might well try to visit a centenarian and discover they’d died before you got there.

How bad a risk is this? Well, if the risk per year is 30%, then one has a 70% chance of surviving a year. To survive a year, you must survive all 365 days; by the multiplication rule, the risk is x where $0.7=x\cdot x\cdot x\cdot ...\ast x\text{[365.25 times]}$ or 0.7 = x^365.25; solving, x = 0.999024.

It takes time to visit a centenarian—it wouldn’t do to be abrupt and see them for only a few minutes, you ought to listen to their stories, and you need to get to a hotel or airport, so let’s assume you visit 1 centenarian per day.

If you visit centenarian A on day 1, and you want to visit centenarian B on day 2, then you can count on a 99.9% chance B is still alive. So far so good. And if you wanted to visit 566 centenarians (let’s imagine you have a regularly-updated master list of centenarians from the Finnish population registry), then you only have to beat the odds 566 times in a row, which is not that hard: 0.999024⁵⁶⁶ = 0.5754023437943274.

But that’s coldblooded of you to objectify those Finnish centenarians! “Any centenarian will do, I don’t care.” What if you picked the current set of 566 centenarians and wanted to visit just them, specifically—with no new centenarians introduced to the list to replace any dead ones.

That’s a little more complicated. When you visit the first centenarian, it’s the same probability: 0.999024. When you visit the second centenarian the odds change since now she (and it’s more often ‘she’ than ‘he’, since remember the exponential and males having shorter mean lifetimes) has to survive 2 days, so it’s $0.999024\cdot 0.999024$ or 0.999024²; for the third, it’s 0.999024³, and so on to #566 who has been patiently waiting and trying to survive a risk of 0.999024⁵⁶⁶, and then you need to multiply to get your odds of beating every single risk of death and the centenarian not leaving for a more permanent rendezvous: $0.999024\cdot {0.999024}^2\cdot {0.999024}^3\cdot ...\cdot {0.999024}^{566}$, which would be ${\prod }_{n=1}^{566}{0.999024}^n$, or in Haskell:

product (map (\x -> 0.999024**x) [1..566])
→ 8.952743340164081e-69

(A little surprisingly, Wolfram Alpha can solve the T_eX expression too.)

Given the use of floating point in that function (567 floating point exponentiations followed by as many multiplications) and the horror stories about floating point, one might worry the answer is wrong & the real probability is much larger. We can retry with an implementation of computable reals, CReal, which can be very slow but should give more precise answers:

:module + Data.Number.CReal
showCReal 100 (product (map (\x -> 0.999024**x) [1..566]))
→ 0.0000000000000000000000000000000000000000000000000000000000000000000089527433401308585720915431195262

Looks good—agrees with the floating point version up to the 11th digit:

8.9527433401 64081e-69
8.9527433401 308585720915431195262

We can also check by rewriting the product equation to avoid all the exponentiation and multiplication (which might cause issues) in favor of a single exponential:

1. $p^1\ast p^2\ast ...p^n$ (as before)

2. = $p^{1+2+...+n}$ (since $\left(x^m\right)\ast \left(x^n\right)=x^{(}m+n)$)

3. = $p^{\frac{n\cdot (1+n)}{2}}$ (by arithmetic progression/Gauss’s famous classroom trick since ${\sum }_1^n=n\cdot \frac{a_1+a_n}{2}$)

4. = ${0.999024}^{\frac{566\cdot (1+566)}{2}}$ (start substituting in specific values)

5. = ${0.999024}^{\frac{320922}{2}}$

6. = 0.999024¹⁶⁰⁴⁶¹

So:

0.999024^160461
→ 8.95274334014924e-69

Or to go back to the longer version:

0.999024**((566*(1 + 566)) / 2)
→ 8.952743340164096e-69

Also close. All probabilities of success are minute.

How fast would you have to be if you wanted to at least try to accomplish the tour with, say, a 50-50 chance?

Well, that’s easy: you can consider the probability of all of them surviving one day and as we saw earlier, that’s 0.999024⁵⁶⁶ = 0.58, and two days would be $({0.999024}^{566}{)}^2=0.33$ So you can only take a little over a day before you’ve probabilistically lost & one of them has died; if you hit all 566 centenarians in 24 hours, that’s ~24 centenarians per hour or ~2 minutes to chat with each one and travel to the next. If you’re trying to collect DNA samples, better hope they’re all awake and able to give consent!

So safe to say, you will probably not be able to manage the Traveling Gerontologist’s tour.

Bayes Nets

library(bnlearn)
# https://www.dropbox.com/s/4nsrszm85m47272/2015-03-22-gwern-weight.csv
weight <- read.csv("selfexperiment/weight.csv")
weight$Date <- NULL; weight$Weight.scale <- NULL
# remove missing data
weightC <- na.omit(weight)
# bnlearn can't handle integers, oddly enough
weightC <- as.data.frame(sapply(weightC, as.numeric))
summary(weightC)
#   Weight.Omron        Weight.BMI        Weight.body.fat    Weight.muscle
#  Min.   :193.0000   Min.   : 26.90000   Min.   :27.00000   Min.   :32.60000
#  1st Qu.:195.2000   1st Qu.: 27.20000   1st Qu.:28.40000   1st Qu.:34.20000
#  Median :196.4000   Median : 27.40000   Median :28.70000   Median :34.50000
#  Mean   :196.4931   Mean   : 28.95409   Mean   :28.70314   Mean   :34.47296
#  3rd Qu.:197.8000   3rd Qu.: 27.60000   3rd Qu.:29.10000   3rd Qu.:34.70000
#  Max.   :200.6000   Max.   : 28.00000   Max.   :31.70000   Max.   :35.50000
#  Weight.resting.metabolism Weight.body.age    Weight.visceral.fat
#  Min.   :1857.000          Min.   :52.00000   Min.   : 9.000000
#  1st Qu.:1877.000          1st Qu.:53.00000   1st Qu.:10.000000
#  Median :1885.000          Median :53.00000   Median :10.000000
#  Mean   :1885.138          Mean   :53.32704   Mean   : 9.949686
#  3rd Qu.:1893.000          3rd Qu.:54.00000   3rd Qu.:10.000000
#  Max.   :1914.000          Max.   :56.00000   Max.   :11.000000
cor(weightC)
#                             Weight.Omron     Weight.BMI Weight.body.fat  Weight.muscle
# Weight.Omron               1.00000000000  0.98858376919    0.1610643221 -0.06976934825
# Weight.BMI                 0.98858376919  1.00000000000    0.1521872557 -0.06231142104
# Weight.body.fat            0.16106432213  0.15218725566    1.0000000000 -0.98704369855
# Weight.muscle             -0.06976934825 -0.06231142104   -0.9870436985  1.00000000000
# Weight.resting.metabolism  0.96693236051  0.95959140245   -0.0665001241  0.15621294274
# Weight.body.age            0.82581939626  0.81286141659    0.5500409365 -0.47408608681
# Weight.visceral.fat        0.41542744168  0.43260100665    0.2798756916 -0.25076619829
#                           Weight.resting.metabolism Weight.body.age Weight.visceral.fat
# Weight.Omron                           0.9669323605    0.8258193963        0.4154274417
# Weight.BMI                             0.9595914024    0.8128614166        0.4326010067
# Weight.body.fat                       -0.0665001241    0.5500409365        0.2798756916
# Weight.muscle                          0.1562129427   -0.4740860868       -0.2507661983
# Weight.resting.metabolism              1.0000000000    0.7008354776        0.3557229425
# Weight.body.age                        0.7008354776    1.0000000000        0.4840752389
# Weight.visceral.fat                    0.3557229425    0.4840752389        1.0000000000

## create alternate dataset expressing the two percentage variables as pounds, since this might fit better
weightC2 <- weightC
weightC2$Weight.body.fat <- weightC2$Weight.Omron * (weightC2$Weight.body.fat / 100)
weightC2$Weight.muscle   <- weightC2$Weight.Omron * (weightC2$Weight.muscle / 100)

pdap <- hc(weightC)
pdapc2 <- hc(weightC2)
## bigger is better:
score(pdap, weightC)
# [1] -224.2563072
score(pdapc2, weightC2)
# [1] -439.7811072
## stick with the original, then
pdap
#   Bayesian network learned via Score-based methods
#
#   model:
#    [Weight.Omron][Weight.body.fat][Weight.BMI|Weight.Omron]
#    [Weight.resting.metabolism|Weight.Omron:Weight.body.fat]
#    [Weight.body.age|Weight.Omron:Weight.body.fat]
#    [Weight.muscle|Weight.body.fat:Weight.resting.metabolism][Weight.visceral.fat|Weight.body.age]
#   nodes:                                 7
#   arcs:                                  8
#     undirected arcs:                     0
#     directed arcs:                       8
#   average markov blanket size:           2.57
#   average neighbourhood size:            2.29
#   average branching factor:              1.14
#
#   learning algorithm:                    Hill-Climbing
#   score:                                 BIC (Gauss.)
#   penalization coefficient:              2.534452101
#   tests used in the learning procedure:  69
#   optimized:                             TRUE
plot(pdap)
## https://i.imgur.com/nipmqta.png

pdap2 <- hc(weightC, blacklist=data.frame(from=c("Weight.body.age", "Weight.resting.metabolism"), to=c("Weight.visceral.fat","Weight.muscle")))

pdap3 <- hc(weightC, blacklist=data.frame(from=c("Weight.body.age", "Weight.resting.metabolism", "Weight.BMI", "Weight.Omron"), to=c("Weight.visceral.fat","Weight.muscle", "Weight.visceral.fat", "Weight.muscle")))

pdap4 <- hc(weightC, blacklist=data.frame(from=c("Weight.body.age", "Weight.resting.metabolism", "Weight.BMI", "Weight.Omron", "Weight.Omron", "Weight.BMI"), to=c("Weight.visceral.fat","Weight.muscle", "Weight.visceral.fat", "Weight.muscle", "Weight.visceral.fat", "Weight.muscle")))

pdap5 <- hc(weightC, blacklist=data.frame(from=c("Weight.body.age", "Weight.resting.metabolism", "Weight.BMI", "Weight.Omron", "Weight.Omron", "Weight.BMI", "Weight.body.age"), to=c("Weight.visceral.fat","Weight.muscle", "Weight.visceral.fat", "Weight.muscle", "Weight.visceral.fat", "Weight.muscle", "Weight.muscle")))
#   Bayesian network learned via Score-based methods
#
#   model:
#    [Weight.body.fat][Weight.muscle|Weight.body.fat][Weight.visceral.fat|Weight.body.fat]
#    [Weight.Omron|Weight.visceral.fat][Weight.BMI|Weight.Omron]
#    [Weight.resting.metabolism|Weight.Omron:Weight.body.fat]
#    [Weight.body.age|Weight.Omron:Weight.body.fat]
#   nodes:                                 7
#   arcs:                                  8
#     undirected arcs:                     0
#     directed arcs:                       8
#   average markov blanket size:           2.57
#   average neighbourhood size:            2.29
#   average branching factor:              1.14
#
#   learning algorithm:                    Hill-Climbing
#   score:                                 BIC (Gauss.)
#   penalization coefficient:              2.534452101
#   tests used in the learning procedure:  62
#   optimized:                             TRUE
plot(pdap5)
## https://i.imgur.com/nxCfmYf.png

## implementing all the prior knowledge cost ~30:
score(pdap5, weightC)
# [1] -254.6061724

fit <- bn.fit(pdap5, weightC); fit
#   Bayesian network parameters
#
#   Parameters of node Weight.Omron (Gaussian distribution)
#
# Conditional density: Weight.Omron | Weight.visceral.fat
# Coefficients:
#         (Intercept)  Weight.visceral.fat
#       169.181651376          2.744954128
# Standard deviation of the residuals: 1.486044472
#
#   Parameters of node Weight.BMI (Gaussian distribution)
#
# Conditional density: Weight.BMI | Weight.Omron
# Coefficients:
#   (Intercept)   Weight.Omron
# -0.3115772322   0.1411044216
# Standard deviation of the residuals: 0.03513413381
#
#   Parameters of node Weight.body.fat (Gaussian distribution)
#
# Conditional density: Weight.body.fat
# Coefficients:
# (Intercept)
# 28.70314465
# Standard deviation of the residuals: 0.644590085
#
#   Parameters of node Weight.muscle (Gaussian distribution)
#
# Conditional density: Weight.muscle | Weight.body.fat
# Coefficients:
#     (Intercept)  Weight.body.fat
#   52.1003347352    -0.6141270921
# Standard deviation of the residuals: 0.06455478599
#
#   Parameters of node Weight.resting.metabolism (Gaussian distribution)
#
# Conditional density: Weight.resting.metabolism | Weight.Omron + Weight.body.fat
# Coefficients:
#     (Intercept)     Weight.Omron  Weight.body.fat
#   666.910582196      6.767607964     -3.886694779
# Standard deviation of the residuals: 1.323176507
#
#   Parameters of node Weight.body.age (Gaussian distribution)
#
# Conditional density: Weight.body.age | Weight.Omron + Weight.body.fat
# Coefficients:
#     (Intercept)     Weight.Omron  Weight.body.fat
#  -32.2651379176     0.3603672788     0.5150134225
# Standard deviation of the residuals: 0.2914301529
#
#   Parameters of node Weight.visceral.fat (Gaussian distribution)
#
# Conditional density: Weight.visceral.fat | Weight.body.fat
# Coefficients:
#     (Intercept)  Weight.body.fat
#    6.8781100009     0.1070118125
# Standard deviation of the residuals: 0.2373649058
## residuals look fairly good, except for Weight.resting.metabolism, where there are some extreme residuals in what looks a bit like a sigmoid sort of pattern, suggesting nonlinearities in the Omron scale's formula?
bn.fit.qqplot(fit)
## https://i.imgur.com/mSallOv.png

library(lavaan)
Weight.model1 <- '
    Weight.visceral.fat ~ Weight.body.fat
    Weight.Omron ~ Weight.visceral.fat
    Weight.BMI ~ Weight.Omron
    Weight.body.age ~ Weight.Omron + Weight.body.fat
    Weight.muscle ~ Weight.body.fat
    Weight.resting.metabolism ~ Weight.Omron + Weight.body.fat
                   '
Weight.fit1 <- sem(model = Weight.model1,  data = weightC)
summary(Weight.fit1)
# lavaan (0.5-16) converged normally after 139 iterations
#
#   Number of observations                           159
#
#   Estimator                                         ML
#   Minimum Function Test Statistic               71.342
#   Degrees of freedom                                 7
#   P-value (Chi-square)                           0.000
#
# Parameter estimates:
#
#   Information                                 Expected
#   Standard Errors                             Standard
#
#                    Estimate  Std.err  Z-value  P(>|z|)
# Regressions:
#   Weight.visceral.fat ~
#     Weight.bdy.ft     0.107    0.029    3.676    0.000
#   Weight.Omron ~
#     Wght.vscrl.ft     2.745    0.477    5.759    0.000
#   Weight.BMI ~
#     Weight.Omron      0.141    0.002   82.862    0.000
#   Weight.body.age ~
#     Weight.Omron      0.357    0.014   25.162    0.000
#     Weight.bdy.ft     0.516    0.036   14.387    0.000
#   Weight.muscle ~
#     Weight.bdy.ft    -0.614    0.008  -77.591    0.000
#   Weight.resting.metabolism ~
#     Weight.Omron      6.730    0.064  104.631    0.000
#     Weight.bdy.ft    -3.860    0.162  -23.837    0.000
#
# Covariances:
#   Weight.BMI ~~
#     Weight.body.g    -0.000    0.001   -0.116    0.907
#     Weight.muscle    -0.000    0.000   -0.216    0.829
#     Wght.rstng.mt     0.005    0.004    1.453    0.146
#   Weight.body.age ~~
#     Weight.muscle     0.001    0.001    0.403    0.687
#     Wght.rstng.mt    -0.021    0.030   -0.700    0.484
#   Weight.muscle ~~
#     Wght.rstng.mt     0.007    0.007    1.003    0.316
#
# Variances:
#     Wght.vscrl.ft     0.056    0.006
#     Weight.Omron      2.181    0.245
#     Weight.BMI        0.001    0.000
#     Weight.body.g     0.083    0.009
#     Weight.muscle     0.004    0.000
#     Wght.rstng.mt     1.721    0.193

zeo <- read.csv("~/wiki/doc/zeo/gwern-zeodata.csv")
zeo$Sleep.Date <- as.Date(zeo$Sleep.Date, format="%m/%d/%Y")

## convert "05/12/2014 06:45" to "06:45"
zeo$Start.of.Night <- sapply(strsplit(as.character(zeo$Start.of.Night), " "), function(x) { x[2] })
## convert "06:45" to 24300
interval <- function(x) { if (!is.na(x)) { if (grepl(" s",x)) as.integer(sub(" s","",x))
                                           else { y <- unlist(strsplit(x, ":")); as.integer(y[[1]])*60 + as.integer(y[[2]]); }
                                         }
                          else NA
                        }
zeo$Start.of.Night <- sapply(zeo$Start.of.Night, interval)
## correct for the switch to new unencrypted firmware in March 2013;
## I don't know why the new firmware subtracts 15 hours
zeo[(zeo$Sleep.Date >= as.Date("2013-03-11")),]$Start.of.Night <- (zeo[(zeo$Sleep.Date >= as.Date("2013-03-11")),]$Start.of.Night + 900) %% (24*60)

## after midnight (24*60=1440), Start.of.Night wraps around to 0, which obscures any trends,
## so we'll map anything before 7AM to time+1440
zeo[zeo$Start.of.Night<420 & !is.na(zeo$Start.of.Night),]$Start.of.Night <- (zeo[zeo$Start.of.Night<420 & !is.na(zeo$Start.of.Night),]$Start.of.Night + (24*60))

zeoSmall <- subset(zeo, select=c(ZQ,Total.Z,Time.to.Z,Time.in.Wake,Time.in.REM,Time.in.Light,Time.in.Deep,Awakenings,Start.of.Night,Morning.Feel))
zeoClean <- na.omit(zeoSmall)
# bnlearn doesn't like the 'integer' class that most of the data-frame is in
zeoClean <- as.data.frame(sapply(zeoClean, as.numeric))

library(bnlearn)
## after a bunch of iteration, blacklisting arrows which violate the prior knowledge
bl <- data.frame(from=c("Morning.Feel", "ZQ", "ZQ", "ZQ", "ZQ", "ZQ", "ZQ", "Time.in.REM", "Time.in.Light", "Time.in.Deep", "Morning.Feel", "Awakenings", "Time.in.Light", "Morning.Feel", "Morning.Feel","Total.Z", "Time.in.Wake", "Time.to.Z", "Total.Z", "Total.Z", "Total.Z"),
                 to=c("Start.of.Night", "Total.Z", "Time.in.Wake", "Time.in.REM", "Time.in.Deep", "Morning.Feel","Start.of.Night", "Start.of.Night","Start.of.Night","Start.of.Night", "Time.to.Z", "Time.to.Z", "Time.to.Z", "Total.Z", "Time.in.Wake","Time.to.Z","Time.to.Z", "Start.of.Night", "Time.in.Deep", "Time.in.REM", "Time.in.Light"))

zeo.hc <- hc(zeoClean, blacklist=bl)
zeo.iamb         <- iamb(zeoClean, blacklist=bl)
## problem: undirected arc: Time.in.Deep/Time.in.REM; since hc inferred [Time.in.Deep|Time.in.REM], I'll copy that for iamb:
zeo.iamb <- set.arc(zeo.iamb, from = "Time.in.REM", to = "Time.in.Deep")
zeo.gs <- gs(zeoClean, blacklist=bl)
## same undirected arc:
zeo.gs <- set.arc(zeo.gs, from = "Time.in.REM", to = "Time.in.Deep")

## Bigger is better:
score(zeo.iamb, data=zeoClean)
# [1] -44776.79185
score(zeo.gs, data=zeoClean)
# [1] -44776.79185
score(zeo.hc, data=zeoClean)
# [1] -44557.6952
## hc scores best, so let's look at it:
zeo.hc
#   Bayesian network learned via Score-based methods
#
#   model:
#    [Start.of.Night][Time.to.Z|Start.of.Night][Time.in.Light|Time.to.Z:Start.of.Night]
#    [Time.in.REM|Time.in.Light:Start.of.Night][Time.in.Deep|Time.in.REM:Time.in.Light:Start.of.Night]
#    [Total.Z|Time.in.REM:Time.in.Light:Time.in.Deep][Time.in.Wake|Total.Z:Time.to.Z]
#    [Awakenings|Time.to.Z:Time.in.Wake:Time.in.REM:Time.in.Light:Start.of.Night]
#    [Morning.Feel|Total.Z:Time.to.Z:Time.in.Wake:Time.in.Light:Start.of.Night]
#    [ZQ|Total.Z:Time.in.Wake:Time.in.REM:Time.in.Deep:Awakenings]
#   nodes:                                 10
#   arcs:                                  28
#     undirected arcs:                     0
#     directed arcs:                       28
#   average markov blanket size:           7.40
#   average neighbourhood size:            5.60
#   average branching factor:              2.80
#
#   learning algorithm:                    Hill-Climbing
#   score:                                 BIC (Gauss.)
#   penalization coefficient:              3.614556939
#   tests used in the learning procedure:  281
#   optimized:                             TRUE

plot(zeo.hc)
## https://i.imgur.com/nD3LXND.png

fit <- bn.fit(zeo.hc, zeoClean); fit
#
#   Bayesian network parameters
#
#   Parameters of node ZQ (Gaussian distribution)
#
# Conditional density: ZQ | Total.Z + Time.in.Wake + Time.in.REM + Time.in.Deep + Awakenings
# Coefficients:
#    (Intercept)         Total.Z    Time.in.Wake     Time.in.REM    Time.in.Deep      Awakenings
# -0.12468522173   0.14197043518  -0.07103211437   0.07053271816   0.21121000076  -0.56476256303
# Standard deviation of the residuals: 0.3000223604
#
#   Parameters of node Total.Z (Gaussian distribution)
#
# Conditional density: Total.Z | Time.in.Wake + Start.of.Night
# Coefficients:
#    (Intercept)    Time.in.Wake  Start.of.Night
# 907.6406157850   -0.4479377278   -0.2680771514
# Standard deviation of the residuals: 68.90853885
#
#   Parameters of node Time.to.Z (Gaussian distribution)
#
# Conditional density: Time.to.Z | Start.of.Night
# Coefficients:
#    (Intercept)  Start.of.Night
# -1.02898431407   0.01568450832
# Standard deviation of the residuals: 13.51606719
#
#   Parameters of node Time.in.Wake (Gaussian distribution)
#
# Conditional density: Time.in.Wake | Time.to.Z
# Coefficients:
#   (Intercept)      Time.to.Z
# 14.7433880499   0.3289378711
# Standard deviation of the residuals: 19.0906685
#
#   Parameters of node Time.in.REM (Gaussian distribution)
#
# Conditional density: Time.in.REM | Total.Z + Start.of.Night
# Coefficients:
#      (Intercept)           Total.Z    Start.of.Night
# -120.62442964234     0.37864195651     0.06275760841
# Standard deviation of the residuals: 19.32560757
#
#   Parameters of node Time.in.Light (Gaussian distribution)
#
# Conditional density: Time.in.Light | Total.Z + Time.in.REM + Time.in.Deep
# Coefficients:
#   (Intercept)        Total.Z    Time.in.REM   Time.in.Deep
#  0.6424267863   0.9997862624  -1.0000587988  -1.0001805537
# Standard deviation of the residuals: 0.5002896274
#
#   Parameters of node Time.in.Deep (Gaussian distribution)
#
# Conditional density: Time.in.Deep | Total.Z + Time.in.REM
# Coefficients:
#   (Intercept)        Total.Z    Time.in.REM
# 15.4961459056   0.1283622577  -0.1187382535
# Standard deviation of the residuals: 11.90756843
#
#   Parameters of node Awakenings (Gaussian distribution)
#
# Conditional density: Awakenings | Time.to.Z + Time.in.Wake + Time.in.REM + Time.in.Light + Start.of.Night
# Coefficients:
#     (Intercept)        Time.to.Z     Time.in.Wake      Time.in.REM    Time.in.Light
# -18.41014329148    0.02605164827    0.05736596152    0.02291139969    0.01060661963
#  Start.of.Night
#   0.01129521977
# Standard deviation of the residuals: 2.427868657
#
#   Parameters of node Start.of.Night (Gaussian distribution)
#
# Conditional density: Start.of.Night
# Coefficients:
# (Intercept)
# 1413.382886
# Standard deviation of the residuals: 64.43144125
#
#   Parameters of node Morning.Feel (Gaussian distribution)
#
# Conditional density: Morning.Feel | Total.Z + Time.to.Z + Time.in.Wake + Time.in.Light + Start.of.Night
# Coefficients:
#     (Intercept)          Total.Z        Time.to.Z     Time.in.Wake    Time.in.Light
# -0.924662971061   0.004808652252  -0.010127269154  -0.008636841492  -0.002766602019
#  Start.of.Night
#  0.001672816480
# Standard deviation of the residuals: 0.7104115719

## some issues with big residuals at the extremes in the variables Time.in.Light, Time.in.Wake, and Time.to.Z;
## not sure how to fix those
bn.fit.qqplot(fit)
# https://i.imgur.com/fmP1ca0.png

library(lavaan)
Zeo.model1 <- '
    Time.to.Z ~ Start.of.Night
    Time.in.Wake ~ Total.Z + Time.to.Z
    Awakenings ~ Time.to.Z + Time.in.Wake + Time.in.REM + Time.in.Light + Start.of.Night
    Time.in.Light ~ Time.to.Z + Start.of.Night
    Time.in.REM ~ Time.in.Light + Start.of.Night
    Time.in.Deep ~ Time.in.REM + Time.in.Light + Start.of.Night
    Total.Z ~ Time.in.REM + Time.in.Light + Time.in.Deep
    ZQ ~ Total.Z + Time.in.Wake + Time.in.REM + Time.in.Deep + Awakenings
    Morning.Feel ~ Total.Z + Time.to.Z + Time.in.Wake + Time.in.Light + Start.of.Night
                   '
Zeo.fit1 <- sem(model = Zeo.model1,  data = zeoClean)
summary(Zeo.fit1)
# lavaan (0.5-16) converged normally after 183 iterations
#
#   Number of observations                          1379
#
#   Estimator                                         ML
#   Minimum Function Test Statistic               22.737
#   Degrees of freedom                                16
#   P-value (Chi-square)                           0.121
#
# Parameter estimates:
#
#   Information                                 Expected
#   Standard Errors                             Standard
#
#                    Estimate  Std.err  Z-value  P(>|z|)
# Regressions:
#   Time.to.Z ~
#     Start.of.Nght     0.016    0.006    2.778    0.005
#   Time.in.Wake ~
#     Total.Z          -0.026    0.007   -3.592    0.000
#     Time.to.Z         0.314    0.038    8.277    0.000
#   Awakenings ~
#     Time.to.Z         0.026    0.005    5.233    0.000
#     Time.in.Wake      0.057    0.003   16.700    0.000
#     Time.in.REM       0.023    0.002   10.107    0.000
#     Time.in.Light     0.011    0.002    6.088    0.000
#     Start.of.Nght     0.011    0.001   10.635    0.000
#   Time.in.Light ~
#     Time.to.Z        -0.348    0.085   -4.121    0.000
#     Start.of.Nght    -0.195    0.018  -10.988    0.000
#   Time.in.REM ~
#     Time.in.Light     0.358    0.018   19.695    0.000
#     Start.of.Nght     0.034    0.013    2.725    0.006
#   Time.in.Deep ~
#     Time.in.REM       0.081    0.012    6.657    0.000
#     Time.in.Light     0.034    0.009    3.713    0.000
#     Start.of.Nght    -0.017    0.006   -3.014    0.003
#   Total.Z ~
#     Time.in.REM       1.000    0.000 2115.859    0.000
#     Time.in.Light     1.000    0.000 2902.045    0.000
#     Time.in.Deep      1.000    0.001  967.322    0.000
#   ZQ ~
#     Total.Z           0.142    0.000  683.980    0.000
#     Time.in.Wake     -0.071    0.000 -155.121    0.000
#     Time.in.REM       0.071    0.000  167.090    0.000
#     Time.in.Deep      0.211    0.001  311.454    0.000
#     Awakenings       -0.565    0.003 -178.407    0.000
#   Morning.Feel ~
#     Total.Z           0.005    0.001    8.488    0.000
#     Time.to.Z        -0.010    0.001   -6.948    0.000
#     Time.in.Wake     -0.009    0.001   -8.592    0.000
#     Time.in.Light    -0.003    0.001   -2.996    0.003
#     Start.of.Nght     0.002    0.000    5.414    0.000

Genome Sequencing Costs

# https://www.genome.gov/about-genomics/fact-sheets/Sequencing-Human-Genome-cost
# http://www.genome.gov/pages/der/sequencing_costs_apr2014.xls
# converted to CSV & deleted cost per base (less precision); CSV looks like:
# https://dl.dropboxusercontent.com/u/182368464/sequencing_costs_apr2014.csv
## Date, Cost per Genome
## Sep-01,"$95,263,072"
## ...
sequencing <- read.csv("sequencing_costs_apr2014.csv")
sequencing$Cost.per.Genome <- as.integer(gsub(",", "", sub("\\$", "", as.character(sequencing$Cost.per.Genome))))
# interpret month-years as first of month:
sequencing$Date <- as.Date(paste0("01-", as.character(sequencing$Date)), format="%d-%b-%y")
head(sequencing)
##         Date Cost.per.Genome
## 1 2001-09-01        95263072
## 2 2002-03-01        70175437
## 3 2002-09-01        61448422
## 4 2003-03-01        53751684
## 5 2003-10-01        40157554
## 6 2004-01-01        28780376

l <- lm(log(Cost.per.Genome) ~ Date, data=sequencing); summary(l)
##
## Coefficients:
##                 Estimate   Std. Error  t value   Pr(>|t|)
## (Intercept) 50.969823683  1.433567932  35.5545 < 2.22e-16
## Date        -0.002689621  0.000101692 -26.4486 < 2.22e-16
##
## Residual standard error: 0.889707 on 45 degrees of freedom
## Multiple R-squared:  0.939559,   Adjusted R-squared:  0.938216
## F-statistic: 699.528 on 1 and 45 DF,  p-value: < 2.22e-16
plot(log(Cost.per.Genome) ~ Date, data=sequencing)
## https://i.imgur.com/3XK8i0h.png
# as expected: linear in log (Moore’s law) 2002–2008, sudden drop, return to Moore’s law-ish ~December 2011?
# but on the other hand, maybe the post-December 2011 behavior is a continuation of the curve
library(segmented)
# 2 break-points / 3 segments:
piecewise <- segmented(l, seg.Z=~Date, psi=list(Date=c(13970, 16071)))
summary(piecewise)
## Estimated Break-Point(s):
##             Est. St.Err
## psi1.Date 12680 1067.0
## psi2.Date 13200  279.8
##
## t value for the gap-variable(s) V:  0 0 2
##
## Meaningful coefficients of the linear terms:
##                 Estimate   Std. Error  t value   Pr(>|t|)
## (Intercept) 35.841699121  8.975628264  3.99322 0.00026387
## Date        -0.001504431  0.000738358 -2.03754 0.04808491
## U1.Date      0.000679538  0.002057940  0.33020         NA
## U2.Date     -0.002366688  0.001926528 -1.22847         NA
##
## Residual standard error: 0.733558 on 41 degrees of freedom
## Multiple R-Squared: 0.962565,  Adjusted R-squared:   0.958
with(sequencing, plot(Date, log(Cost.per.Genome), pch=16)); plot(piecewise, add=T)
## https://i.imgur.com/HSRqkJO.png
# The first two segments look fine, but the residuals are clearly bad for the third line-segment:
# it undershoots (damaging the second segment's fit), overshoots, then undershoots again. Let's try again with more breakpoints:

lots <- segmented(l, seg.Z=~Date, psi=list(Date=NA), control=seg.control(stop.if.error=FALSE, n.boot=0))
summary(segmented(l, seg.Z=~Date, psi=list(Date=as.Date(c(12310, 12500, 13600, 13750,  14140,  14680,  15010, 15220), origin = "1970-01-01", tz = "EST"))))
# delete every breakpoint below t-value of ~|2.3|, for 3 breakpoints / 4 segments:
piecewise2 <- segmented(l, seg.Z=~Date, psi=list(Date=as.Date(c("2007-08-25","2008-09-18","2010-03-12"))))
with(sequencing, plot(Date, log(Cost.per.Genome), pch=16)); plot(piecewise2, add=T)

# the additional break-point is used up on a better fit in the curve. It looks like an exponential decay/asymptote,
# so let's work on fitting that part of the graph, the post-2007 curve:
sequencingRecent <- sequencing[sequencing$Date>as.Date("2007-10-01"),]
lR <- lm(log(Cost.per.Genome) ~ Date, data=sequencingRecent); summary(lR)
piecewiseRecent <- segmented(lR, seg.Z=~Date, psi=list(Date=c(14061, 16071))); summary(piecewiseRecent)
## Estimated Break-Point(s):
##             Est. St.Err
## psi1.Date 14290  36.31
## psi2.Date 15290  48.35
##
## t value for the gap-variable(s) V:  0 0
##
## Meaningful coefficients of the linear terms:
##                 Estimate   Std. Error   t value   Pr(>|t|)
## (Intercept)  1.13831e+02  6.65609e+00  17.10182 2.0951e-13
## Date        -7.13247e-03  4.73332e-04 -15.06865 2.2121e-12
## U1.Date      4.11492e-03  4.94486e-04   8.32161         NA
## U2.Date      2.48613e-03  2.18528e-04  11.37668         NA
##
## Residual standard error: 0.136958 on 20 degrees of freedom
## Multiple R-Squared: 0.995976,  Adjusted R-squared: 0.994971

with(sequencingRecent, plot(Date, log(Cost.per.Genome), pch=16)); plot(piecewiseRecent, add=T)

lastPiece <- lm(log(Cost.per.Genome) ~ Date, data=sequencingRecent[as.Date(15290, origin = "1970-01-01", tz = "EST")<sequencingRecent$Date,]); summary(lastPiece)
## Coefficients:
##                 Estimate   Std. Error  t value   Pr(>|t|)
## (Intercept) 17.012409648  1.875482507  9.07095 1.7491e-05
## Date        -0.000531621  0.000119056 -4.46528  0.0020963
##
## Residual standard error: 0.0987207 on 8 degrees of freedom
## Multiple R-squared:  0.71366,    Adjusted R-squared:  0.677867
with(sequencingRecent[as.Date(15290, origin = "1970-01-01", tz = "EST") < sequencingRecent$Date,],
    plot(Date, log(Cost.per.Genome), pch=16)); abline(lastPiece)

predictDays <- seq(from=sequencing$Date[1], to=as.Date("2030-12-01"), by="month")
lastPiecePredict <- data.frame(Date = predictDays, Cost.per.Genome=c(sequencing$Cost.per.Genome, rep(NA, 305)), Cost.per.Genome.predicted = exp(predict(lastPiece, newdata = data.frame(Date = predictDays))))

nlmR <- nls(log(Cost.per.Genome) ~ SSasymp(as.integer(Date), Asym, r0, lrc), data=sequencingRecent); summary(nlmR)
##
## Parameters:
##          Estimate   Std. Error    t value Pr(>|t|)
## Asym  7.88908e+00  1.19616e-01   65.95328   <2e-16
## r0    1.27644e+08  1.07082e+08    1.19203   0.2454
## lrc  -6.72151e+00  5.05221e-02 -133.04110   <2e-16
##
## Residual standard error: 0.150547 on 23 degrees of freedom
with(sequencingRecent, plot(Date, log(Cost.per.Genome))); lines(sequencingRecent$Date, predict(nlmR), col=2)

# side by side:
with(sequencingRecent, plot(Date, log(Cost.per.Genome), pch=16))
plot(piecewiseRecent, add=TRUE, col=2)
lines(sequencingRecent$Date, predict(nlmR), col=3)
# as we can see, the 3-piece linear fit and the exponential decay fit identically;
# but exponential decay is more parsimonious, IMO, so I prefer that.

predictDays <- seq(from=sequencingRecent$Date[1], to=as.Date("2020-12-01"), by="month")
data.frame(Date = predictDays, Cost.per.Genome.predicted = exp(predict(nlmR, newdata = data.frame(Date = predictDays))))

https://www.unz.com/gnxp/the-intel-of-sequencing/#comment-677904 https://biomickwatson.wordpress.com/2015/03/25/the-cost-of-sequencing-is-still-going-down/

Genome sequencing historically has dropped in price ~18% per year. Consider this simple scenario: if we have a fixed amount of money to spend buying genomes, and we can afford to buy 1 genome in the first year, then the next year we can buy 1.21 genomes, then 1.48 genomes and so on and in 30 years we can afford to buy 385 genomes each year. The number we can afford in year x is:

$$y=\frac{1}{{0.82}^x}$$

sapply(0:30, function(x) 1/(0.82^x))
#  [1]   1.000000000   1.219512195   1.487209994   1.813670724   2.211793566   2.697309227   3.289401497   4.011465240   4.892030780
# [10]   5.965891196   7.275477068   8.872533010  10.820162207  13.195319764  16.091853371  19.624211428  23.931965156  29.185323361
# [19]  35.591857758  43.404704583  52.932566564  64.551910444  78.721842005  96.002246348 117.075910180 142.775500220 174.116463682
# [28] 212.337150832 258.947744917 315.789932826 385.109674178

Genomes are unlike computation, though, as they are data rather than an ephemeral service. Each genome is still useful and accumulates in a database. How many genomes total do we have each year? Quite a lot:

cumsum(sapply(0:30, function(x) 1/(0.82^x)))
#  [1]    1.000000000    2.219512195    3.706722189    5.520392914    7.732186480   10.429495707   13.718897204   17.730362444
#  [9]   22.622393224   28.588284420   35.863761488   44.736294497   55.556456704   68.751776468   84.843629839  104.467841268
# [17]  128.399806424  157.585129785  193.176987543  236.581692126  289.514258690  354.066169134  432.788011139  528.790257487
# [25]  645.866167667  788.641667886  962.758131569 1175.095282401 1434.043027318 1749.832960144 2134.942634322

While initially there’s not much of a pile to concern ourselves with, eventually we have 2000_24ya+ genomes while still only producing <400 genomes that year, a factor of 5 difference. (As it happens, if you consider UKBB at n = 500k produced as a single investment 2012^–₅2017, 23andMe in 2017 is reportedly n = 2-2.5m, so this 5x multiplier is about right.)

23andMe started back in 2007_17ya or so offering $1,435.73^$1,000₂₀₀₇ SNP panels to a few thousand people, growing to ~1m by 8 years later in July 2015. To reproduce that in this model of constant investment we start with a base of 56k SNPs purchased per year, growing according to the cost decrease:

cumsum(sapply(0:7, function(x) (56000*1)/(0.82^x)))
# [1]  56000.0000 124292.6829 207576.4426 309142.0032 433002.4429 584051.7596 768258.2434 992900.2969

What does that yield by 10 years later (2017) or 20 years later (2027)? It yields: 1.6m (1,600,943) and 16.2m (16,212,798) respectively.

Even if we assumed that annual genomes/SNPs leveled off in 2017, the linear increase pushes us into the millions range rapidly:

annualStagnation <- sapply(0:30, function(x) min(334089, (56000*1)/(0.82^x)))
cumsum(annualStagnation)
#  [1]   56000.0000  124292.6829  207576.4426  309142.0032  433002.4429  584051.7596  768258.2434  992900.2969 1266854.0206 1600943.0206
# [11] 1935032.0206 2269121.0206 2603210.0206 2937299.0206 3271388.0206 3605477.0206 3939566.0206 4273655.0206 4607744.0206 4941833.0206
# [21] 5275922.0206 5610011.0206 5944100.0206 6278189.0206 6612278.0206 6946367.0206 7280456.0206 7614545.0206 7948634.0206 8282723.0206
# [31] 8616812.0206
data.frame(Year=2007:2037, total=round(totalStagnation))
# Year   total
# 2007   56000
# 2008  124293
# 2009  207576
# 2010  309142
# 2011  433002
# 2012  584052
# 2013  768258
# 2014  992900
# 2015 1266854
# 2016 1600943
# 2017 1935032
# 2018 2269121
# 2019 2603210
# 2020 2937299
# 2021 3271388
# 2022 3605477
# 2023 3939566
# 2024 4273655
# 2025 4607744
# 2026 4941833
# 2027 5275922
# 2028 5610011
# 2029 5944100
# 2030 6278189
# 2031 6612278
# 2032 6946367
# 2033 7280456
# 2034 7614545
# 2035 7948634
# 2036 8282723
# 2037 8616812

So even if no additional funds per year start getting spent on genomics despite the increasing utility and the cost curve remains the same, the cumulative number of SNPs or whole-genomes will increase drastically over the next 30 years. Genomes on their own have many uses, such as detecting human evolution, allowing better imputation panels, inferring population structure, counting variants, detecting particularly lethal mutations etc, but of course their main use is trait prediction. Given the increases, we would expect large enough n for Hsu’s lasso to undergo phase transition and recover nearly the full SNP heritability (see point-estimates for various traits); the bottleneck increasingly will not be genomes but phenotypic measurements.

Proposal: Hand-Counting Mobile App for More Fluid Group Discussions

Groups use voting for decision-making, but existing vote systems are cumbersome. Hand-raising is faster, but does not scale because hand-counting hands is slow. Advances in machine vision may make it possible for AI to count hands in photos accurately. Combined with a smartphone’s camera, this could yield an app for fast voting in even large groups.

Medium-large (>10 people) groups face a problem in reaching consensus: ballot or pen-and-paper voting is sufficiently slow and clunky that it is too costly to use for anything but the most important discussions. A group is forced to adopt other discussion norms and save a formal vote for only the final decision, and even then the long delay kills a lot of enthusiasm and interest. Voting could be used for many more decisions if it could be faster, and of course all existing group votes would benefit from increased speed. (I am reminded of anime conventions and film festivals where, particularly for short films such as AMVs, one seems to spend more time filling out a ballot & passing them along aisles & the staff painfully counting through each ballot by hand than one actually spends watching the media in question!)

It would be better if voting could be as fluent and easy as simply raising your hand like in a small group such as a classroom—a mechanism which makes it so easy to vote that votes can be held as fast as the alternatives can be spoken aloud and a glance suffices to count (an alert group could vote on 2 or 3 topics in the time it takes to read this sentence). But hand-raising, as great as it is, suffers from the flaw that it does not scale due to the counting problem: a group of 500 people can raise their hands as easily as a group of 50 or 5, but it takes far too long to count ~250 hands: the person counting will quickly tire of the tedium, they will make mistakes counting, and this puts a serious lag on each vote, a lag which increases linearly with the number of voters. (Hands can be easy to approximate if almost everyone votes for or against something, but if consensus is so overwhelming, one doesn’t need to vote in the first place! The hard case of almost-balanced votes is the most important case.)

One might suggest using an entirely different strategy: a website with HTML polls or little clicker gizmos like used in some college lectures to administer quick quizzes. This have the downsides that they require potentially expensive equipment (I used a clicker in one class and I think it cost at least $20, so if a convention wanted to use that in an audience of hundreds, that’s a major upfront cost & my experience was that clickers were unintuitive, did not always work, and slowed things down if anything; a website would only work if you assume everyone has smartphones and is willing to pull them out to use at an instance’s notice and of course that there’s working WiFi in the room, which cannot be taken for granted) and considerable overhead in explaining to everyone how it works and getting them on the same page and making sure every person who wanders in also gets the message. (If anyone is going to be burdened with understanding or using a new system, it should be the handful of conference/festival/group organizers, not the entire audience!) A simpler approach than hands would be specially-printed paper using, for example, QR codes like piCards, which can then be recognized by standard simple computer vision techniques; this is much cheaper than clickers but still requires considerable setup & inconvenience. It’s hard to imagine a film festival running using any system, and difficult to see these systems improving on pen-and-paper ballots which at least are cheap, relatively straightforward, and well-known.

Hand-counting really does seem like the best solution, if only the counting could be fixed. Counting is something computers do fast, so that is the germ of an idea. What if a smartphone could count the votes? You don’t want a smartphone app on the entire audiences’ phones, of course, since that’s even worse than having everyone go to a website to vote; but machine vision has made enormous strides in the 2000s-2010s, reaching human-equivalent performance on challenging image recognition contests like ImageNet. (Machine vision is complicated, but the important thing is that it’s the kind of complicated which can be outsourced to someone else and turned into a dead-easy-to-use app, and the burden does not fall on the primary users—the audience.) What if the organizer had an app which took a photo of the entire audience with lifted arms and counted hands & faces and returned a vote count in a second?

Such an app would be ideal for any cultural, political, or organizational meeting. Now the flow for, eg, a film festival could go: [no explanation given to audience, one just starts] “OK, how many people liked the first short, ‘Vampire Deli’ by Ms Houston?” [everyone raises hand, smartphone flashes, 1s passes] “OK, 140 votes. How many liked the second short, ‘Cthulicious’ by Mr Iouston?” [raises hands, smartphone flashes, 1s passes] “OK… 130 people. Congratulations Ms Houston!” And so on.

Such an app might be considered an infeasible machine vision task, but I believe it could be feasible: facial localization is an old and well-studied image recognition task (and effective algorithms are built into every consumer camera), hands/fingers have very distinct shapes, and both tasks seem easier than the subtle discriminations between, say, various dog breeds demanded of ImageNet contestants.

Specifically, one could implement the machine vision core as follows:

1. multilayer neural networks trained for one task can be repurposed to similar tasks by removing the highest layer and retraining on the new task, potentially reaping great performance gains as the hybrid network has already learned much of what it needs for the second task (“transfer learning”). So one could take a publicly available NN trained for ImageNet (such as AlexNet, available in caffe), remove the top two layers, and retrain on a dataset of audiences; this will perform better since the original NN has already learned how to detect edges, recognize faces, etc

   The simpler task of counting crowds has already shown itself susceptible to deep learning: eg. “Cross-scene Crowd Counting via Deep Convolutional Neural Networks”.

2. raid Flickr and Google Images for CC-licensed photos of audiences raising their arms; then one can manually count how many arms are raised (or outsource to Amazon Mechanical Turk). With the boost from a transferred convolutional deep network, one might get good performance with just a few thousand photos to train with. If each photo takes a minute to obtain and count, then one can create a useful corpus in a week or two of work.

3. train the NN, applying the usual data augmentation tricks to increase one’s meager corpus, trying out random hyperparameters, tweaking the architecture, etc

   (Note that while NNs are very slow and computationally intensive to train, they are typically quite fast to run; the smartphone app would not be training a NN, which is indeed completely infeasible from a CPU and battery life standpoint—it is merely running the NN created by the original developer.)

4. with an accurate NN, one can wrap it in a mobile app framework. The UI, at the simplest, is simply a big button to press to take a photo, feed it into the NN, and display the count. Some additional features come to mind:

   • “headcount mode”: one may not be interested in a vote, but in how many people are in an audience (to estimate how popular a guest is, whether an event needs to move to a new bigger space, etc). If the NN can count faces and hands to estimate a vote count, it can simply report the count of faces instead.

   • the app should save every photo & count, both as an audit trail and also to support post-vote recounts in case of disputes or a desire for a more definitive count

   • the reported count should come with an indication of the NN’s uncertainty/error-rate, so users are not misled by their little handheld oracle and so they can redo a vote if the choice is borderline; Bayesian methods, in which previous votes are drawn upon, might be relevant here.

     • if the original photo could be annotated with graphical notes for each recognized/counted hand & face, this would let the user ‘see’ what the NN is thinking and would help build confidence a great deal

   • it should support manually entering in a vote-count; if the manual count differs, then this indicates the NN made an error and the photo & count should be uploaded to the original developer so it can be added to the corpus and the NN’s performance fixed in future releases of the app

   • smartphone cameras may not be high-resolution or have a sufficiently wide field-of-view to capture the entire audience at once; some sort of “montage mode” should exist so the user can swing the phone across the audience, bursts of shots taken, and the overlapping photos stitched together into a single audience photo which can be then fed into the NN as usual

   • a burst of photos might be superior to a single photo due to smartphone & hand movement blur; I don’t know if it’s best to try to combine the photos, run the NN multiple times and take the median, or feed multiple photos into the NN (perhaps by moving to a RNN architecture?)

   • the full-strength NN might still be too slow and energy-hungry to run pleasantly on a smartphone; there are model compression techniques for simplifying deep NNs to reduce the number of nodes or have fewer layers without losing much performance, which might be useful in this context (and indeed, were originally motivated by wanting to make speech-recognition run better on smartphones)

Given this breakdown, one might estimate building such an app as requiring, assuming one is already reasonably familiar with deep networks & writing mobile apps:

1. 1 week to find an ImageNet NN, learn how to modify it, and set it up to train on a fresh corpus

2. 3 weeks to create a corpus of <5000 photos with manually-labeled hand counts

3. 5 weeks to train the NN (NNs as large as ImageNet NNs take weeks to train; depending on the GPU hardware one has access to and how many tweaks and hyperparameters one tries, 5 weeks could be drastically optimistic; but on the plus side, it’s mostly waiting as the GPUs suck electricity like crazy)

4. 5 weeks to make an intuitive simple app, submitting to an app store, etc

5. These estimates are loose and probably too optimistic (although I would be surprised if it took a good developer more than 6 months to develop this app), but that would suggest >14 weeks or 784 hours of work for a developer, start to finish. (Even at minimum wage, this represents a substantial development cost of >$6k; at more plausible developer salaries, easily >$60k of salary.)

How large is the market for such an app? Groups such as anime conventions or anything on a college campus are cheapskates and would balk at a price higher than $4.99 (even if only 5 or 10 staffers need to buy it and it makes the experience much smoother). There are probably several hundred anime or video game conventions which might use this to vote, so that might be 1000 sales there. There’s easily 13,000 business conventions or conferences in the USA, which might not need voting so much, but would be attracted by a headcount mode to help optimize their event. This suggests perhaps $70k in sales with much less profit after the app store cut & taxes, much of which sales would probably be one-offs as the user reuses it for each conference. So even a wild success, in which most events adopt use of such voting software, would barely recoup the development costs; as a product, it seems this is just too much of a niche unless one could develop it much faster (such as by finding an existing corpus of hands/photos, or be certain of banging out the mobile app in much less than I estimated), find a larger market (theaters for audience participation?), or increase price substantially (10x the price and aim at only businesses?).

Air Conditioner Replacement

Is my old air conditioner inefficient enough to replace? After calculating electricity consumption for it and a new air conditioner, with discounting, and with uncertainty in parameters evaluated by a Monte Carlo method, I conclude that the savings are too small by an order of magnitude to pay for a new replacement air conditioner.

I have an old Whirlpool air conditioner (AC) in my apartment, and as part of insulating and cooling my apartment, I’ve wondered if the AC should be replaced on energy efficiency grounds. Would a new AC save more than it costs upfront? What is the optimal decision here?

Initially I was balked in analysis because I couldn’t figure out what model it was, and thus anything about it like its energy efficiency. (No model number or name appears anywhere visible on it, and I’m not going to rip it out of the wall just to look at hidden parts.)

set.seed(2015-07-26)
simulate <- function() {
    BTUs <- rnorm(1, 10000, 100)
    EER_old <- 10.7 - abs(rnorm(1, 0, 0.5)) # half-normals because efficiencies only get worse, not better
    EER_new <- 11.2 - abs(rnorm(1, 0, 0.5))
    kWh <- rnorm(1, 0.09, 0.01)
    dailyUsage <- rnorm(1, 8, 2)
    months <- sample (4:6, 1)
    minimumSavings <- rnorm(1, 47, 4)

    annualNetSavings <- ((((BTUs / EER_old ) - (BTUs / EER_new)) / 1000) * kWh * dailyUsage * 30 * months) - minimumSavings
    return(annualNetSavings)
}
sims <- replicate(100000, simulate())
summary(sims)
##        Min.     1st Qu.      Median        Mean     3rd Qu.        Max.
## -70.3666500 -46.2051500 -42.3764100 -42.1133700 -38.3134600  -0.7334517
quantile(sims, p=c(0.025, 0.975))
##        2.5%        97.5%
## -53.59989114 -29.13999204

Some Ways of Dealing With Measurement Error

Prompted by a question on LessWrong, some examples of how to analyze noisy measurements in R:

## Create a simulated dataset with known parameters, and then run a ML multilevel model, a ML SEM, &
## a Bayesian multilevel model; with the last, calculate Expected Value of Sample Information (EVSI):

## SIMULATE
set.seed(2015-08-11)
## "There is a variable X, x belongs to [0, 100]."
toplevel <- rnorm(n=1, 50, 25)
## "There are n ways of measuring it, among them A and B are widely used."
## "For any given measurer, the difference between x(A) and x(B) can be up to 20 points."
A <- toplevel + runif(1, min=-10, max=10)
B <- toplevel + runif(1, min=-10, max=10)
c(toplevel, A, B)
# [1] 63.85938385 55.43608379 59.42333264
### the true level of X we wish to recover is '63.85'

## "Between two any measurers, x(A)1 and x(A)2 can differ on average 10 points, likewise with B."
### let's imagine 10 hypothetical points are sample using method A and method B
### assume 'differ on average 10 points' here means something like 'the standard deviation is 10'
A_1 <- rnorm(n=10, mean=A, sd=10)
B_1 <- rnorm(n=10, mean=B, sd=10)

data <- rbind(data.frame(Measurement="A", Y=A_1), data.frame(Measurement="B", Y=B_1)); data
#    Measurement           Y
# 1            A 56.33870025
# 2            A 69.07267213
# 3            A 40.36889573
# 4            A 48.67289213
# 5            A 79.92622603
# 6            A 62.86919410
# 7            A 53.12953462
# 8            A 66.58894990
# 9            A 47.86296948
# 10           A 60.72416003
# 11           B 68.60203507
# 12           B 58.24702007
# 13           B 45.47895879
# 14           B 63.45308935
# 15           B 52.27724328
# 16           B 56.89783535
# 17           B 55.93598486
# 18           B 59.28162022
# 19           B 70.92341777
# 20           B 49.51360373

## MLM

## multi-level model approach:
library(lme4)
mlm <- lmer(Y ~ (1|Measurement), data=data); summary(mlm)
# Random effects:
#  Groups      Name        Variance Std.Dev.
#  Measurement (Intercept)  0.0000  0.000000
#  Residual                95.3333  9.763877
# Number of obs: 20, groups:  Measurement, 2
#
# Fixed effects:
#              Estimate Std. Error  t value
# (Intercept) 58.308250   2.183269 26.70685
confint(mlm)
#                    2.5 %       97.5 %
# .sig01       0.000000000  7.446867736
# .sigma       7.185811525 13.444112087
# (Intercept) 53.402531768 63.213970887

## So we estimate X at 58.3 but it's not inside our confidence interval with such little data. Bad luck?

## SEM

library(lavaan)
X.model <- '        X =~ B + A
                    A =~ a
                    B =~ b'
X.fit <- sem(model = X.model, meanstructure = TRUE, data = data2)
summary(X.fit)
# ...                   Estimate  Std.err  Z-value  P(>|z|)
# Latent variables:
#   X =~
#     B                 1.000
#     A              7619.504
#   A =~
#     a                 1.000
#   B =~
#     b                 1.000
#
# Intercepts:
#     a                58.555
#     b                58.061
#     X                 0.000
#     A                 0.000
#     B                 0.000
## Well, that didn't work well - explodes, unfortunately. Probably still not enough data.

## MLM (Bayesian)

library(R2jags)
## rough attempt at writing down an explicit multilevel model which
## respects the mentioned priors about errors being reasonably small:
model <- function() {
  grand.mean ~ dunif(0,100)

  delta.between.group ~ dunif(0, 10)

  sigma.between.group ~ dunif(0, 100)
  tau.between.group <- pow(sigma.between.group, -2)

  for(j in 1:K){
   # let's say the group-level differences are also normally-distributed:
   group.delta[j] ~ dnorm(delta.between.group, tau.between.group)
   # and each group also has its own standard-deviation, potentially different from the others':
   group.within.sigma[j] ~ dunif(0, 20)
   group.within.tau[j] <- pow(group.within.sigma[j], -2)

   # save the net combo for convenience & interpretability:
   group.mean[j] <- grand.mean + group.delta[j]
  }

  for (i in 1:N) {
   # each individual observation is from the grand-mean + group-offset, then normally distributed:
   Y[i] ~ dnorm(grand.mean + group.delta[Group[i]], group.within.tau[Group[i]])
  }

  }
jagsData <- list(N=nrow(data), Y=data$Y, K=length(levels(data$Measurement)),
             Group=data$Measurement)
params <- c("grand.mean","delta.between.group", "sigma.between.group", "group.delta", "group.mean",
            "group.within.sigma")
k1 <- jags(data=jagsData, parameters.to.save=params, inits=NULL, model.file=model); k1
# ...                      mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
# delta.between.group     4.971   2.945   0.221   2.353   4.967   7.594   9.791 1.008   260
# grand.mean             52.477  11.321  23.453  47.914  53.280  58.246  74.080 1.220    20
# group.delta[1]          6.017  11.391 -16.095   0.448   5.316  10.059  34.792 1.152    21
# group.delta[2]          5.662  11.318 -15.836   0.054   5.009  10.107  33.548 1.139    21
# group.mean[1]          58.494   3.765  50.973  56.188  58.459  60.838  66.072 1.001  3000
# group.mean[2]          58.139   2.857  52.687  56.366  58.098  59.851  63.999 1.003   920
# group.within.sigma[1]  12.801   2.766   8.241  10.700  12.446  14.641  18.707 1.002  1100
# group.within.sigma[2]   9.274   2.500   5.688   7.475   8.834  10.539  15.700 1.002  1600
# sigma.between.group    18.031  21.159   0.553   3.793   9.359  23.972  82.604 1.006  1700
# deviance              149.684   2.877 145.953 147.527 149.081 151.213 156.933 1.001  3000

## VOI

posteriorXs <- k1$BUGSoutput$sims.list[["grand.mean"]]
MSE <- function(x1, x2) { (x2 - x1)^2 }
lossFunction <- function(x, predictions) { mean(sapply(predictions, function(x2) { MSE(x, x2)}))}
## our hypothetical mean-squared loss if we predicted, say, X=60:
lossFunction(60, posteriorXs)
# [1] 184.7087612
## of the possible values for X, 1-100, what value of X minimizes our squared error loss?
losses <- sapply(c(1:100), function (n) { lossFunction(n, posteriorXs);})
which.min(losses)
# [1] 52
## 52 also equals the mean estimate of X, which is good since it's well known that the mean is what minimizes
## the loss when the loss is squared-error so it suggests that I have not screwed up the definitions
losses[52]
[1] 128.3478462

## to calculate EVSI, we repeatedly simulate a few hundred times the existence of a hypothetical 'C' measurement
## and draw n samples from it;
## then we add the C data to our existing A & B data; run our Bayesian multilevel model again on the bigger dataset;,
## calculate what the new loss is, and compare it to the old loss to see how much the new data
## reduced the loss/mean-squared-error.
## Done for each possible n (here, 1-30) and averaged out, this tells us how much 1 additional datapoint is worth,
## 2 additional datapoints, 3 additional datapoints, etc.
sampleValues <- NULL
for (i in seq(from=1, to=30)) {

    evsis <- replicate(500, {
        n <- i

        C <- toplevel + runif(1, min=-10, max=10)
        C_1 <- rnorm(n=n, mean=C, sd=10)
        ## all as before, more or less:
        newData <- rbind(data, data.frame(Measurement="C", Y=C_1))

        jagsData <- list(N=nrow(newData), Y=newData$Y, K=length(levels(newData$Measurement)),
                         Group=newData$Measurement)
        params <- c("grand.mean","delta.between.group", "sigma.between.group", "group.delta", "group.mean",
                    "group.within.sigma")
        jEVSI <- jags(data=jagsData, parameters.to.save=params, inits=NULL, model.file=model)

        posteriorTimesEVSI <- jEVSI$BUGSoutput$sims.list[["grand.mean"]]
        lossesEVSI <- sapply(c(1:100), function (n) { lossFunction(n, posteriorTimesEVSI);})

        oldOptimum <- 128.3478462 # losses[52]
        newOptimum <- losses[which.min(lossesEVSI)]
        EVSI <- newOptimum - oldOptimum
        return(EVSI)
        }
        )

    print(i)

    print(mean(evsis))
    sampleValues[i] <- mean(evsis)
}
sampleValues
#  [1] 13.86568780 11.07101087 14.15645538 13.05296681 11.98902668 13.86866619 13.65059093 14.05991443
#  [9] 14.80018511 16.36944874 15.47624541 15.64710237 15.74060632 14.79901214 13.36776390 15.35179426
# [17] 14.31603459 13.70914727 17.20433606 15.89925289 16.35350861 15.09886204 16.30680175 16.27032067
# [25] 16.30418553 18.84776433 17.86881713 16.65973397 17.04451609 19.17173439

## As expected, the gain in reducing MSE continues increasing as data comes in but with diminishing returns;
## this is probably because in a multilevel model like this, you aren't using the _n_ datapoints to estimate X
## directly so much as you are using them to estimate a much smaller number of latent variables, which are then
## the _n_ used to estimate X. So instead of getting hyperprecise estimates of A/B/C, you need to sample from additional
## groups D/E/F/... Trying to improve your estimate of X by measuring A/B/C many times is like trying to estimate
## IQ precisely by administering a WM test a hundred times.

## If we wanted to compare with alternatives like instead sampling n data points from C and a D, it's easy to modify
## the EVSI loop to do so: generate `D <- toplevel + runif(1, min=-10, max=10); D_1 <- rnorm(n=n, mean=D, sd=10)`
## and now `rbind` D_1 in as well. At a guess, after 5-10 samples from the current group, estimates of X will be improved more
## by then sampling from a new group.

## Or the loss function could be made more realistic. It's unlikely one is paid by MSE, and if one adds in how much
## money each sample costs, with a realistic loss function, one could decide exactly how much data is optimal to collect.

## To very precisely estimate X, when our measurements are needed to measure at least 3 latent variables,
## requires much more data than usual.

## In general, we can see the drawbacks and benefits of each approach. A canned MLM
## is very fast to write but doesn't let us include prior information or easily run
## additional analyses like how much additional samples are worth. SEM works poorly
## on small samples but is still easy to write in if we have more complicated
## models of measurement error. A full-blown modeling language like JAGS is quite
## difficult to write in and MCMC is slower than other approaches but handles small
## samples without any errors or problems and offers maximal flexibility in using
## the known prior information and then doing decision-theoretic stuff. Overall for
## this problem, I think JAGS worked out best, but possibly I wasn't using LAVAAN
## right and that's why SEM didn't seem to work well.

Value of Information: Clinical Prediction Instruments for Suicide

https://slatestarcodex.com/2015/08/31/magic-markers/

I agree. When criticizing the study for claiming the blood levels added predictive power and it’s not clear they did, this is solely a statistical claim and can be done in a vacuum. But when one then goes on to pan the predictive power of the underlying clinical prediction instruments as useless in all circumstances, based on just the prediction stats:

So when people say “We have a blood test to diagnose suicidality with 92% accuracy!”, even if it’s true, what they mean is that they have a blood test which, if it comes back positive, there’s still less than 50-50 odds the person involved is suicidal. Okay. Say you’re a psychiatrist. There’s a 48% chance your patient is going to be suicidal in the next year. What are you going to do? Commit her to the hospital? I sure hope not. Ask her some questions, make sure she’s doing okay, watch her kind of closely? You’re a psychiatrist and she’s your depressed patient, you would have been doing that anyway. This blood test is not really actionable. And then remember that this isn’t the blood test we have. We have some clinical prediction instruments that do this…But having “a blood test for suicide” won’t be very useful, even if it works.

One is implicitly making some strong cost-benefit claims here and stepping from statistics (‘what are the probabilities?’) to decision theory (‘given these probabilities, how should I act?’). They are not identical: no AUC graph will ever tell you if a model’s predictions are useful or not, and there is no universal threshold where 92% specificity/sensitivity is totally useless but 95% would make a difference—these clinical prediction instruments might be useless indeed, but that will depend on costs, base rates, and available actions. (I tried to make this point to Coyne on Twitter earlier but I don’t think he understood what I was getting at & he blew me off.)

Discontinuities come from our actions; our inferences are incremental. There are some contexts where a tiny 1% improvement in AUC might be worth a lot (Wall Street) and there are some contexts where sensitivity or specificity of 99% is still useless because it won’t change your actions at all (I’m currently comparing my riding lawn mower to a robotic lawn mower, and thus far, it doesn’t matter how precise my parameters are, the robotic lawn mowers are, to my disappointment, just too expensive right now). I think p-values have shown us how well arbitrary thresholds work out in practice (and remember where they came from in the first place! decision rules set per problem—Gosset, in optimizing a brewery, did not have the pathologies we have with p<0.05 fetishism.) I also don’t believe your choices are really that restricted: you mean if you were absolutely convinced that your patient was about to commit suicide, there is absolutely nothing you could do besides treat them like any other depressive? That seems unlikely. But whatever, even if commitment is the only alternative, there is still a value to the information provided by a clinical prediction instrument, and we can calculate it, and you should if you want to rule it out as having any value, in the same way that in criticizing a study as weak, it’s better to ignore the p-values and just work out the right posterior and demonstrate directly how little evidence it contains.

Let’s try this as an example, it’s not hard or terribly complex (just tedious). So we have a ward of 100 depressive patients where we are interested in preventing suicide; our prior probability is that 7.5% or ~7 of them will commit suicide. The value of a life has been given a lot of different valuations, but $10 million is a good starting point.

Action 1:

What are our costs or losses? We could say that we expect a loss of 7.5*$10m or −$75m, and if we stand by and do no treatment or intervention whatsoever, we spend no more money and so the total loss is

0 + 0.075 * 100 * 10,000,000 = −$75,000,000

Action 2:

Let’s say they all stay by default for one week and this costs a net $1,000 a day; let’s say further that, since commitment is the mentioned alternative, while committed a suicide attempt will fail. And since we know that suicides are so often spontaneous and major depression comes and goes, a frustrated suicide attempt doesn’t simply mean that they will immediately kill themselves as soon as they get out. This 7% comes from a followup period of a year, so the probability any will attempt suicide in the next week might be 0.075/52 or 0.001442307692. So this gives us our default setup: we have 100 patients staying for 7 days at a net cost of $1,000 a day or $700,000 total, and by having them stay, we stop an expected average of 0.14 suicides and thus we prevent an expected loss of 0.14 * $10m = $1,440,000, for a total loss of treatment-cost minus treatment-gain plus remaining-loss:

$700,000 - (0.14 * $10m) - $10m * 100 * (0.075-(0.075/52)) = −$74,257,692.

Note that this loss is smaller than in the scenario in which we don’t do any commitment at all; since one week of suicide-watch reduced the suicide loss more than it cost, this is not surprising.

Specifically, the benefit is:

action1 - action2 = gain to switching 75000000 - 74257692 = $742,308

Not fantastic, but it’s in the right order of magnitude (you can’t expect more from a low base-rate event and a treatment with such a low probability of making a difference, after all) so it looks plausible, and it’s still more than zero. We can reject the action of not committing them at all as being inferior to committing them for one week.

Action 3:

What if we were instead choosing between one week and committing them for a full year—thus catching the full 7.5% of suicides during the 1-year followup? Does that work? First, the loss from this course of action:

((100*365.2*1000) - (0 * 10000000) - (10000000 * 100 * (0.075-(0.075/1)))) = −$36,520,000

Since there are no suicides, we avoid the default loss of −$75m, but we still have to spend $36,520,000 to pay for the long-term commitment. However, the benefit to the patients has increased dramatically since we stop so many more suicides:

action 2 - action 3 = $35,637,692.31

(We go from a loss of −$74m to a loss of −$36m.) So we see action 3 is even better than action 2 for the patients. Of course, we can’t extrapolate out any further than 1 year, because that’s what our followup number is, and we don’t know how the suicide risk falls after the 1 year point—if it drops to ~0, then further commitment is a terrible idea. So I’m not going to calculate out any further. (Since this is all linear stuff, the predicted benefit will increase smoothly over the year and so there’s no point in calculating out alternatives like 1 month, 3 months, 6 months, 9 months, etc.) What’s that, action 3 is totally infeasible and no one would ever agree to this—the patients would scream their heads off and the health insurance companies would never go for it—even if we could show that long commitments do reduce the suicide rate enough to justify the costs? And, among other things, I’ve oversimplified in assuming the 7% risk is evenly distributed over the year rather than a more plausible distribution like exponentially decreasing from Day 1, so likely commitment stops being a good idea more like month 3 or something? Yeah, you’re probably right, so let’s go back to using action 2’s loss as our current best alternative.

Now, having set out some of the choices available, we can find out how much better information is worth. First, let’s ask what the Expected Value of Perfect Information is: if we were able to take our 100 patients and exactly predict which 7 were depressive and would commit suicide this year in the absence of any intervention, where our choice is between committing them for one week or not at all. Given such information we can eject the 93 who we now know were never a suicide risk, and we hold onto the 7 endangered patients, and we have a new loss of the commitment cost of 7 people for a week vs the prevented loss of the chance they will try to commit suicide that week of this year:

((771000) - (0.14 * 10000000) - (10000000 * 7 * (1-(1/52)))) = −$70,004,846

How much did we gain from our perfect information? About $4m:

74257692 - 70004846 = $4,252,846

(This passes our sanity checks: additional information should never hurt us, so the amount should be >=$0, but we are limited by the intervention to doing very little, so the ceiling should be a low amount compared to the total loss, which this is.)

So as long as the perfect information did not cost us more than $4m or so, we would have net gained from it: we would have been able to focus commitment on the patients at maximal risk. So suppose we had a perfect test which cost $1,000 a patient to run, and we wanted to know if the gained information was valuable enough to bother with using this expensive test; the answer in this case is definitely yes: with 100 patients, it’ll cost $100,000 to run the test but it’ll save $4.25m for a net profit of $4.15m. In fact, we would be willing to pay per-patient costs up to $42k, at which point we hit break-even (4252846 / 100).

OK, so that’s perfect information. What about imperfect information? Well, imperfect is a lot like perfect information, just, y’know—less so. Let’s consider this test: with the same prior, a negative on it means the patient now has P=0.007 to commit suicide that year, and a positive means P=0.48, and the sensitivity/specificity at 92%. (Just copying that from OP & ButYouDisagree, since those sound plausible.) So when we run the test on our patients, we find of the 4 possible outcomes:

• 85.1 patients are non-suicidal and the test will not flag them

• 7.4 are non-suicidal but the test will flag them

• 6.9 are suicidal and the test will flag them

• 0.6 are suicidal but the test will not flag them

So if we decide whether to commit or not commit solely based on this test, we will send home 85.1 + 0.6 = 85.7 patients (and indeed 0.6/85.7=0.007), and we will retain the remaining 7.4 + 6.9 = 14.3 patients (and indeed, 6.9/14.3=0.48). So our loss is the wrongly ejected patient of 0.6 suicides plus the cost of committing 14.3 patients (both safe and at-risk) for a week in exchange for the gain of a small chance of stopping the suicide of the 6.9 actually at risk:

(10000000*85.7*0.007) + (14.3*7*1000) + (10000000 * (0.4814.3) (1-(1/52))) = −$73,419,100

How much did we gain from our imperfect information? About $0.8m:

74257692 - 73419100 = $838,592

or $8,385.92 per patient. (This passes our sanity check: greater than $0, but much less than the perfect information. The exact amount may seem lame, but as a fraction of the value of perfect information, it’s not too bad: the test gets us 20% - 838592 / 4252846 - of the way to perfection.)

And that’s our answer: the test is not worth $0—it’s worth $8k. And once you know what the cost of administering the test is, you simply subtract it and now you have the Net Expected Value of Information for this test. (I can’t imagine it costs $8k to administer what this sounds like, so at least in this model, the value is highly likely >$0.)

By taking the posterior of the test and integrating all the estimated costs and benefits into a single framework, we can nail down exactly how much value these clinical instruments could deliver if used to guide decision-making. And if you object to some particular parameter or assumption, just build another decision-theory model and estimate the new cost. For example, maybe commitment actually costs, once you take into account all the disruption to lives and other such side-effects, not $1,000 but net of $5,000 per day, what then? Then the gain halves to $438,192, etc. And if it costs $1,0000 then the test is worth nothing because you won’t commit anyone ever because it’s just way too expensive, and now you know it’s worth $0; or if commitment is so cheap that it’s more like $100 a day, then the test is also worth $0 because you would just commit everyone (since breakeven is then a suicide probability way below 7%, all the way at ~0.4% which is still below the 0.7% which the test can deliver, so the test result doesn’t matter for deciding whether to commit, so it’s worth $0), or if you adopt a more reasonable value of life like $20m, the value of perfect information shoots up (obviously, since the avoided loss doubles) but the value of imperfect information drops like a stone (since now that one suicidal patient sent home blows away your savings from less committing) and the test becomes worthless; and playing with the formulas, you can figure out the various ranges of assumptions in which the test has positive value and estimate how much it has under particular parameters, and of course if parameters are uncertain, you can cope with that uncertainty by embedding this in a Bayesian model to get posterior distributions of particular parameters incorporating all the uncertainty.

So to sum up: there are no hard thresholds in decision-making and imposing them can cost us better decision-making, so to claim additional information is worthless, more analysis needed, and this analysis must be done with respect to the available actions & their consequences, which even under the somewhat extreme conditions here of very weak interventions & low base-rates, suggests that the value of this information is positive.

Bayesian Model Averaging

## original: "Bayesian model choice via Markov chain Monte Carlo methods" Carlin & Chib 1995 https://www.ma.imperial.ac.uk/~das01/MyWeb/SCBI/Papers/CarlinChib.pdf
## Kobe example & data from: "A tutorial on Bayes factor estimation with the product space method", Lodewyckx et al 2011 https://ejwagenmakers.com/2011/LodewyckxEtAl2011.pdf
## Lodewyckx code can be downloaded after registration & email from https://ppw.kuleuven.be/okp/software/scripts_tut_bfepsm/

## "Table 2: Observed field goals (y) and attempts (n) by Kobe Bryant during the NBA seasons of 1999 to 2006."
kobe <- read.csv(stdin(),header=TRUE)
Year, y,   n,    y.n
1999, 554, 1183, 0.47
2000, 701, 1510, 0.46
2001, 749, 1597, 0.47
2002, 868, 1924, 0.45
2003, 516, 1178, 0.44
2004, 573, 1324, 0.43
2005, 978, 2173, 0.45
2006, 399,  845, 0.47

library(runjags)
model1 <- "model{
      # 1. MODEL INDEX
      # Model index is 1 or 2.
      # Prior probabilities based on argument prior1.
      # Posterior probabilities obtained by averaging
      # over postr1 and postr2.
      M ~ dcat(p[])
      p[1] <- prior1
      p[2] <- 1-prior1
      postr1 <- 2-M
      postr2 <- 1-postr1

      # 2. MODEL LIKELIHOOD
      # For each year, successes are Binomially distributed.
      # In M1, the success rate is fixed over years.
      # In M2, the success rate is year-specific.
      for (i in 1:n.years){
       successes[i] ~ dbin(pi[M,i], attempts[i])

       pi[1,i] <- pi.fixed
       pi[2,i] <- pi.free[i]
      }

      # 3. MODEL 1 (one single rate)
      # The fixed success rate is given a Beta prior and pseudoprior.
      # Whether it is a prior or pseudoprior depends on the Model index.
      pi.fixed ~ dbeta(alpha.fixed[M],beta.fixed[M])
      alpha.fixed[1] <- alpha1.prior
      beta.fixed[1] <- beta1.prior
      alpha.fixed[2] <- alpha1.pseudo
      beta.fixed[2] <- beta1.pseudo

      # 4. MODEL 2 (multiple independent rates)
      # The year-specific success rate is given a Beta prior and pseudoprior.
      # Whether it is a prior or pseudoprior depends on the Model index.
      for (i in 1:n.years){
       pi.free[i] ~ dbeta(alpha.free[M,i],beta.free[M,i])
       alpha.free[2,i] <- alpha2.prior
       beta.free[2,i] <- beta2.prior
       alpha.free[1,i] <- alpha2.pseudo[i]
       beta.free[1,i] <- beta2.pseudo[i]
      }
      # predictive interval for hypothetical 2007 data in which Kobe makes 1000 shots:
      successes.new.1 ~ dbin(pi.fixed, 1000)
      successes.new.2 ~ dbin(pi.free[n.years], 1000)

#      success.new.weighted ~ dcat(M)
  }"
# 'prior1' value from paper
data <- list("prior1"=0.000000007451, "n.years"= length(kobe$Year), "successes"=kobe$y, "attempts"=kobe$n,
             "alpha1.prior"=1, "beta1.prior"=1, "alpha2.prior"=1, "beta2.prior"=1,
             "alpha1.pseudo"=1, "beta1.pseudo"=1, "alpha2.pseudo"=rep(1,8), "beta2.pseudo"=rep(1,8) )
# inits <- function() { list(mu=rnorm(1),sd=30,t=as.vector(apply(mailSim,1,mean))) }
params <- c("pi.free", "pi.fixed", "postr1", "postr2", "M", "successes.new.1", "successes.new.2")
j1 <- run.jags(model=model1, monitor=params, data=data, n.chains=getOption("mc.cores"), method="rjparallel", sample=500000); j1
# JAGS model summary statistics from 4000000 samples (chains = 8; adapt+burnin = 5000):
#
#                  Lower95  Median Upper95    Mean      SD Mode      MCerr MC%ofSD SSeff
# pi.free[1]        0.3145 0.46864 0.98709 0.47383 0.11553 ---0.00041958     0.4 75810
# pi.free[2]       0.10099 0.46447 0.77535 0.47005  0.1154 ---0.00042169     0.4 74887
# pi.free[3]       0.19415  0.4692 0.86566  0.4741 0.11457 ---0.00040171     0.4 81342
# pi.free[4]      0.020377 0.45146 0.69697 0.45867 0.11616 ---0.00042696     0.4 74023
# pi.free[5]      0.024472 0.43846  0.7036 0.44749 0.11757 ---0.00043352     0.4 73548
# pi.free[6]      0.076795 0.43325 0.74944 0.44318 0.11684 ---0.00043892     0.4 70863
# pi.free[7]       0.06405 0.45033 0.73614 0.45748 0.11541 ---0.00041715     0.4 76543
# pi.free[8]       0.30293 0.47267 0.97338 0.47708 0.11506 ---0.00040938     0.4 79000
# pi.fixed        0.039931 0.45756 0.97903 0.49256 0.26498 ---0.00099537     0.4 70868
# postr1                 0       0       1 0.15601 0.36287    0    0.15113    41.6     6
# postr2                 0       1       1 0.84399 0.36287    1    0.15113    41.6     6
# M                      1       2       2   1.844 0.36287    2    0.15113    41.6     6
# successes.new.1        0     463     940  492.57  265.28  454    0.99543     0.4 71019
# successes.new.2      300     473     971  477.05  116.03  473     0.4152     0.4 78094
getLogBF <- function(prior0, postr0) { log((postr0/(1-postr0)) / (prior0/(1-prior0))) }
getLogBF(0.000000007451, 0.15601)
# [1] 17.02669704
## analytic BF: 18.79; paper's MCMC estimate: 18.80; not sure where I lost 1.8 of the BF.

Dealing With All-Or-Nothing Unreliability of Data

A question was asked of me: what should one infer if one is given what would be definitive data if one could take it at face value—but one suspects this data might be totally 100% incorrect? An example would be if one wanted to know what fraction of people would answer ‘yes’ to a particular question, and one had a very small poll (n = 10) suggesting 90% say yes, but then one was also given the results from a much larger poll (n = 1000) saying 75% responded yes—but this poll was run by untrustworthy people, people that, for whatever reason, you believe might make something up half the time. You should be able to learn something from this unreliable poll, but you can’t learn everything from it because you would be burned half the time.

If not for this issue of unreliability, this would be an easy binomial problem: specify a uniform or Jeffreys prior on what percentage of people will say yes, add in the binomial data of 9⁄10, and look at the posterior. But what do we do with the unreliability joker?

## install.packages("devtools")
## devtools::install_github("rasmusab/bayesian_first_aid")
library(BayesianFirstAid)
b <- bayes.binom.test(oldData$Yes, oldData$N); b
# ...number of successes = 9, number of trials = 10
# Estimated relative frequency of success:
#   0.85
# 95% credible interval:
#   0.63 0.99
# The relative frequency of success is more than 0.5 by a probability of 0.994
# and less than 0.5 by a probability of 0.006

library(runjags)
model_string <- "model {
  x ~ dbinom(theta, n)
  theta ~ dbeta(1, 1) }"
model <- autorun.jags(model_string, monitor="theta", data=list(x=oldData$Yes, n=oldData$N)); model
# JAGS model summary statistics from 20000 samples (chains = 2; adapt+burnin = 5000):
#
#       Lower95  Median Upper95    Mean      SD Mode     MCerr MC%ofSD SSeff    AC.10   psrf
# theta 0.63669 0.85254  0.9944 0.83357 0.10329 ---0.0007304     0.7 20000 0.011014 1.0004

oldData <- data.frame(Yes=9, N=10)
simulatePoll <- function(n, pr)  { rbinom(1, size=n, p=pr); }
poll_abc <- replicate(100000, {
    # draw from our uniform prior
    p <- runif(1,min=0,max=1)
    # simulate a hypothetical poll dataset the same size as our original
    newData <- data.frame(Yes=simulatePoll(oldData$N, p), N=oldData$N)
    # were they equal? if so, save sample as part of posterior
    if (all(oldData == newData)) { return(p) }
   }
  )
resultsABC <- unlist(Filter(function(x) {!is.null(x)}, poll_abc))
summary(resultsABC)
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
# 0.3260816 0.7750520 0.8508855 0.8336383 0.9117471 0.9991691
hist(resultsABC)
# https://i.imgur.com/fn3XYQW.png

oldData2 <- data.frame(Yes=c(9,750), N=c(10,1000)); oldData2
#   Yes    N
# 1   9   10
# 2 750 1000
simulateHonestPoll <- function(n, pr)  { rbinom(1, size=n, p=pr); }
simulateFauxPoll <- function(n, pr, switchp) { if(sample(c(TRUE, FALSE), 1, prob=c(switchp, 1-switchp))) { rbinom(1, size=n, p=pr); } else { round(runif(1, min=0, max=n)); }}
poll_abc <- replicate(1000000, {
 priorp <- runif(1,min=0,max=1)
 switch <- 0.5
 n1 <- 10
 n2 <- 1000
 data1 <- data.frame(Yes=simulateHonestPoll(n1, priorp), N=n1)
 data2 <- data.frame(Yes=simulateFauxPoll(n2, priorp, switch), N=n2)
 newData <- rbind(data1, data2)
 if (all(oldData2 == newData)) { return(priorp) }
 }
)
resultsABC <- unlist(Filter(function(x) {!is.null(x)}, poll_abc))
summary(resultsABC)
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
# 0.5256471 0.7427098 0.7584650 0.7860109 0.8133581 0.9765648
hist(resultsABC)
# https://i.imgur.com/atMz0jg.png

library(runjags)
model1 <- "model {
  for (i in 1:N) {
   y[i] ~ dbinom(theta[i], n[i])
   theta[i] <- thetaOfClust[ clust[i] ]
   clust[i] ~ dcat(pi[])
  }
  pi[1]  <- switch[1]
  pi[2]  <- switch[2]
  thetaOfClust[1] ~ dbeta(1,1)
  thetaOfClust[2] ~ dunif(0,1)
 }"
j1 <- autorun.jags(model1, monitor=c("theta"), data = list(N=nrow(oldData2), y=oldData2$Yes, n=oldData2$N, switch=c(0.5, 0.5), clust=c(1,NA))); j1
# ...      Lower95  Median Upper95    Mean       SD Mode      MCerr MC%ofSD SSeff    AC.10   psrf
# theta[1] 0.70582 0.75651 0.97263 0.77926  0.07178  ---  0.001442       2  2478  0.12978 1.0011
# theta[2] 0.72446 0.75078 0.77814 0.75054 0.013646  ---0.00009649     0.7 20000 0.009458      1
plot(j1)
# https://i.imgur.com/EaqR0dD.png

switch <- 0.5
oldData3 <- data.frame(Yes=c(9,(750*switch)), N=c(10,(1000*switch)))
b2 <- bayes.binom.test(sum(oldData3$Yes), sum(oldData3$N)); b2
#
#   Bayesian First Aid binomial test
#
# data: sum(oldData3$Yes) and sum(oldData3$N)
# number of successes = 384, number of trials = 510
# Estimated relative frequency of success:
#   0.75
# 95% credible interval:
#   0.71 0.79
# The relative frequency of success is more than 0.5 by a probability of >0.999
# and less than 0.5 by a probability of <0.001

Dysgenics Power Analysis

The dysgenics hypothesis argues that due to observed reproductive patterns where the highly educated or intelligent tend to have fewer offspring, genotypic IQ (the upper bound on phenotypic IQs set by genes and the sort of thing measured by a polygenic score). If dysgenics is true, then it is an extremely important phenomenon, as important as many things that get far more attention like lead remediation; but to paraphrase Richard Hamming¹, just because a problem is important does not mean it is worth working on or researching or discussing if there is no chance of making progress—if the data is hopelessly compromised by many systematic biases which would cause false positives or if the data is too scanty to overcome random error or analyses so flexible that they could deliver any answer the partisan wishes.

Phenotypic data will, in all probability, never allow for a clear & decisive answer to the question of whether dysgenics exists or matters, as long-term comparisons are roughly as credible as noting that global piracy rates have declined while global warming increases, or paracetamol consumption rates have increased in tandem with Alzheimer’s rates; only direct examination of genetics will deliver the decisive answer. It would be nice to have an idea of how much genetic data we would need to overcome random error (and hence, whether it’s possible to make progress in the near future), which we can answer by doing some statistical power analyses.

Changes over time in genetics could be due to changes within a particular race or population (for example, in all white Englishmen), or could be due to population movements like one group replacing or migrating or merging into another (population genetics has revealed innumerable complex examples historically). The latter is possible thanks to the increasing availability of ancient DNA, often made public for researchers; so one could observe very long-term trends with cumulatively large effects (implying that small samples may suffice), but this approach has serious issues in interpretation and questions about how comparable intelligence variants may be across groups or throughout human evolution. With the former, there is less concern about interpretation due to greater temporal and ethnic homogeneity—if a GWAS on white northern Europeans in 2013_11ya turns up intelligence variants and produces a useful polygenic score, it will almost certainly work on samples of white northern Europeans in 1900_124ya too—but because the time-scale is so short the effect will be subtler and harder to detect. Nevertheless, a result within a modern population would be much more credible, so we’ll focus on that.

How subtle and hard to detect an effect are we talking about here? Woodley2012 summarizes a number of estimates:

Early in the 20^th century, negative correlations were observed between intelligence and fertility, which were taken to indicate a dysgenic fertility trend (eg. Cattell, 1936_88ya; Lentz, 1927_97ya; Maller, 1933_91ya; Sutherland, 1929_95ya). Early predictions of the rate of dysgenesis were as high as between 1 and 1.5 IQ points per decade (Cattell, 1937_87ya, 1936_88ya)…In their study of the relationship between intelligence and both completed and partially completed fertility, van Court & Bean 1985_39ya reported that the relationships were predominantly negative in cohorts born between the years 1912_112ya and 1982_42ya…Vining (1982_42ya) was the first to have attempted an estimation of the rate of genotypic IQ decline due to dysgenesis with reference to a large national probability cohort of US women aged between 24 and 34 years in 1978_46ya. He identified significant negative correlations between fertility and IQ ranging from −.104 to −.221 across categories of sex, age and race, with an estimated genotypic IQ decline of one point a generation. In a 10year follow-up study using the same cohort, Vining (1995_29ya) re-examined the relationship between IQ and fertility, now that fertility was complete, finding evidence for a genotypic IQ decline of .5 points per generation. Retherford & Sewell 1988_36ya examined the association between fertility and IQ amongst a sample of 9000 Wisconsin high-school graduates (graduated 1957_67ya). They found a selection differential that would have reduced the phenotypic IQ by .81 points per generation under the assumption of equal IQs for parents and children. With an estimate of .4 for the additive heritability of IQ, they calculated a more modest genotypic decline of approximately .33 points. The study of Ree & Earles 1991_33ya, which employed the NLSY suggests that once the differential fertility of immigrant groups is taken into consideration, the phenotypic IQ loss amongst the American population may be greater than .8 of a point per generation. Similarly, in summarizing various studies, Herrnstein & Murray (1994_30ya) suggest that “it would be nearly impossible to make the total [phenotypic IQ decline] come out to less than one point per generation. It might be twice that.” (p. 364). Loehlin (1997_27ya) found a negative relationship between the fertility of American women aged 35-44 in 1992_32ya and their educational level. By assigning IQ scores to each of six educational levels, Loehlin estimated a dysgenesis rate of .8 points in one generation. Significant contributions to the study of dysgenesis have been made by Lynn, 1996_28ya (see also: 2011_13ya) whose book Dysgenics: Genetic deterioration in modern populations provided the first estimates of the magnitude of dysgenesis in Britain over a 90 year period, putting the phenotypic loss at .069 points per year (about 1.7 points a generation assuming a generational length of 25 years). In the same study, Lynn estimated that the genotypic IQ loss was 1.64 points per generation between 1920_104ya and 1940_84ya, which reduced to .66 points between 1950_74ya and the present. Subsequent work by Lynn has investigated dysgenesis in other populations. For example Lynn (1999_25ya) found evidence for dysgenic fertility amongst those surveyed in the 1994_30ya National Opinion Research Center survey, which encompassed a representative sample of American adults, in the form of negative correlations between the intelligence of adults aged 40+ and the number of children and siblings. Lynn estimates the rate of dysgenesis amongst this cohort at .48 points per generation. In a more recent study, Lynn and van Court (2004_20ya) estimated that amongst the most recent US cohort for which fertility can be considered complete (ie. those born in the years 1940^–₉1949_75ya), IQ has declined by .9 points per generation. At the country level, Lynn & Harvey 2008_16ya have found evidence of a global dysgenesis of around .86 points between 1950_74ya and 2000_24ya, which is projected to increase to 1.28 points in the period 2000^–₅₀2050. This projection includes the assumption that 35% of the variance in cross-country IQ differences is due to the influence of genetic factors. A subsequent study by Meisenberg (2009_15ya), found that the fertility differential between developed and developing nations has the potential to reduce the phenotypic world population IQ mean by 1.34 points per decade (amounting to a genotypic decline of .47 points per decade assuming Lynn & Harvey’s 35% estimate). This assumes present rates of fertility and pre-reproductive mortality within countries. Meisenberg (2010_14ya) and Meisenberg & Kaul 2010_14ya have examined the factors through which intelligence influences reproductive outcomes. They found that amongst the NLSY79 cohort in the United States, the negative correlation between intelligence and fertility is primarily associated with g and is mediated in part by education and income, and to a lesser extent by more “liberal” gender attitudes. From this Meisenberg has suggested that in the absence of migration and with a constant environment, selection has the potential to reduce the average genotypic IQ of the US population by between .4, .8 and 1.2 points per generation.

All of these estimates are genetic selection estimates: indirect estimates inferred from IQ being a heritable trait and then treating it as a natural selection/breeding process, where a trait is selected against based on phenotype and how fast the trait decreases in each succeeding generation depends on how genetic the trait is and how harsh the selection is. So variation in these estimates (quoted estimates per generation range .3–3+) is due to sampling error, differences in populations or time periods, expressing the effect by year or generation, the estimate used for heritability, reliability of IQ estimates, and whether additional genetic effects are taken into account—for example, Woodley et al 2015 finds -.262 points per decade from selection, but in another paper argues that paternal mutation load must be affecting intelligence by ~-0.84 in the general population, giving a total of -1 per decade.

Dysgenics effects should be observable by looking at genomes & SNP data with known ages/birth-years and looking for increases in total mutations or decreases in intelligence-causing SNPs, respectively.

power.t.test(delta=0.455, power=0.8, alternative="one.sided")
#      Two-sample t test power calculation
#
#               n = 60.4155602
# ...

power.prop.test(p1=0.5, p2=0.497725, power=0.8, alternative="one.sided")
#      Two-sample comparison of proportions power calculation
#
#               n = 597,272.2524
# ...

power.t.test(delta=0.012006003, power=0.8)
#      Two-sample t test power calculation
#
#               n = 108904.194
# ...

power.t.test(delta=((43-42.742)/4.63675), power=0.8)
#      Two-sample t test power calculation
#
#               n = 5071.166739
# ...

power.t.test(n=1000, power=0.8)
#     Two-sample t test power calculation
#
#              n = 1000
#          delta = 0.1253508704
hits=437;
mean1=hits*0.5; sd1=sqrt(hits*0.5*0.5);
mean2=hits*0.497; sd2=sqrt(hits*0.497*(1-0.497));
d=(mean1-mean2)/mean(c(sd1,sd2)); d
# [1] 0.1254283986

effects <- sort(rexp(10000)/1, decreasing=TRUE)
genomeOld <- function() { ifelse(sample(c(FALSE,TRUE), prob=c(0.5, 0.5), 10000, replace = TRUE), 0, effects) }
mean(replicate(10000, sum(genomeOld())))
# [1] 5000.270218
sd(replicate(10000, sum(genomeOld())))
# [1] 69.82652816
genomeNew <- function() { ifelse(sample(c(FALSE,TRUE), prob=c(0.497, 1-0.497), 10000, replace = TRUE), 0, effects) }

mean(replicate(10000, sum(genomeOld() - genomeNew())))
# [1] 29.75354558

polyNew <- replicate(1000, sum(genomeNew()[1:4]))
polyOld <- replicate(1000, sum(genomeOld()[1:4]))
t.test(polyOld, polyNew, alternative="greater")
#   Welch Two Sample t-test
#
# data:  polyOld and polyNew
# t = 0.12808985, df = 1995.8371, p-value = 0.8980908
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#  -0.7044731204  0.8029267301
# sample estimates:
#   mean of x   mean of y
# 17.72741040 17.67818359

t.test(replicate(1000, sum(genomeOld()[1:86])), replicate(1000, sum(genomeNew()[1:86])))
#
#     Welch Two Sample t-test
#
# t = 1.2268929, df = 1997.6307, p-value = 0.2200074
# alternative hypothesis: true difference in means is not equal to 0
# 95% confidence interval:
#  -0.8642674547  3.7525210076
# sample estimates:
#   mean of x   mean of y
# 244.5471658 243.1030390

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

genomeOld <- function(efft) { ifelse(sample(c(FALSE,TRUE), prob=c(0.5, 0.5), length(efft), replace = TRUE), 0, efft) }
genomeNew <- function(efft) { ifelse(sample(c(FALSE,TRUE), prob=c(0.497, 1-0.497), length(efft), replace = TRUE), 0, efft) }

simulateStudy <- function(n, hits) {
        effects <- sort(rexp(10000)/1, decreasing=TRUE)[1:hits]
        polyOld <- replicate(n, sum(genomeOld(effects)))
        polyNew <- replicate(n, sum(genomeNew(effects)))
        t <- t.test(polyOld, polyNew, alternative="greater")
        return(data.frame(N=n, P=t$p.value, PO.mean=mean(polyOld), PO.sd=sd(polyOld), PN.mean=mean(polyNew), PN.sd=sd(polyNew))) }

hits <- 86
parallelStudies <- function(n, itr) { ldply(mclapply(1:itr, function(x) { simulateStudy(n, hits); })); }

sampleSizes <- seq(500, 5000, by=100)
iters <- 300
powerExponential <- ldply(lapply(sampleSizes, function(n) { parallelStudies(n, iters) })); summary(powerExponential)
#        N              P                  PO.mean             PO.sd             PN.mean
#  Min.   : 500   Min.   :0.000000000   Min.   :222.5525   Min.   :23.84966   Min.   :221.2894
#  1st Qu.:1600   1st Qu.:0.002991554   1st Qu.:242.8170   1st Qu.:26.46606   1st Qu.:241.3242
#  Median :2750   Median :0.023639517   Median :247.2059   Median :27.04467   Median :245.7044
#  Mean   :2750   Mean   :0.093184735   Mean   :247.3352   Mean   :27.06300   Mean   :245.8298
#  3rd Qu.:3900   3rd Qu.:0.107997575   3rd Qu.:251.7787   3rd Qu.:27.64103   3rd Qu.:250.2157
#  Max.   :5000   Max.   :0.997322161   Max.   :276.2614   Max.   :30.67000   Max.   :275.7741
#      PN.sd
#  Min.   :23.04527
#  1st Qu.:26.45508
#  Median :27.04299
#  Mean   :27.05750
#  3rd Qu.:27.63241
#  Max.   :30.85065
powerExponential$Power <- powerExponential$P<0.05
powers <- aggregate(Power ~ N, mean, data=powerExponential); powers
# 1   500 0.2133333333
# 2   600 0.2833333333
# 3   700 0.2833333333
# 4   800 0.3133333333
# 5   900 0.3033333333
# 6  1000 0.3400000000
# 7  1100 0.4066666667
# 8  1200 0.3833333333
# 9  1300 0.4133333333
# 10 1400 0.4166666667
# 11 1500 0.4700000000
# 12 1600 0.4600000000
# 13 1700 0.4666666667
# 14 1800 0.4733333333
# 15 1900 0.5233333333
# 16 2000 0.5366666667
# 17 2100 0.6000000000
# 18 2200 0.5900000000
# 19 2300 0.5600000000
# 20 2400 0.6066666667
# 21 2500 0.6066666667
# 22 2600 0.6700000000
# 23 2700 0.6566666667
# 24 2800 0.7133333333
# 25 2900 0.7200000000
# 26 3000 0.7300000000
# 27 3100 0.7300000000
# 28 3200 0.7066666667
# 29 3300 0.7433333333
# 30 3400 0.7133333333
# 31 3500 0.7233333333
# 32 3600 0.7200000000
# 33 3700 0.7766666667
# 34 3800 0.7933333333
# 35 3900 0.7700000000
# 36 4000 0.8100000000
# 37 4100 0.7766666667
# 38 4200 0.8000000000
# 39 4300 0.8333333333
# 40 4400 0.8466666667
# 41 4500 0.8700000000
# 42 4600 0.8633333333
# 43 4700 0.8166666667
# 44 4800 0.8366666667
# 45 4900 0.8666666667
# 46 5000 0.8800000000
qplot(N, Power, data=powers)  + stat_smooth()

hits <- 86
sampleSizes <- seq(8000, 30000, by=1000)
iters <- 100
genome <- function(effects) {
  t <- sample(c(1:(2015 - 1950)), 1)
  decreasedFrequency <- 0.5 - (((0.5-0.497)/(2015 - 1950)) * t)
  geneFlips <- sample(c(FALSE,TRUE), prob=c(decreasedFrequency, 1-decreasedFrequency), replace = TRUE, length(effects))
  geneValues <- ifelse(geneFlips, effects, 0)
  return(data.frame(Year=1950+t,
                    PolygenicScore=sum(geneValues)))
  }
simulateStudy <- function(n, hits) {
        effects <- sort(rexp(10000)/1, decreasing=TRUE)[1:hits]
        d <- ldply(replicate(n, genome(effects), simplify=FALSE))
        l <- lm(PolygenicScore ~ Year, data=d)
        p <- anova(l)$`Pr(>F)`[1]
        return(data.frame(N=n, P=p, PO.mean=predict(l, newdata=data.frame(Year=1950)),
                                           PN.mean=predict(l, newdata=data.frame(Year=2015)))) }
parallelStudies <- function(n, itr) { ldply(mclapply(1:itr, function(x) { simulateStudy(n, hits); })); }
powerExponentialDistributed <- ldply(lapply(sampleSizes, function(n) { parallelStudies(n, iters) })); summary(powerExponential)
powerExponentialDistributed$Power <- powerExponentialDistributed$P<0.05
powers <- aggregate(Power ~ N, mean, data=powerExponentialDistributed); powers
#        N Power
# 1   8000  0.27
# 2   9000  0.32
# 3  10000  0.35
# 4  11000  0.33
# 5  12000  0.41
# 6  13000  0.34
# 7  14000  0.41
# 8  15000  0.48
# 9  16000  0.55
# 10 17000  0.62
# 11 18000  0.55
# 12 19000  0.60
# 13 20000  0.69
# 14 21000  0.61
# 15 22000  0.65
# 16 23000  0.63
# 17 24000  0.71
# 18 25000  0.67
# 19 26000  0.71
# 20 27000  0.74
# 21 28000  0.70
# 22 29000  0.79
# 23 30000  0.83
qplot(N, Power, data=powers)  + stat_smooth()

power.t.test(d=(130/206), power=0.8)
#      Two-sample t test power calculation
#
#               n = 40.40035398
# ...

Power Analysis for Racial Admixture Studies of Continuous Variables

Admixture studies examine racial phenotypic differences in traits such as blood pressure by comparing people with ancestry from multiple groups, and correlating differences in ancestry percentage with differences in the phenotype. So, for example, African-Americans have higher blood-pressure than white Americans, and most African-Americans have an average white ancestry of something like 20-25% (see later); if having 26% white ancestry predicts slightly lower blood pressure while 24% predicts higher, that suggests the difference is (as is currently believed) genetic; and this logic can be used to narrow down to specific chromosome regions, and has contributed to study of racial differences in disease.

One application would be to thorny questions like potential group differences in non-medical traits like intelligence. The standard admixture design, requiring a few thousand subjects spanning the full range, might not necessarily work here here because of the claimed environmental effects. A proposed resolution to the question is to do an admixture study comparing African-American siblings. Siblings are highly genetically related on average (50%) but in a randomized fashion due to recombination; so two siblings, including fraternal twins, born to the same parents in the same family in the same neighborhood going to the same schools, will nevertheless have many different variants, and will differ in how related they are—the average is 50% but it could be as low as 45% or high as 55%. So given two siblings, they will differ slightly in their white ancestry, and if indeed white ancestry brings with it more intelligence variants, then the sibling with a higher whiter percentage ought to be slightly more intelligent on average, and this effect will have to be causal, as the inheritance is randomized and all other factors are equal by design. (A result using ancestry percentages measured in the general population, outside families, would be able to make far more powerful comparisons by comparing people with ~0% white ancestry to those with anywhere up to 100%, and require small sample sizes, and such analyses have been done with the expected result, but are ambiguous & totally unconvincing, as the correlation of greater whiteness with intelligence could easily be due to greater SES or greater blackness could be a marker for recent immigration or any of a number of confounds that exist.) This has historically been difficult or impossible since how does one measure the actual ancestry in siblings? But with the rise of cheap genotyping, precise measure of actual (rather than average) ancestry can be done for <$100, so that is no longer an obstacle.

588 * 0.5*(1-0.5)
# [1] 147

sqrt((588 * 0.5*(1-0.5)))
# [1] 12.12435565
12/147
# [1] 0.08163265306

L=35; sqrt(1/(16*L) - (1/(3*L^2)))
# [1] 0.03890508247

n <- 10000
siblings <- data.frame(
sibling1AncestryPercentage = rnorm(n, mean=25, sd=1),
sibling1NonancestryIQ = rnorm(n, mean=0, sd=12),

sibling2AncestryPercentage = rnorm(n, mean=25, sd=1),
sibling2NonancestryIQ = rnorm(n, mean=0, sd=12))

siblings$sibling1TotalIQ <- with(siblings, sibling1NonancestryIQ + sibling1AncestryPercentage*(0.008*15))
siblings$sibling2TotalIQ <- with(siblings, sibling2NonancestryIQ + sibling2AncestryPercentage*(0.008*15))
siblings$siblingAncestryDifference <- with(siblings, sibling1AncestryPercentage - sibling2AncestryPercentage)
siblings$siblingIQDifference <- with(siblings, sibling1TotalIQ - sibling2TotalIQ )

summary(siblings)
# ...
# siblingAncestryDifference siblingIQDifference
# Min.   :-5.370128122      Min.   :-68.2971343
# 1st Qu.:-0.932086950      1st Qu.:-11.7903864
# Median : 0.002384529      Median : -0.2501536
# Mean   : 0.007831583      Mean   : -0.4166863
# 3rd Qu.: 0.938513265      3rd Qu.: 11.0720667
# Max.   : 5.271052675      Max.   : 67.5569825
summary(lm(siblingIQDifference ~ siblingAncestryDifference, data=siblings))
# ...Coefficients:
#                             Estimate Std. Error  t value  Pr(>|t|)
# (Intercept)               -0.4192761  0.1705125 -2.45892 0.0139525
# siblingAncestryDifference  0.3306871  0.1220813  2.70874 0.0067653
#
# Residual standard error: 17.05098 on 9998 degrees of freedom
# Multiple R-squared:  0.000733338,    Adjusted R-squared:  0.0006333913
# F-statistic: 7.337294 on 1 and 9998 DF,  p-value: 0.006765343
confint(lm(siblingIQDifference ~ siblingAncestryDifference, data=siblings))
#                                    2.5 %         97.5 %
# (Intercept)               -0.75351500523 -0.08503724643
# siblingAncestryDifference  0.09138308561  0.56999105507

admixtureTest <- function(n, alpha=0.05, ancestryEffect=0.008) {
 siblings <- data.frame(
     sibling1AncestryPercentage =pmax(0, rnorm(n, mean=25, sd=1)),
     sibling1NonancestryIQ = rnorm(n, mean=0, sd=12),


     sibling2AncestryPercentage = pmax(0,rnorm(n, mean=25, sd=1)),
     sibling2NonancestryIQ = rnorm(n, mean=0, sd=12))

 siblings$sibling1TotalIQ <- with(siblings, sibling1NonancestryIQ + sibling1AncestryPercentage*(ancestryEffect*15))
 siblings$sibling2TotalIQ <- with(siblings, sibling2NonancestryIQ + sibling2AncestryPercentage*(ancestryEffect*15))
 siblings$siblingAncestryDifference <- with(siblings, sibling1AncestryPercentage - sibling2AncestryPercentage)
 siblings$siblingIQDifference <- with(siblings, sibling1TotalIQ - sibling2TotalIQ )

 p <- summary(lm(siblingIQDifference ~ siblingAncestryDifference, data=siblings))$coefficients[8]
 return(p<alpha)
}

power <- function(n, iters=10000, n.parallel=8) {
    library(parallel)
    library(plyr)
    mean(unlist(mclapply(1:n.parallel, function(i) {
        replicate(iters/n.parallel, admixtureTest(n)) }))) }

# powers <- sapply(seq(100, 10000, by=10), power)
power(100)
# [1] 0.0502
power(500)
# [1] 0.0535
power(1500)
# [1] 0.0642
power(15000)
# [1] 0.2251
power(50000)
# [1] 0.6077

admixture <- read.csv(stdin(), header=TRUE, colClasses=c("factor", "numeric", "numeric", "integer"))
Study,Mean,SD,N
"Halder et al 2008",0.143,0.133,136
"Ducci et al 2009",0.07,0.09,864
"Signorello et al 2010",0.071,0.08,379
"Bryc et al 2010",0.185,0.0465,365
"Nassir et al 2012",0.225,0.147,11712
"Bryc et al 2014",0.24,0.17,5269
"Kirkegaard et al 2016",0.17,0.11,140

# what is the standard error/precision of a population SD? http://davidmlane.com/hyperstat/A19196.html
admixture$SD.SE <- (0.71*admixture$SD) / sqrt(admixture$N)
summary(admixture)
#                    Study        Mean                 SD                  N               SD.SE
#  Bryc et al 2010      :1   Min.   :0.0700000   Min.   :0.0465000   Min.   :  136.0   Min.   :0.0009644066
#  Bryc et al 2014      :1   1st Qu.:0.1070000   1st Qu.:0.0850000   1st Qu.:  252.5   1st Qu.:0.0016954481
#  Ducci et al 2009     :1   Median :0.1700000   Median :0.1100000   Median :  379.0   Median :0.0021739221
#  Halder et al 2008    :1   Mean   :0.1577143   Mean   :0.1109286   Mean   : 2695.0   Mean   :0.0034492579
#  Kirkegaard et al 2016:1   3rd Qu.:0.2050000   3rd Qu.:0.1400000   3rd Qu.: 3066.5   3rd Qu.:0.0047591374
#  Nassir et al 2012    :1   Max.   :0.2400000   Max.   :0.1700000   Max.   :11712.0   Max.   :0.0080973057
#  Signorello et al 2010:1

library(metafor)
r.mean <- rma(yi=Mean, sei=SD/sqrt(N), measure="SMD", ni=N, data=admixture); r.mean
# Random-Effects Model (k = 7; tau^2 estimator: REML)
#
# tau^2 (estimated amount of total heterogeneity): 0.0046 (SE = 0.0027)
# tau (square root of estimated tau^2 value):      0.0680
# I^2 (total heterogeneity / total variability):   99.82%
# H^2 (total variability / sampling variability):  566.51
#
# Test for Heterogeneity:
# Q(df = 6) = 3477.2614, p-val < .0001
#
# Model Results:
#
# estimate       se     zval     pval    ci.lb    ci.ub
#   0.1578   0.0258   6.1187   <.0001   0.1072   0.2083
predict(r.mean)
#   pred     se  ci.lb  ci.ub  cr.lb  cr.ub
# 0.1578 0.0258 0.1072 0.2083 0.0153 0.3003

r.sd <- rma(yi=SD, sei=SD.SE, measure="SMD", ni=N, data=admixture); r.sd
# Random-Effects Model (k = 7; tau^2 estimator: REML)
#
# tau^2 (estimated amount of total heterogeneity): 0.0018 (SE = 0.0011)
# tau (square root of estimated tau^2 value):      0.0425
# I^2 (total heterogeneity / total variability):   99.77%
# H^2 (total variability / sampling variability):  440.67
#
# Test for Heterogeneity:
# Q(df = 6) = 3819.2793, p-val < .0001
#
# Model Results:
#
# estimate       se     zval     pval    ci.lb    ci.ub
#   0.1108   0.0162   6.8587   <.0001   0.0792   0.1425
predict(r.sd)
#    pred     se  ci.lb  ci.ub  cr.lb  cr.ub
#  0.1108 0.0162 0.0792 0.1425 0.0216 0.2001

par(mfrow=c(2,1))
forest(r.mean, slab=admixture$Study)
forest(r.sd, slab=admixture$Study)

adopteeTest <- function(n, alpha=0.05, ancestryEffect=0.008, populationAncestryMean=0.1440, populationAncestrySD=0.1008, c=0.1) {
 unrelatedSiblingSD <- sqrt(15^2 * (1-c)) # subtract 10% for same shared-environment
 siblings <- data.frame(
     sibling1AncestryPercentage = pmax(0, rnorm(n, mean=populationAncestryMean*100, sd=populationAncestrySD*100)),
     sibling1NonancestryIQ = rnorm(n, mean=0, sd=unrelatedSiblingSD),

     adopteeAncestryPercentage = pmax(0, rnorm(n, mean=populationAncestryMean*100, sd=populationAncestrySD*100)),
     adopteeNonancestryIQ = rnorm(n, mean=0, sd=unrelatedSiblingSD))

 siblings$sibling1TotalIQ <- with(siblings, sibling1NonancestryIQ + sibling1AncestryPercentage*(ancestryEffect*15))
 siblings$adopteeTotalIQ  <- with(siblings, adopteeNonancestryIQ + adopteeAncestryPercentage*(ancestryEffect*15))
 siblings$siblingAncestryDifference <- with(siblings, sibling1AncestryPercentage - adopteeAncestryPercentage)
 siblings$siblingIQDifference       <- with(siblings, sibling1TotalIQ - adopteeTotalIQ )

 p <- summary(lm(siblingIQDifference ~ siblingAncestryDifference, data=siblings))$coefficients[8]
 return(p<alpha)
}

power <- function(n, sd, iters=10000, n.parallel=8) {
    library(parallel)
    library(plyr)
    mean(unlist(mclapply(1:n.parallel, function(i) {
    replicate(iters/n.parallel, adopteeTest(n, populationAncestrySD=sd)) }))) }

ns <- seq(100, 10000, by=100)
powerLow  <- sapply(ns, function(n) { power(n, sd=0.0216)})
powerMean <- sapply(ns, function(n) { power(n, sd=0.1108)})
powerHigh <- sapply(ns, function(n) { power(n, sd=0.2001)})

library(ggplot2); library(gridExtra)
pl <- qplot(ns, powerLow)  + coord_cartesian(ylim = c(0,1))
pm <- qplot(ns, powerMean) + coord_cartesian(ylim = c(0,1))
ph <- qplot(ns, powerHigh) + coord_cartesian(ylim = c(0,1))
grid.arrange(pl, pm, ph, ncol=1)

Operating on an Aneurysm

In the excellent neurosurgery memoir Do No Harm: Stories of Life, Death, and Brain Surgery (Henry Marsh 2014), chapter 2 “Aneurysm”, there is a passage on weighing the costs of action and inaction:

“A thirty-two-year-old woman,” he said tersely. “For surgery today. Had some headaches and had a brain scan.” As he talked a brain scan flashed up on the wall.

…“It’s an unruptured aneurysm, seven millimetres in size,” Fiona—the most experienced of the registrars—said. “So there’s a point zero five per cent risk of rupture per year according to the international study published in 1998_26ya.” “And if it ruptures?” “Fifteen per cent of people die immediately and another thirty per cent die within the next few weeks, usually from a further bleed and then there’s a compound interest rate of four per cent per year.”

…If we did nothing the patient might eventually suffer a haemorrhage which would probably cause a catastrophic stroke or kill her. But then she might die years away from something else without the aneurysm ever having burst. She was perfectly well at the moment, the headaches for which she had had the scan were irrelevant and had got better. The aneurysm had been discovered by chance. If I operated I could cause a stroke and wreck her—the risk of that would probably be about four or five per cent. So the acute risk of operating was roughly similar to the life-time risk of doing nothing. Yet if we did nothing she would have to live with the knowledge that the aneurysm was sitting there in her brain and might kill her any moment.

Reading this, I was a little surprised by Marsh’s evaluation given those specific numbers. Intuitively, it did not seem to me that a single risk of 5% was anywhere near as bad as a lifelong risk of 0.5%, for a 32 year old woman who would probably live another 50 years—the one number is 10x bigger than the other, but the other number is 50x bigger, and a quick heuristic for the total probability of many independent small probabilities is to just sum them up, suggesting that the risk of the untreated aneurysm was much worse (50*0.005=0.25, and 0.25>0.05). So I thought after I finished reading the book, I would work it out a little more accurately.

MCTS

An implementation in R of a simple Monte Carlo tree search algorithm (using Thompson sampling rather than a UCT) implemented with data.tree. This MCTS assumes binary win/loss (1/0) terminal rewards with no intermediate rewards/costs so it cannot be used to solve general MDPs, and does not expand leaf nodes in the move tree passed to it. (I also suspect parts of it are implemented wrong though it reaches the right answer in a simple Blockworld problem and seems OK in a Tic-Tac-Toe problem. I have since understood how memoization works with backwards induction and would probably probably drop the painful explicit tree manipulation in favor of the indirect recursive aproach.)

library(data.tree)
## MCTS helper functions:
playOutMoves <- function(move, state, actions) {
  for (i in 1:length(actions)) {
     state <- move(state, actions[i])$State
    }
    return(state)
    }
playOutRandom <- function(move, state, actions, timeout=1000, verbose=FALSE) {
 action <- sample(actions, 1)
 turn <- move(state, action)
 if(verbose) { print(turn); };
 if (turn$End || timeout==0) { return(turn$Reward) } else {
                               playOutRandom(move, turn$State, actions, timeout=timeout-1, verbose) }
 }

createTree <- function(plys, move, moves, initialState, tree=NULL) {
 if (is.null(tree)) { tree <- Node$new("MCTS", win=0, loss=0) }
 if (plys != 0) {
  for(i in 1:length(moves)) {
    x <- tree$AddChild(moves[i], win=0, loss=0)
    createTree(plys-1, move, moves, initialState, tree=x)
  }
 }
 # cache the state at each leaf node so we don't have to recompute each move as we later walk the tree to do a rollout
 tree$Do(function(node) { p <- node$path; node$state <- playOutMoves(move, initialState, p[2:length(p)]); }, filterFun = isLeaf)
 return(tree)
}

mcts <- function (tree, randomSimulation, rollouts=1000) {
 replicate(rollouts, {
   # Update posterior sample for each node based on current statistics and use Thompson sampling.
   # With a beta uniform prior (Beta(1,1)), update on binomial (win/loss) is conjugate with simple closed form posterior: Beta(1+win, 1+n-win).
   # So we sample directly from that posterior distribution for Thompson sampling
   tree$Do(function(node) { node$Thompson <- rbeta(1, 1+node$win, 1+(node$win+node$loss)-node$win) })
   # find & run 1 sample:
   node <- treeWalk(tree)
   rollout <- randomSimulation(node$state)
   if(rollout==1) { node$win <- node$win+1; } else { node$loss <- node$loss+1; }

   # propagate the new leaf results back up tree towards root:
   tree$Do(function(x) { x$win  <- Aggregate(x, "win",  sum); x$loss <- Aggregate(x, "loss", sum) }, traversal = "post-order")
  })
}

## walk the game tree by picking the branch with highest Thompson sample down to the leaves
## and return the leaf for a rollout
treeWalk <- function(node) {
    if(length(node$children)==0) { return(node); } else {
        children <- node$children
         best <- which.max(sapply(children, function(n) { n$Thompson; } ))
        treeWalk(children[[best]]) } }

mctsDisplayTree <- function(tree) {
    tree$Do(function(node) { node$P <- node$win / (node$win + node$loss) } )
    tree$Sort("P", decreasing=TRUE)
    print(tree, "win", "loss", "P", "Thompson")
    }

## Blockworld simulation
## 0=empty space, 1=agent, 2=block, 3=goal point
blockActions <- c("up", "down", "left", "right")
blockInitialState <- matrix(ncol=5, nrow=5, byrow=TRUE,
                       data=c(0,0,0,0,1,
                              0,2,0,0,2,
                              0,0,0,2,0,
                              0,2,0,0,0,
                              0,0,0,0,3))
blockMove <- function(state, direction) {
   if(state[5,5] == 2) { return(list(State=state, Reward=1, End=TRUE)) }
   position <- which(state == 1, arr.ind=TRUE)
   row <- position[1]; col <- position[2]
   rowNew <- 0; colNew <- 0
   switch(direction,
     # if we are at an edge, no change
     up    = if(row == 1) { rowNew<-row; colNew<-col; } else { rowNew <- row-1; colNew <- col; },
     down  = if(row == 5) { rowNew<-row; colNew<-col; } else { rowNew <- row+1; colNew <- col; },
     left  = if(col == 1) { rowNew<-row; colNew<-col; } else { rowNew <- row;   colNew <- col-1; },
     right = if(col == 5) { rowNew<-row; colNew<-col; } else { rowNew <- row;   colNew <- col+1; }
   )
   # if there is not a block at the new position, make the move
   if (state[rowNew,colNew] != 2) {
      state[row,col] <- 0
      state[rowNew,colNew] <- 1
      return(list(State=state, Reward=0, End=FALSE))
       } else {
               state[rowNew,colNew] <- 1
               state[row,col] <- 0
               switch(direction,
                # if the block is at the edge it can't move
                up    = if(rowNew == 1) { } else { state[rowNew-1,colNew] <- 2 },
                down  = if(rowNew == 5) { } else { state[rowNew+1,colNew] <- 2 },
                left  = if(colNew == 1) { } else { state[rowNew,colNew-1] <- 2 },
                right = if(colNew == 5) { } else { state[rowNew,colNew+1] <- 2 } )
                # a block on the magic 5,5 point means a reward and reset of the playing field
                if(state[5,5] == 2) { return(list(State=state, Reward=1, End=TRUE)) } else { return(list(State=state, Reward=0, End=FALSE)) }
                }
}

## Blockworld examples:
# blockMove(blockInitialState, "left")
# blockMove(blockInitialState, "down")
# blockMove(blockInitialState, "right")$State
# blockMove(blockMove(blockInitialState, "right")$State, "down")
# blockMove(blockMove(blockMove(blockInitialState, "down")$State, "down")$State, "down")
# playOutMoves(blockMove, blockInitialState, c("down", "down", "down"))
# playOutRandom(blockMove, blockInitialState, blockActions)

tree <- createTree(2, blockMove, blockActions, blockInitialState)
mcts(tree, function(state) { playOutRandom(blockMove, state, blockActions) })
mctsDisplayTree(tree)

tree2 <- createTree(3, blockMove, blockActions, blockInitialState)
mcts(tree2, function(state) { playOutRandom(blockMove, state, blockActions) })
mctsDisplayTree(tree2)

## Tic-Tac-Toe
tttActions <- 1:9
tttInitialState <- matrix(ncol=3, nrow=3, byrow=TRUE, data=0)
tttMove <- function(state, move) {
   move <- as.integer(move)
   # whose move is this? Player 1 moves first, so if the number of pieces are equal, it must be 1's turn:
   player <- 0;  if(sum(state == 1) == sum(state == 2)) { player <- 1 } else { player <- 2}

   # check move is valid:
   if(state[move] == 0) { state[move] <- player }

   ## enumerate all possible end-states (rows, columns, diagonals): victory, or the board is full and it's a tie
   victory <- any(c(
       all(state[,1] == player),
       all(state[1,] == player),
       all(state[,2] == player),
       all(state[2,] == player),
       all(state[,3] == player),
       all(state[3,] == player),
       all(as.logical(c(state[1,1], state[2,2], state[3,3]) == player)),
       all(as.logical(c(state[1,3], state[2,3], state[3,1]) == player))
   ))
   tie <- all(state != 0)

   # if someone has won and the winner is player 1, then a reward of 1
   if(victory) { return(list(State=state, Reward=as.integer(player==1), End=TRUE)) } else {
    if(tie) { return(list(State=state, Reward=0, End=TRUE)) } else {
      return(list(State=state, Reward=0, End=FALSE)) }
     }
}

## Tic-Tac-Toe examples:
# tttMove(tttMove(tttMove(tttInitialState, 5)$State, 9)$State, 2)
# playOutMoves(tttMove, tttInitialState, c(5, 9, 2))
# playOutRandom(tttMove, tttInitialState, tttActions, verbose=TRUE)

treeTTT <- createTree(2, tttMove, tttActions, tttInitialState)
mcts(treeTTT, function(state) { playOutRandom(tttMove, state, tttActions) })
mctsDisplayTree(treeTTT)
## hypothetical: if opponent plays center (5), what should be the reply?
treeTTT2 <- createTree(2, tttMove, tttActions, tttMove(tttInitialState, 5)$State)
mcts(treeTTT2, function(state) { playOutRandom(tttMove, state, tttActions) })
mctsDisplayTree(treeTTT2)

DeFries-Fulker Power Analysis

DeFries-Fulker (DF) extremes analysis

• DeFries & Fulker1985

• DeFries et al 1987_37ya, “Evidence for a genetic aetiology in reading disability of twins” /doc/genetics/heritable/1987-defries-2.pdf

• DeFries & Fulker 1988_36ya, “Multiple regression analysis of twin data: Etiology of deviant scores versus individual differences”

• “A Model-Fitting Implementation of the DeFries-Fulker Model for Selected Twin Data” Purcell & Sham 2003_21ya /doc/genetics/heritable/2003-purcell.pdf

• LaBuda et al 1986

• DeFries et al 1991_33ya, “Colorado Reading Project: An update”

• Gillespie & Neale 2006_18ya “A Finite Mixture Model for Genotype and Environment Interactions: Detecting Latent Population Heterogeneity” https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.534.6298&rep=rep1&type=pdf

• Purcell et al 2001_23ya, “Comorbidity between verbal and non-verbal cognitive delays in 2-year-olds: A bivariate twin analysis”

• https://www.sciencedirect.com/science/article/pii/S0160289614001676 “Thinking positively: The genetics of high intelligence”, Shakeshaft et al 2015

• Schwartz & Connolly2017, “An Investigation of Genetic and Environmental Influences Across The Distribution of Self-Control”

mixture model

generateSiblingPair <- function(ID=TRUE) {
   ## Population mean 100, SD 15; let's make family means distributed normally too;
   ## heritability 0.8, shared environment 0.1, siblings share half of genes on average + shared environment
   ## so a pair of siblings has 1 - (0.8*0.5+0.1) = 0.5 of the variance of the general population.
   parental <- mean(rnorm(1,mean=100,sd=15*0.8), rnorm(1,mean=100,sd=15*0.8))
   siblings <- rnorm(2, mean=parental, sd=15*(1 - (0.8*0.5+0.1)))
   ## Siblings will tend to vary this much, unless they are, lamentably, one of the, say,
   ## 5% struck by mutational lightning and reduced to an IQ of, let's say, 80
   if(ID) { siblings <- ifelse(rbinom(2,1,prob=0.05), siblings,rnorm(2, mean=80, sd=15)) }
   return(c(max(siblings), min(siblings)))
}
generateSiblingPairs <- function(n,ID=TRUE) { as.data.frame(t(replicate(n, generateSiblingPair(ID=ID)))) }
## dataset with lightning:
df <- round(rescale(generateSiblingPairs(1000000, ID=TRUE), mean=5, sd=2))
## floor/ceiling at 0/9 for everyone:
df[df$V1>9,]$V1 <- 9
df[df$V1<1,]$V1 <- 1
df[df$V2>9,]$V2 <- 9
df[df$V2<1,]$V2 <- 1

## dataset without:
df2 <- round(rescale(generateSiblingPairs(1000000, ID=FALSE), mean=5, sd=2))
df2[df2$V1>9,]$V1 <- 9
df2[df2$V1<1,]$V1 <- 1
df2[df2$V2>9,]$V2 <- 9
df2[df2$V2<1,]$V2 <- 1

par(mfrow=c(2,1))
hist(df$V1 - df$V2)
hist(df2$V1 - df2$V2)

## mixture modeling:
library(flexmix)
## check k=1 vs k=2 on df1, where k=2 is ground truth:
g1.1 <- flexmix(I(V1-V2) ~ 1, k=1, data=df)
g1.2 <- flexmix(I(V1-V2) ~ 1, k=2, data=df)
summary(g1.1); summary(g1.2)

## check k=1 vs k=2 on df2, where k=1 is ground truth:
g2.1 <- flexmix(I(V1-V2) ~ 1, k=1, data=df2)
g2.2 <- flexmix(I(V1-V2) ~ 1, k=2, data=df2)
summary(g2.1); summary(g2.2)

Inferring Mean IQs From SMPY/TIP Elite Samples

The properties of statistical distributions can be very different from the properties of specific subsets of those distributions in counterintuitive ways. A point drawn from an extreme will exhibit “regression to the mean”, a phenomenon which routinely trips people up. Another common example is that a small difference in means for many distributions can lead to large differences in extreme subsets.

For example, male and female average heights differ by a relatively small amount, inches at most. So in a random sample, plenty of women will be taller than men, and vice versa. However, if instead ask the sex of the tallest person in the sample, it will often be male, and the larger the sample, the more certain we can be that it will be male, and that the top X% by height will be male. Likewise, if we wanted to start a basketball league and recruited the tallest 100 people in the country, this small mean difference will show up as our entire basketball league turning out to be male. (And since height is highly heritable, we may find out that many of them are related!) What seemed like a small difference become a large one; we could have worked it out in advance if we had thought about it.

Reasoning from the general to the particular turned out to be tricky in this case because we were dealing with extreme values rather than random samples—1 basketball player chosen by height from thousands of people. Many things of great interest turn out to be like that: we are interested in the extremes much more than the expectation. Running a 2-hour marathon is an extreme on athleticism; winning the Nobel is an extreme on scientific accomplishment; being enlisted in the NBA is an extreme on height; being admitted to MIT/Stanford/Harvard is an extreme on intelligence; murdering someone is an extreme on violence; winning an Academy Award is an extreme on acting success. When we ask questions like, “why does the world record in this sport keep being shattered” or “why are so many NBA players related” or “how good can we expect the best chess player to be in 10 years” or “does this racial composition prove bias” or “how much more important are the best authors in literature than obscurer figures” or “why do so few women win the Field Medal”, we’re asking extreme value questions whose answers may be counterintuitive—and the answer may be as simple as the shape of distributions, and a slightly lower mean here or a slightly higher standard deviation there. (Working backwards from a sample selected for passing a threshold to a mean can be called “the method of limits” or “the method of thresholds”.)

The study “When Lightning Strikes Twice: Profoundly Gifted, Profoundly Accomplished”, Makel et al 2016 describes the accomplishments of the Duke TIP sample, 259 children selected for their intelligence by taking the highest-scorers out of 425,000 adolescents taking the SAT (usually <13yo) starting in 1981_43ya, representing the top 0.01% of the test-takers. The TIP sample parallels the better-known SMPY sample, which also selected extremely intelligent adolescents, who were included in a longitudinal sample. It’s frequently suggested, based on anecdotal evidence or some biased convenience samples, that more intelligence may not be better; extremely intelligent people may be unhealthy, neurotic, insane, isolated, lonely, discriminated against by society and their peers, and doomed to failure; or if things are not quite that dire, as all studies show things improving up to 130, then at around that point greater intelligence may stop making any difference, and there be little difference between someone with an IQ of 130 and 160. This is difficult to study cross-sectionally, because once you start talking about as extreme as 0.01%, it is difficult to recruit any subjects at all, and your sample will be biased in unknown ways; if you only look at successful people, you are missing the hypothetical homeless bum living out of a trash can who is a troubled and misunderstood genius. To solve these problems, you want to filter through hundreds of thousands of people so you can select the very brightest possible, and you want to find them as early as possible in life, before they have had any chance to fail or succeed, and track them longitudinally as they grow up. This is what the SMPY & TIP studies do, and the results are that the subjects are spectacularly successful in life; great intelligence is not harmful and the returns to greater intelligence are not zero even as high as 1 in 10,000.

Makel et al 2016 also reports the ethnic breakdown of the TIP and SMPY samples: 72% white, 22% Asian, 6% not reported or other. This distribution might seem remarkable given that subjects taking the SAT in 1981_43ya were born ~1970, when the USA was ~77% white, ~11% black, and ~0.7% Asian, so white are slightly under-represented, blacks are very under-represented (even if we assume all 6% are black, then that’s still half), and Asians are 31x (!) overrepresented.

## SMPY/TIP sample size & ethnic percentages: https://pbs.twimg.com/media/Cj9DXwxWEAEaQYk.jpg
tip <- 259; smpy <- 320 ## total: 579
white <- ((0.65*tip) + (0.78*smpy)) / (tip+smpy)
asian <- ((0.24*tip) + (0.20*smpy)) / (tip+smpy)
white; asian
# [1] 0.7218480138
# [1] 0.2178929188

# https://drjamesthompson.blogspot.com/2016/06/some-characteristics-of-eminent-persons.html
# > The data on ethnicity are rather sparse, but we can do a little bit of work on them by looking at US Census
# > figures for the 1970s when most of these children were born: White 178,119,221...Asia 1,526,401...So, in the
# > absence of more detailed particulars about the Other category, Asians win the race by a country mile. If we
# > simplify things by considering only Whites, Blacks and Asians the US in 1970 then the country at that time was
# > 88% White, 11% Black, and less than 1% Asian. The actual results of eminent students are 77% White, 0% Black,
# > 22% Asian. No need for a Chi square.
#
# Asian is 0.7%: 1526401  / (178119221 / 0.80)
whiteRR <- white / 0.77; asianRR <- asian / 0.007
whiteRR; asianRR
# [1] 0.937464953
# [1] 31.12755983

Of course, races in the USA have long differed by mean intelligence, with the rule of thumb being Asians ~105 IQ, whites ~100, and blacks ~90. So the order is expected—but still, 31x! Are the results being driven by some sort of pro-Asian bias or otherwise bizarre?

But this is an extreme sample. 1-in-10,000 is far out on the tails: 3.71SDs.

-qnorm(1/10000)
# [1] 3.719016485

Maybe this is normal. Can we work backwards from the overrepresentations to what differences would have generated them?

Yes, we can, even with this small sample which is so extreme and unrepresentative of the general population. This is because it is an order statistics problem: we know the order represented by the sample and so can work back to parameters of the distribution the order statistics are being generated by. Since IQ is a normal distribution, we know the overrepresentation RR, and the exact cutoff/limit used in the sample, we can convert the limit to a standard deviations, and then find the normal distribution $N(100+x,15)$ which is RR (31) times the normal distribution $N(100,15)$ at a standard deviations.

We can compare using two pnorms and shifting the second by a SDs. So for example, shifting by 15 IQ points or 1 SD would lead to 84x overrepresentation

pnorm(qnorm(1/10000)) / pnorm(qnorm(1/10000) - (15/15))
# [1] 84.39259519

We would like to solve for the shift which leads to an exact overrepresentation like 31.127; an optimization routine like R’s optim function can do that, but it requires an error to minimize, so minimizing pnorm()/pnorm(x) doesn’t work since it just leads to negative infinity, nor will RR == pnorm()/pnorm(x) work, because it evaluates to 0 for all values of x except the exact right one . Instead, we minimize the squared error between the ratio predicted by a particular x and our observed RR. This works:

## An optimization routine which automatically finds for us the IQ increase which most closely matches the RR:
solver <- function(RR, cutoff=10000) {
    optim(1,
        function(IQ_gain) { (RR - (pnorm(qnorm(1/cutoff)) / pnorm(qnorm(1/cutoff)-(IQ_gain/15))))^2 },
        )$par }

100 + solver(whiteRR)
# [1] 99.75488281
100 + solver(asianRR)
# [1] 111.8929688

So our inferred white & Asian populations means are: 99.8 and 111.9. These are relatively close to the expected values.

This approach can be used to infer other things as well. For example, the SMPY/TIP papers have not, as far as I’ve seen, mentioned what fraction of the white subjects were ethnic Jewish; since they are so over-represented in areas like Nobel prizes, we would expect many of the SMPY/TIP white students to have been Jewish. Using an estimate of the Jewish population in 1970_54ya and estimates of their mean IQ, we can work forward to what fraction of SMPY/TIP subjects might be Jewish. The 1970–1971_53ya National Jewish Population Study estimated “5,800,000 persons (of whom 5,370,000 were Jews) living in Jewish households” out of a total US population of 205 million, or 2.8% of the total population or ~3.6% of the white population. So of the ~418 white subjects, ~15 would be expected to be Jewish under the null hypothesis of no difference. The majority of American Jews are of Ashkenazi descent³, for whom intelligence estimates are debated but tend to range 105-115 (with occasional samples suggesting even higher values, like Levinson1957). In the Barbe1964 Ohio sample (IQ ~143), 8% were Jewish⁴; in Terman’s (ratio IQ >140) 1920s sample in SF/LA, 10% were Jewish; Hollingworth’s1930s sample (>180) turned up 51⁄55 or 90% Jewish⁵; Byrns1936’s 1931_93ya Wisconsin state sample found 18% of the Jewish sample to be in the top decile vs 10% American; in the Hunter College Elementary School sample 1948^–₁₂1960_64ya (>140, mean 157) in New York City, 62% were Jewish (Subotnik et al 1989, Subotnik et al 1993⁶). Given estimates of the Jewish population of children in those specific times and places, one could work backwards to estimate a Jewish mean.

We can calculate the fraction of the white sample being Jewish for each possible mean IQ:

proportion <- function (gain, cutoff=10000) {
   (pnorm(qnorm(1/cutoff)) / pnorm(qnorm(1/cutoff)-(gain/15))) }
possibleIQs <- seq(5, 15, by=0.5)
data.frame(Advantage=possibleIQs, Fraction.of.white=(sapply(possibleIQs, proportion) * 15) / 418)
   Advantage Fraction.of.white
1        5.0      0.1415427303
2        5.5      0.1633099334
3        6.0      0.1886246225
4        6.5      0.2180947374
5        7.0      0.2524371552
6        7.5      0.2924980125
7        8.0      0.3392769622
8        8.5      0.3939561508
9        9.0      0.4579348680
10       9.5      0.5328710150
11      10.0      0.6207307813
12      10.5      0.7238482059
13      11.0      0.8449966589
14      11.5      0.9874747049
15      12.0      1.1552093388
16      12.5      1.3528802227
17      13.0      1.5860693342
18      13.5      1.8614413902
19      14.0      2.1869615788
20      14.5      2.5721585555
21      15.0      3.0284424112

Judging from earlier samples with high cutoffs, I’d guess SMPY/TIP has at least a majority Jewish, giving a mean IQ of ~110; this is similar to estimates based on regular samples & estimation. This result is also similar to La Griffe du Lion’s 2003_21ya threshold analysis estimating a mean IQ of 112 based on Ashkenazi overrepresentation among USSR championship chess players, 111 based on Western Fields Medal awards, and 110 based on the USA/Canada Putnam competition. But if the mean IQ was as high as 112, then almost every single white subject would be Jewish in every sampling, which seems implausible and like something so striking that anyone writing or involved with SMPY/TIP would have to have mentioned at some point—right?

For the same reason, the original estimate of 112 for the Asians strikes me as on the high side. This could be due to problems in the data like underestimating the Asian population at the time—perhaps the Southeast/Midwest states that TIP samples from were more than 0.7% Asian—or it could be due to sampling error (only n = 579, after all).

Working backwards doesn’t immediately provide any measurement of precision or confidence intervals. Presumably someone has worked out analytic formulas which come with standard errors and confidence intervals, but I don’t know it. Instead, since the selection process which generated our data is straightforward (population mean → millions of samples → take top 1-in-10000s → calculate overrepresentation), I can again use Approximate Bayesian computation (ABC) to turn a simulation of the data generating process into a method of Bayesian inference on the unknown parameters (population means) and get credible intervals.

What sort of confidence do we have in these estimates given that these RRs are based only on? We can simulate SMPY/TIP-like selection by taking the hypothetical means of the two groups, generating ~3 million simulates (579 * 10000) each, selecting the top 1⁄10000th⁷, taking the RRs and then solving for the mean IQ. If we provide a prior on the means and we hold onto only the means which successfully generate SMPY/TIP-like fractions of 72% & 21%, this becomes ABC with the saved means forming the posterior distribution of means. (It would likely be faster to use MCMC like JAGS, but while JAGS provides truncated normal distributions which one could sample from quickly, and the necessary pnorm/qnorm functions, but it’s not clear to me how one could go about estimating the overperformance ratio and the binomial.⁸) For my priors, I believe that the rule of thumbs of 100⁄105 are accurate and highly unlikely to be more than a few points off, so I use a very weak prior of populations means being $N(100/105,4)$.

In exact ABC, we would keep only data which exactly matched 72%/22%, but that would require rejecting an extremely large number of samples. Here we’ll loosen it to ±2% tolerance:

simulateTIPSMPY <- function() {
    ## informative priors: IQs are somewhere close to where we would estimate based on other datasets
    whiteMean <- round(rnorm(1, mean=100, sd=5), digits=2)
    asianMean <- round(rnorm(1, mean=105, sd=5), digits=2)

    iqCutoff <- 100 + -qnorm(1/10000) * 15

    whites <- rnorm(0.770 * 579 * 10000, mean=whiteMean, sd=15)
    whiteSample <- max(1, sum(ifelse(whites>iqCutoff, 1, 0)))

    asians <- rnorm(0.007 * 579 * 10000, mean=asianMean, sd=15)
    asianSample <- max(1, sum(ifelse(asians>iqCutoff, 1, 0)))

    ## white+Asian = 92% of original total sample, so inflate by that much to preserve proportions: 1.08
    totalSample <- (whiteSample+asianSample) * (1 + (1-(white+asian)))

    whiteFraction <- round(whiteSample / totalSample, digits=2)
    asianFraction <- round(asianSample / totalSample, digits=2)
    # print(paste("samples: ", c(whiteSample, asianSample), "fractions: ", c(whiteFraction, asianFraction)))

    tolerance <- 0.02
    if ((abs(whiteFraction - 0.7218480138) < tolerance) && (abs(asianFraction - 0.2178929188) < tolerance)) {
      return(data.frame(White=whiteMean, Asian=asianMean))
    }
    }
library(parallel); library(plyr)
simulateSamples <- function(n.sample=10000, iters=getOption("mc.cores")) {
    ## because of rejection sampling, no run is guaranteed to produce a sample so we loop:
    results <- data.frame()
    while (nrow(results) < n.sample) {
        simResults <- ldply(mclapply(1:iters, function(i) { simulateTIPSMPY()  } ))
        results <- rbind(results, simResults)
        # print(paste("Samples so far: ", nrow(results)))
    }
    return(results) }
posteriorSamples <- simulateSamples()

mean(posteriorSamples$White < posteriorSamples$Asian)
# [1] 1
## we have relatively few samples, so get a better posterior estimate by shuffling the posterior samples & comparing many times:
mean(replicate(1000, mean(c(sample(posteriorSamples$White) < sample(posteriorSamples$Asian)))))
# [1] 0.9968822
quantile(probs=c(0.025, 0.975), posteriorSamples$White, na.rm=TRUE)
#     2.5%     97.5%
# 89.49975 101.38050
quantile(probs=c(0.025, 0.975), posteriorSamples$Asian, na.rm=TRUE)
#      2.5%     97.5%
# 101.37000 116.74075
par(mfrow=c(2,1))
hist(posteriorSamples$White, main="Posterior white mean IQ estimated from SMPY/TIP cutoff & ratio", xlab="IQ")
hist(posteriorSamples$Asian, main="Posterior Asian mean", xlab="IQ")