Response to selection in teosinte assimilation

Data

Download data frrom Dolores's growouts. 360ppm and 215ppm branch lengths (3 years data lumped together).

system("wget -O b360.txt https://gist.github.com/rossibarra/5964105/raw/71adcff6492367aa7608ca561df3c4a6f4c2716f/bl360")
bl360 <- scan("b360.txt")
system("wget -O b215.txt https://gist.github.com/rossibarra/5964126/raw/980f82bc928e45529c12ccb0dc74075e74a6eea4/bl215")
bl215 <- scan("b215.txt")

Estimates \( \sigma^2_P \approx0.03 \) and \( \mu\approx 0.42 \) for branch length at 360ppm. Also look at data.

Vp = var(bl215)
Vp

## [1] 0.03434

bmu = mean(bl215)
bmu

## [1] 0.3707

hist(bl215, xlab = "branch length 215ppm", main = "")

Vp = var(bl360)
Vp

## [1] 0.0298

bmu = mean(bl360)
bmu

## [1] 0.4165

hist(bl360, xlab = "branch length 360ppm", main = "")

Estimate response to selection

\( R=h^2s \) is easy enough. No idea what \( h^2 \) or \( s \) are, so try lots. But \( s=\frac{s_f + s_m}{2} \) and because we phenotype after flowering, there is no selection on males (\( s_m=0 \)) and \( s=\frac{s_f}{2} \). To get \( s \) we decide how many shortest plants we choose. We'll try \( n=3,6,9,12... \) from our population of 24 plants and estimate \( s \) as if we were doing the equivalent sampling from the data (32 plants). Note that this \( n \) corresponds roughly to \( N_e \)!

h2 = 1:10/10  # 0-100% h^2 in 5%
prop = c(1:4)/8  # prop adults chosen
Ne = prop * 24  # Ne in our theoretical experimental pop
sorted_bl360 = bl360[order(bl360)]  # sort observed data
prop = prop * 32  # equivalent n in observed data

get_s <- function(x) {
    mean(bmu - sorted_bl360[1:x])
}  # gets mean of shortest branch length teos of size x
s = sapply(prop, function(x) get_s(x))/2  # divide by two because of s=sf/2 above

resp <- function(h2, s) {
    h2 * s
}
R = sapply(h2, function(x) resp(h2 = x, s = s))
rownames(R) = Ne
colnames(R) = h2

library(lattice)
wireframe(R, scales = list(arrows = F), xlab = expression(N[e]), ylab = expression(h^2))

Calculate sampling variance in Response to selection

Following p. 210 in Falconer & Mackay to estimate sampling variance in response to selection.

\( \sigma^2_R=\sigma^2_P*\frac{t*h^2}{N_e}+\frac{1}{M} \)

Assume \( t=1 \) generations. Total plants measured \( M=24 \).

We will have two replicate pops of \( M=24 \) but assume here single pop of 24 plants.

Vr <- function(Ne, h2, Vp = 0.03, t, M = 24) {
    Vp * (t * h2/Ne + 1/M)
}

Have to use range of \( h^2 \) as above, and now substitute in \( N_e \) from above too.

est_vr = sapply(Ne, function(x) Vr(h2 = h2, Ne = x, t = 1))
rownames(est_vr) = h2
colnames(est_vr) = Ne
wireframe(est_vr, scales = list(arrows = F), screen = list(z = 230, x = -70, 
    y = 0), zlab = expression(sigma[R]^2), xlab = expression(h^2), ylab = expression(N[e]))

Confidence intervals

With \( h^2 \) and \( s \) (and our bad estimate of \( N_e \) we can get \( R \) and \( \sigma^2_R \) which let's us calcualte confidence interval.

Response to selection. Same \( h^2 \) and \( N_e \) as before.

resp <- function(h2, s) {
    h2 * s
}
R = sapply(h2, function(x) resp(h2 = x, s = s))
rownames(R) = s
colnames(R) = h2

Get \( \sigma^2_R \)

# Uses Vr function from above
est_vr = sapply(h2, function(x) Vr(h2 = x, Ne = Ne, t = 1))
rownames(est_vr) = Ne
colnames(est_vr) = h2

Get lower bound of Confidence Interval (Rmin) assuming Confidence Interval is \( R \pm 2* \sigma_R \). We want \( Rmin>0 \) in order to guarantee positive response to selection.

sd_R = est_vr^0.5
Rmin = R - 2 * sd_R
rownames(Rmin) = Ne
colnames(Rmin) = h2
wireframe(Rmin, zlab = list("Rmin", rot = 90), scales = list(arrows = F), xlab = expression(N[e]), 
    ylab = expression(h^2))

Higher \( N_e \)

Now try with assumption that \( N_e \) is fourfold higher,but all else stays same:

fNe = Ne * 4
est_vr = sapply(h2, function(x) Vr(h2 = x, Ne = fNe, t = 1))  #var. response
rownames(est_vr) = fNe
colnames(est_vr) = h2
sd_R = est_vr^0.5
Rmin = R - 2 * sd_R
rownames(Rmin) = fNe
colnames(Rmin) = h2
wireframe(Rmin, scales = list(arrows = F), zlab = list("Rmin_4N_e", rot = 90), 
    xlab = expression(N[e]), ylab = expression(h^2), screen = list(z = 230, 
        x = -70, y = 0))

Works for high \( h^2 \).

Try over multiple generations

Assume \( t=5 \). Here it gets a bit odd, as we don't know how strong selection (difference between final and initial mean) we can really do. Let's assume 3x original difference, but \( N_e \) stays same.

R = sapply(h2, function(x) resp(h2 = x, s = 3 * s))  #re-estimate response with 3X s 
est_vr = sapply(h2, function(x) Vr(h2 = x, Ne = Ne, t = 5))  #var. response now over 5 gens
rownames(est_vr) = Ne
colnames(est_vr) = h2
sd_R = est_vr^0.5
Rmin = R - 2 * sd_R
rownames(Rmin) = Ne
colnames(Rmin) = h2
wireframe(Rmin, zlab = list("Rmin_5gens", rot = 90), scales = list(arrows = F), 
    xlab = expression(N[e]), ylab = expression(h^2))

Conclusion

Not so great outlook. \( Rmin<0 \) for almost all \( h^2 \) unless we dramatically increase \( N_e \) or can select hard for same \( N_e \). Results look even worse (essentially never \( Rmin<0 \) always) for data from 215ppm.