Iterative optimization of alternative glm family I'm setting up an alternative response function to the commonly used exponential function in poisson glms, which is called softplus and defined as $\frac{1}{c} \log(1+\exp(c \eta))$, where $\eta$ corresponds to the linear predictor $X\beta$.
NOW I have to find a way to optimize the usual regression coefficients $\beta$
and the additional parameter $c$ of this new response function simultaneously and don't know how to solve it. (So I'm not asking for code review but theoretical suggestions how to manage this kind of iterative/simultaneous optimization)
I tried to write a log-lik function, score function and then setting up a Newton Raphson optimization (using a while loop) but I don't know how to seperate the updating of $c$ in an outer step and updating $\beta$ in an inner step.. 
Are there any suggestions? Is this procedure wrong? 
# Response function:
sp <- function(eta, c = 1 ) {  
  return(log(1 + exp(abs(c * eta)))/ c) 
} 

# Log Likelihood
l.lpois <- function(par, y, X){
  beta <- par[1:(length(par)-1)]
  c <- par[length(par)]
  l <- rep(NA, times = length(y))
  for (i in 1:length(l)){
    l[i] <- y[i] * log(sp(X[i,]%*%beta, c)) - sp(X[i,]%*%beta, c) 
  }
  l <- sum(l)
  return(l)
}

# Score function
score <- function(y, X, par){
  beta <- par[1:(length(par)-1)]
  c <- par[length(par)]

  s <- matrix(rep(NA, times = length(y)*length(par)), ncol = length(y))
  for (i in 1:length(y)){
    s[,i] <- c(X[i,], 1) * (y[i] * plogis(c * X[i,]%*%beta) / sp(X[i,]%*%beta, c) -     plogis(c * X[i,]%*%beta))
  }
  score <- rep(NA, times = nrow(s))
  for (j in 1:length(score)){
    score[j] <- sum(s[j,])
  }
  return(score)
}

# Optimization function
opt <- function(y, X, b.start, eps=0.0001, maxiter = 1e5){
  beta <- b.start[1:(length(b.start)-1)]
  c <- b.start[length(b.start)]

  b.old <- b.start
  i <- 0
  conv <- FALSE

  while(conv == FALSE){ 

    eta <- X%*%b.old[1:(length(b.old)-1)]
    s <- score(y, X, b.old)
    h <- numDeriv::hessian(l.lpois,b.old,y=y,X=X)

    invh <- solve(h)

    # update 
    b.new <- b.old + invh %*% s                                                         

    i <- i + 1

    # Test 
    if(any(is.nan(b.new))){                                                             
      b.new <- b.old                                                                
      warning("convergence failed")
      break 
    } 

    # convergence reached?
    if(sqrt(sum((b.new - b.old)^2))/sqrt(sum(b.old^2)) < eps | i >= maxiter){ 
      conv <- TRUE
    }
    b.old <- b.new
  }
  eta <- X%*%b.new[1:(length(b.new)-1)]

  # covariance
  invh  <- solve(numDeriv::hessian(l.lpois,b.new,y=y,X=X)) 


  fitted <- sp(eta, b.new[length(b.new)])

  result <- list("coefficients" = c(beta = b.new),
                 "fitted.values" = fitted,
                 "covariance" = invh)
}

# convergence fails ..
n <- 100
x <- runif(n, 0, 1)
Xdes <- cbind(1, x) 
eta <- 1 + 2 * x
y <- rpois(n, sp(eta, c = 1))



opt(y,Xdes,c(0,1,1)) 

 A: Parametric link functions
The "softplus" model you have proposed is called a parametric link function in generalized linear model theory. Technically, the link function is the inverse of your softplus function.
I can think of three major ways to estimate parametric link functions.
I'll list them in decreasing order of computational speed and programming difficulty.
Nested Fisher scoring
Mathematical algorithms to estimate parametric link functions have been given by Scallan et al (1984) and Smyth (1996). In Example 4.4 of Smyth (1996), I gave a very efficient algorithm for estimating the linear predictor parameters (your $\beta$) in an inner loop and the link parameter (your $c$) in an outer loop.
The idea of inner and outer loops is called nested iterations in my paper and my algorithm is called nested Fisher scoring.
Generalized nonlinear models
Parametric link functions are a special case of the larger class of Generalized nonlinear models (GNMs).
You will find a number of publications on that topic if you do a Google Scholar search.
Rather than programming a super-efficient algorithm, for which you really do need to be an R expert, you might find it simpler to use the gnm package on CRAN, which provides a code approach for GNMs.
Profile likelihood
The third way to estimate $c$ is take a profile likelihood approach.
For each fixed values of $c$ you have an ordinary GLM with a nonstandard link function.
For given $c$, fit the GLM and compute the log-likelihood or AIC.
As a function of $c$, this called a profile likelihood.
You can then use the optimize function in R to maximize the profile likelihood with respect to the single parameter $c$.
You can implement the GLM fit using R's built-in glm function.
Just copy the poisson(link="log") family and edit the link and link derivative entries to be appropriate for your softplus function.
Then you can specify your own softplus link when you run glm.
All of these approach are much better than trying to work with numerical derivatives.
References
Scallan, A., Gilchrist, R. and Green, M. (1984). Fitting parametric link functions
in generalised linear models. Comput. Statist. Data Anal., 2, 37–49.
Smyth, G. K. (1996). Partitioned algorithms for maximum likelihood and other nonlinear estimation. Statistics and Computing, 6, 201–216.
http://www.statsci.org/smyth/pubs/partitio.pdf
