Fit logistic growth model like "nlsList" but allow user to specify a design matrix for each parameter

logistic_funcs.R

### New growth function
logistic_growth = function(x, phi1, phi2, phi3){
  phi1/(1 + exp(-(x-phi2)/phi3))
} 

ll_logistic = function(par, y, xvar, phi1_dm, phi2_dm, phi3_dm, sigma_dm, sigma=1, predict=FALSE, scale=-2){
  n1 = ncol(phi1_dm)
  n2 = ncol(phi2_dm)
  n3 = ncol(phi3_dm)
  if(!missing(sigma_dm)) n4 = ncol(sigma_dm)
  beta1 = par[1:n1]
  beta2 = par[(n1+1):(n1+n2)]
  beta3 = par[(n1+n2+1):(n1+n2+n3)]
  if(!missing(sigma_dm)) betas = par[(n1+n2+n3+1):(n1+n2+n3+n4)]
  phi1 = exp(phi1_dm%*%beta1)
  phi2 = phi2_dm%*%beta2
  phi3 = exp(phi3_dm%*%beta3)
  if(!missing(sigma_dm)) sigma=exp(sigma_dm%*%betas)
  lg_mean = logistic_growth(xvar, phi1, phi2, phi3)
  if(predict){
    return(lg_mean)
  } else{
    return(scale*sum(dnorm(y, lg_mean, sigma, log = TRUE)))
  }
}

check_dm = function(X){
  if(ncol(X)==1){
    return(X)
  } else{
    zeros = apply(X, 2, mean)==0
    const = apply(X,2,var)==0
    X = X[,!(zeros&const)]
    dup_col = duplicated(t(X))
    return(X[,!dup_col])
  }
}

proj = function(X,y){
  solve(t(X)%*%X, t(X)%*%y)
}

