R和nls的函数错误——“必须有相同数量的行”

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?

我辛辛苦苦地干了一天，还没弄清这件事的底细。什么好主意吗?

1 个解决方案

#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.

通过返回一个形状的对象，它干扰了使用qr()的梯度计算。在可能的情况下尽量远离单一维度的数组，只使用简单的向量。

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

现在，有人说，解决这个问题似乎会导致另一个问题。似乎当它试图求数值导数时，它会跑进NaN/Inf值。您可以将此代码添加到return(theSurvs)行之上，以查看在发生这种情况时调用的参数。

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

#1

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