R and nls with functions error - 'must have same number of rows'

Question

I am wanting to find optimum parameters for a function that determines survival from hazard rates (force of mortality).

eh <- function(kappa,lambda,time,delay){
  FoM <- -log(1-(1-exp(-kappa*(time+delay)))*(exp(-lambda*time)))  
  return(FoM)  
}

survival <- function(k, l, t, d){
  theSurvs <- array(1, c(1, length(t)))
  S = array(1)
  Qs <- array(0)
  for(j in 1:length(t)){
    S[1] = 1
    for(i in 1:(t[j]+1)){
      theEHs <- eh(k, l, i-1, d)
      theQs <- hazards2Qs(theEHs)
      S[i+1] <- S[i]-S[i]*theQs
    }
    theSurvs[j] <- S[i]
  }
  return(theSurvs)
}

hazards2Qs <- function(hazards){
  Qs<- 1-exp(-hazards)
  return(Qs)
}

timeX = 0  # time cycles of arbitrary length
kappaX = 0.001
lambdaX = 0.05
delayX = 50

theYears <- floor(runif(20)*100)
theYears

EHs <- eh(kappaX, lambdaX, theYears, delayX)


theS <- survival(kappaX, lambdaX, theYears, delayX)
theS

theData <- data.frame(theYears, t(theS))
theData
survival(kappaX, lambdaX, theYears, delayX)

mod <- nls(~ survival(kappa, lambda, theYears, delay), start=list(lambda=0.0015, kappa=0.0013, delay=48), data=theData, trace=1)

When I run it I get this error:

3.647752 :   0.0015  0.0013 48.0000
Error in qr.qty(QR, resid) : 
  'qr' and 'y' must have the same number of rows

traceback():

4: stop("'qr' and 'y' must have the same number of rows")
3: qr.qty(QR, resid)
2: (function () 
   {
       if (npar == 0) 
           return(0)
       rr <- qr.qty(QR, resid)
       sqrt(sum(rr[1L:npar]^2)/sum(rr[-(1L:npar)]^2))
   })()

> dim(theS)
[1]  1 20
> dim(theYears)
NULL
> dim(theData)
[1] 20  2
> typeof(theData)
[1] "list"
> typeof(theS)
[1] "double"
> typeof(theYears)
[1] "double"

I have been slogging for a day and not got to the bottom of this. Any ideas?

Answer 1

After spending a long time trying to figure out what was going on myself, I think the problem is how you're shaping your survival function results. I'm not sure why you're trying to return an array, but you should be returning a vector here. So just change the first line to

#OLD: theSurvs <- array(1, c(1, length(t)))
theSurvs <- rep.int(1, length(t)

which means you'll also have to change

#OLD: theData <- data.frame(theYears, t(theS))
theData <- data.frame(theYears, theS)

By returning a shaped object like that, it was interfering with the gradient calculation which uses qr(). Try to stay away from single dimension arrays when possible and just use simple vectors.

Now, that being said, fixing that problem seems to lead to another one. It seems that when it's trying to take the numerical derivative, it's running into NaN/Inf values. You can add this code above your return(theSurvs) line to see the parameters that are called when this happens

if(any(!is.finite(theSurvs))) {
    dput(c(k,l,d))
    dput(t)
    print(theSurvs)     
}