R functions used in demonstration on 12/8
These functions can be cut-and-pasted into your R environment.
The function driftandm simulates a stepping-stone model of 100 populations
evolving according to a phenotypic optimum selection model with genetic drift
and migration. N is local population size, m is the migration rate (m/2 from
each adjacent stepping-stone), f is the fraction of the way to the optimum
that each population moves under selection.
The optimum has a linear cline with phenotype 10 in the middle of the
species and a slope of b per population.
driftandm <- function(N,m,f,b) {
# genetic drift, selection, and migration in a linear cline of optima
y <- 10.0+b*(seq(1:100)-50.5); # optimum phenotypes
x <- y; # start out with means sitting on optima
for(t in 1:1000) {
x <- x + rnorm(100, mean=0, sd=sqrt(1/N)); # genetic drift of means
xl <- x[1:99]; # migration code starts here ...
xr <- x[2:100];
x[1] <- (1-m/2)*x[1];
x[100] <- x[100]*(1-m/2);
x[2:99] <- (1-m)*x[2:99];
x[1:99] <- x[1:99]+(m/2)*xr;
x[2:100] <- x[2:100]+(m/2)*xl; # ... and ends here
x <- x+f*(y-x); # optimum selection: go f of way to optimum
}
plot(seq(1:100),x,ylim=c(0,20),ylab="character mean",xlab="population");
lines(seq(1:100),x);
lines(seq(1:100),y);
}
The function regressm runs 100 replicates under the above model. It then
plots the confidence intervals from a simple regression of population means
on population position, and draws a line for the 95% interval. If the
interval fails to include a slope of zero, it is colored red and drawn with
a double-thickness line.
To see whether the intervals include b = 0 the correct fraction of the time,
run this function with the argument b, the true slope, being zero. When
m is zero, the interval should have the correct coverage, but when m is
greater than zero, it should not owing to residuals being correlated owing to
migration.
regressm <- function(N,m,f,b) {
y <- 10+b*(seq(1:100)-50.5);
plot(c(0,100),c(0,0),type="l",xlim=c(1,100),ylim=c(-0.1,0.1),xlab="replicate",
ylab="confidence interval on slope");
for (reps in 1:100) {
x <- y; # start right on top of optima
for(t in 1:1000) {
x <- x + rnorm(100, mean=0, sd=sqrt(1/N)); # genetic drift
xl <- x[1:99]; # migration ...
xr <- x[2:100];
x[1] <- (1-m/2)*x[1];
x[100] <- x[100]*(1-m/2);
x[2:99] <- (1-m)*x[2:99];
x[1:99] <- x[1:99]+(m/2)*xr;
x[2:100] <- x[2:100]+(m/2)*xl; # ... done migration
x <- x+f*(y-x); # selection: move fraction f of way to optimum
}
fit <- lm(x~seq(1:100));
coef <- summary(fit)$coefficients[2,1] ;
err <- summary(fit)$coefficients[2,2];
b <- coef + c(-1,1)*err*qt(0.975, 98);
if ((b[1] > 0.0) | (b[2] < 0.0))
lines(c(reps,reps),c(b[1],b[2]),col="red",lwd=2,ylim=c(-1,1),ylab="character mean",xlab="population")
else
lines(c(reps,reps),c(b[1],b[2]),ylim=c(-1,1),ylab="character mean",xlab="population")
}
}