fitdiff: Compute invasion fitness difference of learning strategy

Description Usage Arguments Details Author(s) Examples

Description

Used with sim_birds to calculate selection gradients on social learning strategies.

Usage

1
fitdiff(simdat,simdat2,mutant=1)

Arguments

simdat

Result of sim_birds that includes a mutant individual

simdat2

Result of sim_birds that includes only common-type strategies

mutant

Position of mutant in simdat

Details

This function simulates groups of birds learning under the assumed statistical model. It produces data suitable to feed directly into the statistical model, for purpose of validating the analysis code. It can also be used to explore a variety of learning dynamics.

Author(s)

Richard McElreath

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# calculate gradient over s and conf
s_seq <- c( 0 , 0.05 , 0.1 , 0.15 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 )
l_seq <- c( 1 , 1.25 , 1.5 , 1.75 , 2.0 , 2.5 , 3.0 )
seq <- expand.grid( s=s_seq , l=l_seq , dl=NA , dlse=NA , ds=NA , dsse=NA )

N_bird <- 10
tmax <- 100
dhi_set <- rep(c(2.5,1),each=tmax/2)
dlo_set <- rep(c(1,2.5),each=tmax/2)

# function to compute common-type fitness
fsx0 <- function(ks,kl) sim_birds( N_bird , tmax , 
        s=rep( ks[1] , N_bird ) , 
        gamma=rep( 0.6 ,N_bird ) , 
        conf=rep( kl[1] , N_bird ) , 
        pay=rep(0,N_bird) , 
        d1=dhi_set , d2=dlo_set )
# function to compute mutant fitness in lambda* (kl)
fsx1 <- function(ks,kl) sim_birds( N_bird , tmax , 
        s=c( ks[1] , rep( ks[1] ,N_bird-1) ) , 
        gamma=c( 0.6 , rep( 0.6 ,N_bird-1) ) , 
        conf=c( kl[2] , rep( kl[1] ,N_bird-1) ) , 
        pay=c( 0 , rep(0,N_bird-1) ) , 
        d1=dhi_set , d2=dlo_set )
# function to compute mutant fitness in s* (ks)
fsx2 <- function(ks,kl) sim_birds( N_bird , tmax , 
        s=c( ks[2] , rep( ks[1] ,N_bird-1) ) , 
        gamma=c( 0.6 , rep( 0.6 ,N_bird-1) ) , 
        conf=c( kl[1] , rep( kl[1] ,N_bird-1) ) , 
        pay=c( 0 , rep(0,N_bird-1) ) , 
        d1=dhi_set , d2=dlo_set )
# compute fitness difference between individuals
fitdiff <- function(simdat,simdat2,mutant=1) {
    tmax <- nrow(simdat)
    n <- ncol(simdat)
    #non_mutant <- (1:n)[-mutant]
    non_mutant <- 1:n # common-types in simdat2
    d1 <- attr(simdat,"d1")
    d2 <- attr(simdat,"d2")
    p <- 0
    p2 <- 0
    for ( i in 1:tmax ) {
        p <- p + (simdat[i,]*d1[i])
        p <- p + ((1-simdat[i,])*d2[i])
        p2 <- p2 + (simdat2[i,]*d1[i])
        p2 <- p2 + ((1-simdat2[i,])*d2[i])
    }
    Wmutant <- mean(p[mutant])
    Wnmutant <- mean(p2[non_mutant])
    return( Wmutant - Wnmutant )
}

# my multi-core replicate wrapper
mcrepsilent <- function (n, expr, mc.cores = 2) 
{
    result <- simplify2array(mclapply(1:n, eval.parent(substitute(function(i, 
        ...) {
        expr
    })), mc.cores = mc.cores))
    result
}

# this is the top-level function to compute fitness differential
frep <- function(ks,kl,i) {
    n_rep <- 1e4
    # lambda gradient
    xl <- mcrepsilent( n_rep , fitdiff( fsx1(ks,kl) , fsx0(ks,kl) ) , mc.cores=63 )
    # s gradient
    xs <- mcrepsilent( n_rep , fitdiff( fsx2(ks,kl) , fsx0(ks,kl) ) , mc.cores=63 )
    print( concat("[",i,"] ",ks[1]," ",kl[1]) )
    return( c( mean(xl) , sd(xl)/sqrt(n_rep) , mean(xs) , sd(xs)/sqrt(n_rep) ) )
}

# DO NOT RUN THIS, UNLESS YOU HAVE CLUSTER OR A LOT OF TIME
# we used a 64 core cluster and it took 8 or 9 hours
# instead use data(wythamgradient) to load the resulting object into namespace
if ( FALSE ) {
time_start <- Sys.time()
for ( i in 1:nrow(seq) ) {
    seq[i,3:6] <- frep( c(seq[i,1],seq[i,1]+0.1) , c(seq[i,2],seq[i,2]+0.5) , i )
}
time_end <- Sys.time()
( time_end - time_start ) # report elapsed time

}

# write.csv( seq , file="gradientx.csv" , row.names=FALSE )
# can also use data(wythamgradient) to load this result into namespace
data(wythamgradient)
seq <- wythamgradient

#####################
# FIGURE 6
# make contour plot
blank(w=3)
par(mfrow=c(1,3))

n <- length(unique(seq$s))
m <- matrix( seq$dl , nrow = n )
contour( unique(seq$s) , unique(seq$l) , m , xlab="social learning weight" , ylab="conformity strength" )
mtext( "dl" )
lev0 <- contourLines(unique(seq$s) , unique(seq$l) , m , levels=0 )
isol <- lev0
for( i in 1:length(lev0) ) lines( lev0[[i]] , col="red" , lwd=2 )

m <- matrix( seq$ds , nrow = n )
contour( unique(seq$s) , unique(seq$l) , m , xlab="social learning weight" , ylab="conformity strength" )
mtext("ds")
lev0 <- contourLines(unique(seq$s) , unique(seq$l) , m , levels=0 )
isos <- lev0
for( i in 1:length(lev0) ) lines( lev0[[i]] , col="blue" , lwd=2 )

# make vector field plot
plot(NULL, type = "n", xlim=range(seq$s) , ylim=range(seq$l) , xlab="social learning" , ylab="conformity strength" )
fz <- function(x) sign(x)*abs(x)^(0.5)/20
arrows(seq[,1], seq[,2], seq[,1] + fz(seq[,5]), seq[,2] + fz(seq[,3]) , length=0.05 )

# add isoclines
for( i in 1:length(isol) ) lines( isol[[i]] , col="red" , lty=2 )
for( i in 1:length(isos) ) lines( isos[[i]] , col="blue" , lty=2 )

rmcelreath/wythamewa documentation built on May 14, 2019, 2:56 p.m.