############################################################################# # # # Exemplary code for the estimation of plant neighbourhood models # # as applied in the Journal of Ecology paper # # # # "Neighbourhood effects on plant reproduction: an experimental-analytical # # framework and its application to the invasive Senecio inaequidens" # # (Authors: Susanne Lachmuth, Colette Henrichmann, Juliane Horn, # # Jörn Pagel, Frank M. Schurr) # # # ############################################################################# # Installing the NHM package from source files in NHM_1.1.0.tar.gz install.packages('NHM_1.1.0.tar.gz', repos = NULL, type="source") # (depends on packages: methods, MASS, Matrix, coda, mnormt) # Loading the NHM package library(NHM) ############################################################################# # Data preparation # Loading data for target and source individuals target.data <- read.table("target.data.txt",header=TRUE,sep="\t") source.data <- read.table("source.data.txt",header=TRUE,sep="\t") # Selection of target individuals # Neighbourhood models target.ind <- which(target.data$treatment == 'open') target.data <- target.data[target.ind,] # Log-transformation and scaling of explanatory variables # in target.data target.data$sc.log.capitula <- log(target.data$capitula+1)/sd(log(target.data$capitula+1)) target.data$sc.log.shoots <- log(target.data$shoots)/sd(log(target.data$shoots)) target.data$sc.log.height <- log(target.data$height)/sd(log(target.data$height)) # in source.data source.data$sc.log.capitula <- log(source.data$capitula+1)/sd(log(source.data$capitula+1)) source.data$sc.log.shoots <- log(source.data$shoots+1)/sd(log(source.data$shoots+1)) source.data$sc.log.height <- log(source.data$height+1)/sd(log(source.data$height+1)) source.data$sc.log.diameter <- log(source.data$diameter+1)/sd(log(source.data$diameter+1)) source.data$sc.log.cover <- log(source.data$cover+1)/sd(log(source.data$cover+1)) # Transforming the variable 'round' into factor type target.data$round <- factor(target.data$round) # Adding a column of zeros to the target data # (see specification of 'response' model below) target.data$zeros <- 0 # The runNHM() functions can handle data from multiple sites and expects the data # for target and source individuals each as a list with one element holding the # data of one site. # In this example, however, all the data come from one site (named 'Zwenkau'): target.list <- list(Zwenkau = target.data) source.list <- list(Zwenkau = source.data) ############################################################################# # Specification and parameter estimation of the neighbourhood model # List of model specifications nhm.model <- list() # The function runNHM() requires a specification of the hierachical model structure # that it is provided as an argument of the format # list( observation = list(data, family, link), local = list(model, link), # response = list(model, link), effect = list(model, link), # distance = list(model, link), max.range ) # In this list, the element 'observation' defines the measured response variable ('data') # of the model, its error distribution ('family') and the link function ('link'): nhm.model$observation <- list(data = c('fertile.seeds', 'infertile.seeds'), family = 'binomial', link = 'logit') # The element 'local' defines a linear model formula for the effect of target plant traits # on the measured response variable (term 'b' in the main text): nhm.model$local <- list(model = ~ sc.log.height + sc.log.shoots + sc.log.capitula + round) # The element 'response' defines a linear model formula for the effect of target plant traits # on the strength with which a target plant responds to neighbourhood effects (term 'r'). # Here, response stength is fixed at a constant exp(0) => 1: nhm.model$response <- list(model = ~ zeros - 1, link ='log' ) # The element 'effect' defines a linear model formula for the effect of source plant traits # on the strength of their neighbourhood effect (term 'e' in the main text): nhm.model$effect <- list(model = ~ sc.log.height + sc.log.shoots + sc.log.capitula + sc.log.diameter + sc.log.cover + spec.fact - 1) # The element 'distance' defines a linear model formula for the effect of source plant traits # on the distance decay of their neighbourhood effect (term 'f' in the main text). # Here, distance decay is estimated constant across all source plants nhm.model$distance <- list(model = ~ 1, link ='log') # The element 'max.range' defines a maximum distance between target and source plant # for which neighbourhood effects are considered. # This is defined mainly for computational efficiency: nhm.model$max.range <- 12 # Specifications of MCMC parameters mcmc.para <- list(ni = 400000, nc = 2, nb = 0, nt = 10) # ni: total number of MCMC iteration # nc: number of independent MCMC chains # nb: number of iterations discarded as burn-in # nt: thinning of the stored MCMC chains # Running MCMC posterior sampling with runNHM() nhm.out <- runNHM(Target = target.list, Source = source.list, Model = nhm.model, mcmc.para = mcmc.para) ############################################################################# # The output of runNHM() is an object of class 'NHMsample'. It has a slot @Samples # that contains the MCMC samples as a mcmc.list object, to which methods of # the coda package can be applied: # Discarding of burn-in iterations using window.mcmc() nhm.sample.burn <- window(nhm.out@Samples$Zwenkau, start = 300000) # Ploting MCMC chains with plot.mcmc() plot(nhm.sample.burn) # (Note that occassionally we observed numerically degenerated MCMC chains where the sampling # algorithm apparently got stuck, in which case the runNHM() function had to be restarted.) # Posterior summaries using summary.mcmc() summary(nhm.sample.burn) #############################################################################