7.3 S-plus version

fitmodel does, as already mentioned, use nlminb to carry out the actual numerical maximization. In addition to the name of the function that is to be optimized, nlminb also needs vectors specifying starting parameter values and upper and lower box constraints. The fitmodel function assumes that these vectors will be returned through additional calls to FUN with one of the optional optional argument start, lower, or upper set true. The following call illustrates the correct behaviour of steppingstone:

   > steppingstone(upper=T)
   [1] 1 1

The source code to achieve this in the case of the steppingstone model is as follows:

   steppingstone <- function(theta,n,upper=F,lower=F,start=F)
   {
      if (upper) {return(c(1.0,   1.0 )) }
      if (lower) {return(c(0.0001,0.0 )) }
      if (start) {return(c(0.1,   0.1 )) }

      M <- matrix(0,ncol=n,nrow=n)
      diag(M) <- 1 - theta[2]
      M[abs(row(M)-col(M))==1] <- theta[2]/2
      M[1,1] <- M[n,n]  <- 1 - theta[2]/2
      M <- (1-theta[1])*M
      return(M)
   }

Only when none of the optional arguments are true (which is default) is the main code executed and a migration matrix returned. Although this perhaps isn't very good programming practice, the main advantage of implementing things this way is that good starting and box constraint vectors are defined once and for all, and that functions carrying out, for example, bootstrapping and likelihood ratio tests can retrieve this information only knowing the name of the migration matrix submodel. This will prove very useful when working with several alternative models.