CODEWITHSUNDEEP
roptimizationiterationglm
parastoo
parastoo

Reputation: 51

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$

I already managed optimization by setting parameter $c$ to arbitrary fixed values and only searching for $\hat{\beta}$.

BUT now for the next step I have to optimize this parameter $c$ as well (iteratively changing between updated $\beta$ and current $c$).

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? 

# 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)
}

# Running 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))

Upvotes: 0

Views: 62

Answers (1)

Paweł Chabros
Paweł Chabros

Reputation: 2399

You have 2 bugs:

line 25:

(y[i] * plogis(c * X[i,]%*%beta) / sp(X[i,]%*%beta, c) - plogis(c * X[i,]%*%beta))

this returns matrix so you must convert to numeric:

as.numeric(y[i] * plogis(c * X[i,]%*%beta) / sp(X[i,]%*%beta, c) - plogis(c * X[i,]%*%beta))

line 23: ) is missing:

you have:

s <- matrix(rep(NA, times = length(y)*length(par), ncol = length(y))

while it should be:

s <- matrix(rep(NA, times = length(y)*length(par)), ncol = length(y))

Upvotes: 1

Related Questions