R/merlin_p_gompertz.R

In merlin: Mixed Effects Regression for Linear, Non-Linear and User-Defined Models

#predict fucntions for family(gompertz)

merlin_p_gompertz_h <- function(gml,t=NULL)
{
    if (length(t)==0) t = merlin_util_timevar(gml)
    gam = gml$par[gml$dapbindex[[gml$modelind]][1]]
    return(exp(merlin_util_xzb(gml,t) + gam * t))
}

merlin_p_gompertz_ch <- function(gml,t=NULL)
{
    if (length(t)==0) t = merlin_util_timevar(gml)

    if (gml$NI[[gml$modelind]]) {
        N   = gml$Nobs[gml$modtouse]
        gq 	= merlin_get_GK()
        qp	= sweep((t %*% (gq[[1]]/2)),1,(t/2),"+")
        if (gml$Nlevels & !gml$fixedonly) {
              result = matrix(0,nrow=N,ncol=gml$Ndim[gml$Nlevels])
        }
        else  result = matrix(0,nrow=N,ncol=1)
        for (q in 1:15) {
            result = result + sweep(merlin_p_gompertz_h(gml,qp[,q]),1,gq[[2]][q] * t/2,"*")
        }
        return(result)
    }
    else {
        gam = gml$par[gml$dapbindex[[gml$modelind]][1]]
        return(sweep(exp(merlin_util_xzb(gml,t)),1,(1/gam)*(exp(gam*t)-1),"*"))
    }
}

merlin_p_gompertz_logch <- function(gml,t=NULL)
{
    if (length(t)==0) t = merlin_util_timevar(gml)
    return(log(merlin_p_gompertz_ch(gml,t)))
}

merlin_p_gompertz_s <- function(gml,t=NULL)
{
    if (length(t)==0) t = merlin_util_timevar(gml)
    return(exp(-merlin_p_gompertz_ch(gml,t)))
}