fit_logistic = function(formula, model_list = list(phi1=~1, phi2=~1, phi3=~1, sigma_group=~1), data, quiet=FALSE){
  require(numDeriv)
  xvar = model.frame(formula, data=data)[,2]
  y = model.frame(formula, data=data)[,1]
  phi1_dm = check_dm(model.matrix(model_list$phi1, data))
  phi2_dm = check_dm(model.matrix(model_list$phi2, data))
  phi3_dm = check_dm(model.matrix(model_list$phi3, data))
  sigma_dm = check_dm(model.matrix(model_list$sigma_group, data))
  # start
  phi1 = max(y) + 1
  y_tilde = log((y/phi1)/(1 - y/phi1))
  cc = coef(lm(y_tilde~phi2_dm + phi3_dm:xvar - 1))
  cc_2 = cc[grep("phi2_dm",names(cc))]
  cc_3 = cc[grep("phi3_dm",names(cc))]
  phi2 = -(phi2_dm%*%cc_2)/(phi3_dm%*%cc_3)
  beta2 = proj(phi2_dm, phi2)
  phi3 = 1/(phi3_dm%*%cc_3)
  beta3 = proj(phi3_dm, log(phi3))
  beta1 = proj(phi1_dm, log(rep(phi1, length(y))))
  start=c(beta1, beta2, beta3)
  # Fitting
  if(!quiet) message("Getting start values for the parameters ...")
  start = optim(start, ll_logistic, y=y, xvar=xvar, 
                phi1_dm=phi1_dm, phi2_dm=phi2_dm, phi3_dm=phi3_dm,
                method="SANN", control = list(maxit=10000))
  start = optim(start$par, ll_logistic, y=y, xvar=xvar, 
                phi1_dm=phi1_dm, phi2_dm=phi2_dm, phi3_dm=phi3_dm,
              method="BFGS", control = list(maxit=10000))
  start = optim(start$par, ll_logistic, y=y, xvar=xvar, 
                phi1_dm=phi1_dm, phi2_dm=phi2_dm, phi3_dm=phi3_dm, 
                method="SANN", control = list(maxit=10000))
  start = optim(start$par, ll_logistic, y=y, xvar=xvar, 
                phi1_dm=phi1_dm, phi2_dm=phi2_dm, phi3_dm=phi3_dm, 
              method="BFGS", control = list(maxit=10000))
  sigma_dm = check_dm(model.matrix(model_list$sigma_group, data=data))
  fitted = ll_logistic(start$par, y=y, xvar, phi1_dm, phi2_dm, phi3_dm, predict=TRUE)
  betas = proj(sigma_dm, log(abs(y-fitted)))
  start=c(start$par, betas)
  if(!quiet) message("Optimizing likelihood ...")
  mle = optim(start, ll_logistic, y=y, xvar=xvar, 
                phi1_dm=phi1_dm, phi2_dm=phi2_dm, phi3_dm=phi3_dm, 
                sigma_dm=sigma_dm,
                method="BFGS", control = list(maxit=10000))
  H = hessian(ll_logistic, x=mle$par, y=y, xvar=xvar, 
              phi1_dm=phi1_dm, phi2_dm=phi2_dm, phi3_dm=phi3_dm, 
              sigma_dm=sigma_dm)
  vcov=2*solve(H)
  n4 = ncol(sigma_dm)
  betas = tail(mle$par, n4)
  sigma2 = as.vector(exp(2*sigma_dm%*%betas))
  Jac = jacobian(ll_logistic, mle$par, y=y, xvar=xvar, 
             phi1_dm=phi1_dm, phi2=phi2_dm, phi3=phi3_dm, predict=TRUE)
  n1 = ncol(phi1_dm)
  n2 = ncol(phi2_dm)
  n3 = ncol(phi3_dm)
  n4 = ncol(sigma_dm)
  beta_idx = list( 1:n1, (n1+1):(n1+n2), (n1+n2+1):(n1+n2+n3), (n1+n2+n3+1):(n1+n2+n3+n4))
  results=vector(mode="list",length=2)
  names(results)=c("beta", "phi")
  results$beta = vector("list",4)
  results$phi = vector("list",3)
  names(results$beta) = c("phi_1","phi_2","phi_3", "sigma")
  names(results$phi) = c("phi_1","phi_2","phi_3")
  for(i in 1:4){
    idx=beta_idx[[i]]
    beta=as.vector(mle$par[idx])
    # Beta table
    se=as.vector(sqrt(diag(vcov)[idx]))
    z=as.vector(beta/se)
    p_val = 1-pnorm(z)
    results$beta[[i]] = data.frame(beta=beta, SE=se, Z=z, P_val=p_val)
  }
  # phi 1
  idx = beta_idx[[1]]
  results$phi$phi_1 = model.frame(model_list$phi1, data)
  results$phi$phi_1$phi_1 = exp(phi1_dm%*%mle$par[idx])
  results$phi$phi_1$SE = exp(phi1_dm%*%mle$par[idx])*sqrt(diag(phi1_dm%*%vcov[idx,idx]%*%t(phi1_dm)))
  results$phi$phi_1 = unique(results$phi$phi_1)
  # phi 2
  idx = beta_idx[[2]]
  results$phi$phi_2 = model.frame(model_list$phi2, data)
  results$phi$phi_2$phi_2 = phi2_dm%*%mle$par[idx]
  results$phi$phi_2$SE = sqrt(diag(phi2_dm%*%vcov[idx,idx]%*%t(phi2_dm)))
  results$phi$phi_2 = unique(results$phi$phi_2)
  # phi 3
  idx = beta_idx[[3]]
  results$phi$phi_3 = model.frame(model_list$phi3, data)
  results$phi$phi_3$phi_3 = exp(phi3_dm%*%mle$par[idx])
  results$phi$phi_3$SE = exp(phi3_dm%*%mle$par[idx])*sqrt(diag(phi3_dm%*%vcov[idx,idx]%*%t(phi3_dm)))
  results$phi$phi_3 = unique(results$phi$phi_3)
  fitted=data.frame(fitted=fitted, SE=sqrt(diag(Jac%*%vcov%*%t(Jac))))
  parameter_design = list(phi1_dm=phi1_dm, phi2_dm=phi2_dm,phi3_dm=phi3_dm,sigma_dm=sigma_dm)
  out=list(model_list=model_list, data=data, fitted=fitted, results=results,
           par=mle$par, vcov=vcov, beta_idx=beta_idx, parameter_design=parameter_design,
           logLik = ll_logistic(mle$par, y=y, xvar=xvar, 
                                phi1_dm=phi1_dm, phi2=phi2_dm, phi3=phi3_dm,
                                sigma_dm = sigma_dm, scale=1)
  )
           
  class(out) = c("logistic_model", class(out))
  return(out)
}

LRT = function(fit_full, fit_red){
  x2 = -2*(fit_red$logLik-fit_full$logLik)
  df = length(fit_full$par)-length(fit_red$par)
  pval=pchisq(x2, df=df, lower.tail = FALSE)
  cat("LRT: chisq = ", x2, "; df = ", df, "; P = ", round(pval,4), "\n")
}