1

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?

Chris Martin
  • 123
  • 2
  • 9
  • 1
    Without sample data this error isn't [reproducible](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). It it a lot easier to get help if you have a complete running example that returns the same error. See the link for suggestions for including a sample data set. – MrFlick Jun 01 '14 at 06:30
  • Sorry, I failed to copy one of the functions over correctly. Now put back in and it should run to reproduce the error mentioned. – Chris Martin Jun 01 '14 at 08:27
  • The sample data is generated by the code as dummy data that goes into the data frame 'theData'. – Chris Martin Jun 01 '14 at 08:28
  • R version 3.1.0, run in Ubuntu 12.04 and using RStudio. NO packages added. – Chris Martin Jun 01 '14 at 08:36

1 Answers1

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)     
}
MrFlick
  • 195,160
  • 17
  • 277
  • 295