Initialization of stackingPredictiveFamily

[mfasiolo]

Currently, all the coefficients are initialized at 0, which means that every expert gets the same weight.

One important point is that family$initialize is used in initial.spg, so getting a good initialization is important for getting decent initial parameter estimates.

One possible initialization function that works only when discrete = FALSE can be found below. The idea is to do something similar to what is done in the multinomial family. The code would need to be adapted to work when discrete = TRUE, for example by looking at what Simon did in the gaulss family.

 initialize <- expression({

       if (is.null(start)){

       logP <- family$getLogP()
       ybest <- as.numeric(apply(logP, 1, which.max)) - 1

       ## should E be used unscaled or not?..
       use.unscaled <- if (!is.null(attr(E,"use.unscaled"))) TRUE else FALSE
      if (is.null(start)) {

        jj <- attr(x,"lpi")
        start <- rep(0,ncol(x))
        for (k in 1:length(jj)) { ## loop over the linear predictors
          yt1 <- 6 * as.numeric(ybest == k) - 3 
          x1 <- x[,jj[[k]],drop=FALSE]
          e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
          if (use.unscaled) {
            qrx <- qr(rbind(x1,e1))
            x1 <- rbind(x1,e1)
            startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
            startji[!is.finite(startji)] <- 0
          } else startji <- penReg(x1,e1,yt1)
          start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
        } ## lp loop
      }
    }

  }) ## initialize

[capezza]

an alternative, maybe, could be

1. to estimate a first model with the intercepts only and putting all the other coefficients to zero: this should improve all zero solution, or even better
2. estimate the intercepts and all the parametric terms included, putting the coefficients corresponding to smooth effects to zero, or maybe even better
3. estimate the model assuming everything is linear (ie replace smooths with parametric terms), so no smoothing parameters need to be estimated. then, as initial regression coefficients you set the ones you estimated and for smooths you set the coefficients that make the smooth effect linear with slope corresponding to the coefficient estimated.

[mfasiolo]

Maybe the easiest solution would be 1. This would require optimizing the objective:

obj <-  - sum(log( rowSums(exp(logP) * trans(beta0)) ))

over beta0 to get the intercepts (where trans() is the usual parametrization). But then it would be necessary to do regress

 jj <- attr(x,"lpi")
        start <- rep(0,ncol(x))
        for (k in 1:length(jj)) { ## loop over the linear predictors
          yt1 <- rep(b0[k], nobs)
          x1 <- x[,jj[[k]],drop=FALSE]
          e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
            qrx <- qr(rbind(x1,e1))
            x1 <- rbind(x1,e1)
            startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
            startji[!is.finite(startji)] <- 0
          start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
        } ## lp loop

To work out what the coefficients of each linear predictor should be (that is, to avoid having to find out which coefficient is the intercept and so on).

[mfasiolo]

For Christian:

small question:

yt1 <- 6 * as.numeric(ybest == k) - 3

why 6 and 3? if i understand well, you should get 6 if true and 3 if false. why not 1/0 instead?

in the following also i don't understand

yt1 <- rep(b0[k], nobs)

you are building a response with the estimated intercepts?

Anyway, i have thought to a possibility that could be closer to what you suggest. if i understand also on the basis of other families/methods, the objective should be to just use penalized least squares and qr decomposition to ensure to get a solution in one iteration, but this solution has to be reasonable in order to be a good initial value.
if, as you suggets, we do ybest <- as.numeric(apply(logP, 1, which.max)) - 1, we search WHICH is the max among the experts, and we treat the problem as a classification problem and we loose info about how much difference we have among experts, ie how much better the best expert is wrt the others.
since we have a convex combination, at each obs, the best we can get is the score of the best expert. then, why don't we use as response in the least squares the MAX value instead oh the WHICH.MAX? then the betas estimated should be the ones minimizing the squared distance between estimated log-score and the best log-score.
here however i am not precise. we have several linear predictors. to be more precise, notice that

alpha1 = 1 / sum
alpha2 = exp nu2 / sum
...
alphak = exp nuk / sum,

where sum = 1 + exp nu1 + ... + exp nuk and nuj is the j-th linear predictor

then nuj = log(alphaj/alpha1)

i.e. nuj is the log-ratio between the j-th expert and the first reference expert, then the j-th linear predictor just compares the j-th expert with the first.

then for the j-th linear predictor we can solve the penalized least squares problem
sum_i ||y_ij - nu_ij||^2 + penalty, where y_ij = max(p_i1, p_ij)

probably p_i1 and p_ij should be on the log-scale, i am not sure.

does it make sense to you?

[capezza]

sorry i have deleted the last comment, i am still elaborating on the final part because i made a mistake

[mfasiolo]

Regarding the 3 and 6 magic constants appearing in the multinomial initialization, the idea is to use:

yt1 <- 2 * const * as.numeric(ybest == k) - const

where we used 3 in this case. I don't know why, I tried to optimize the penalized likelihood over mult, but I am not sure it is worth it.

Regarding this:

yt1 <- rep(b0[k], nobs)

Yes, we are building a response with the estimated intercepts. It sounds stupid, but it is needed to set the whole linear predictor to b0[k] (when you don't where the intercept is among the coefficients of that linear predictor).

I think that maybe an idea would be trying to impose

nu_j = log(p_j) - log(p_1) 

Or possibly

nu_j = const * ( log(p_j) - log(p_1) ) 

Where const is a sort of learning rate, as in Opera, and we could optimized over it?

[capezza]

nu_j = log(p_j) - log(p_1)

yeah i was thinking to this now

[capezza]

actually i was thinking to put y = log(p_j) - log(p_1) and then to min |y-nuj|^2 + pen

[mfasiolo]

Something like the following?

 initialize <- expression({

       if (is.null(start)){

       logP <- family$getLogP()

       ## should E be used unscaled or not?..
       use.unscaled <- if (!is.null(attr(E,"use.unscaled"))) TRUE else FALSE
      if (is.null(start)) {

        jj <- attr(x,"lpi")
        start <- rep(0,ncol(x))
        for (k in 1:length(jj)) { ## loop over the linear predictors
          yt1 <- logP[ , k + 1] - logP[ , 1]
          x1 <- x[,jj[[k]],drop=FALSE]
          e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
          if (use.unscaled) {
            qrx <- qr(rbind(x1,e1))
            x1 <- rbind(x1,e1)
            startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
            startji[!is.finite(startji)] <- 0
          } else startji <- penReg(x1,e1,yt1)
          start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
        } ## lp loop
      }
    }

  }) ## initialize

[capezza]

exactly, or alternatively i was thinking to linearize max(p1,pj)-1/(1+expnuj)*p1 - expnuj/(1+expnuj)*pj wrt nuj around 0 and then solve sum of squares of this linearized difference

[capezza]

but this linearization idea doesn t look good. probably it s better the other option
yt1 <- logP[ , k + 1] - logP[ , 1]

[mfasiolo]

I agree, I am more clear about what the logP[ , k + 1] - logP[ , 1] solution does.

[capezza]

I tried this initialization on a simple example. The problem is that, by varying a coefficient k, if you do yt1 <- k * (logP[ , k + 1] - logP[ , 1]) you get different results.

Maybe it is better something like this:

yt1 <- ifelse(logP[ , k + 1] - logP[ , 1] > 0, 10, ifelse(logP[ , k + 1] - logP[ , 1] < 0, - 10, 0))
why? ideally we would like the linear predictor to be +infinite if logP[ , k + 1] > logP[ , 1] and -infinite if logP[ , k + 1] < logP[ , 1], such that, in the comparison of the two experts 1 and k only, we would have alpha1=0 alphak=1 if logP[ , k + 1] > logP[ , 1] and alpha1=1 alphak=0 if logP[ , k + 1] < logP[ , 1], but in a least squares problem we can't take yt1 arbitrarily large in absolute value, both because numerically we can have problems and because by varying the absolute value the result changes.
So we can take a fixed positive value if logP[ , k + 1] > logP[ , 1] and a negative value if logP[ , k + 1] < logP[ , 1].
I chose 10, why? because

a <- c(.99, .999, .9999, .9999, .99999)
log(a / (1 - a))
[1]  4.595120  6.906755  9.210240  9.210240 11.512915

in other words, 10 corresponds to the case alpha1 is in practice 0 and alphak is in practice 1.
I don't know if this makes sense, or if it's better to use yt1 <- k * (logP[ , k + 1] - logP[ , 1]) and find a way to set/estimate k

[mfasiolo]

One idea could be to optimize the penalized likelihood over k, as done in the example below
(which works only when we have 2 experts and 1 linear predictor, but could be extended). It adds a bit to the computational burden, but this should be marginal for non-trivial models.

I might be worth using something like this when unscaled == FALSE, in which case we are actually trying to get initial estimates of the regression coefficients.

But in initial.spg, the initialize() expression gets called with unscaled == TRUE, and there the idea is trying to get a good initial hessian. In particular, have a look at the code for gamFactory::penReg (same as mgcv::pen.reg) and you'll see that it is actually changing the penalty to try to balance the likelihood with the penalty (similarly to Wood's book second edition, section 6.7).

initialize <- expression({
       
    n <- rep(1, nobs)
    if (is.null(start)){
      
      logP <- family$getLogP()
      
      ## should E be used unscaled or not?..
      use.unscaled <- if (!is.null(attr(E,"use.unscaled"))) TRUE else FALSE
      
      if (is.null(start)) {
        
        # Objective function is penalized likelihood, mult is the multiplier logP[ , k + 1] - logP[ , 1]
        objFun <- function(mult){
          
          jj <- attr(x,"lpi")
          start <- rep(0,ncol(x))
          for (k in 1:length(jj)) { ## loop over the linear predictors
            yt1 <- mult * ( logP[ , k + 1] - logP[ , 1] )
            x1 <- x[,jj[[k]],drop=FALSE]
            e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
            if (use.unscaled) {
              qrx <- qr(rbind(x1,e1))
              x1 <- rbind(x1,e1)
              startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
              startji[!is.finite(startji)] <- 0
            } else startji <- penReg(x1,e1,yt1)
            start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
          } ## lp loop
          
          w <- cbind(1 / (1 + exp(x %*% start)), drop(exp(x %*% start)/(1 + exp(x %*% start))))
          out <-  - sum(log( rowSums(exp(logP) * w) )) + sum( (E%*%start) ^ 2 )
          return(out)
          
        }

        # Optimize, then initialize using the optimal mult
        tmpOpt <- optimize(f = objFun, interval = c(0.01, 100))
        mult <- tmpOpt$minimum
        
        jj <- attr(x,"lpi")
        start <- rep(0,ncol(x))
        for (k in 1:length(jj)) { ## loop over the linear predictors
          yt1 <- mult * ( logP[ , k + 1] - logP[ , 1] )
          x1 <- x[,jj[[k]],drop=FALSE]
          e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
          if (use.unscaled) {
            qrx <- qr(rbind(x1,e1))
            x1 <- rbind(x1,e1)
            startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
            startji[!is.finite(startji)] <- 0
          } else startji <- penReg(x1,e1,yt1)
          start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
        } ## lp loop
      }

    }
        
  }) ## initialize

[capezza]

I wrote the code and adapted it such that it works also in the discrete==TRUE case.
However, I get different results if I estimate a model using discrete==TRUE and discrete==FALSE. I have tried with an example with 1e6 observations, but since I just simulated some white noise I need to try with a more relevant example.

initialize <- expression({
    
    logP <- family$getLogP()
    P <- family$getP()
    
    ## should E be used unscaled or not?..
    use.unscaled <- if (!is.null(attr(E,"use.unscaled"))) TRUE else FALSE
    
    if (is.null(start)) {
      
      jj <- attr(x,"lpi")

      # Objective function is penalized likelihood, 
      # mult is the multiplier of ( logP[ , k + 1] - logP[ , 1] )
      objFun <- function(mult) {
        
        if (is.list(x)) { ## discrete case
          start <- rep(0, max(unlist(jj)))
          for (k in 1:length(jj)) { ## loop over the linear predictors
            yt1 <- mult * ( logP[ , k + 1] - logP[ , 1] )
            R <- suppressWarnings(
              chol(XWXd(x$Xd,w=rep(1,length(y)),k=x$kd,ks=x$ks,ts=x$ts,dt=x$dt,v=x$v,
                        qc=x$qc,nthreads=1,drop=x$drop,lt=x$lpid[[k]]) + 
                     crossprod(E[,jj[[k]]]),
                   pivot=TRUE))
            Xty <- XWyd(x$Xd,rep(1,length(y)),yt1,x$kd,x$ks,x$ts,x$dt,x$v,
                        x$qc,x$drop,lt=x$lpid[[k]])
            piv <- attr(R,"pivot")
            rrank <- attr(R,"rank")
            startji <- rep(0,ncol(R))
            if (rrank<ncol(R)) {
              R <- R[1:rrank,1:rrank]
              piv <- piv[1:rrank]
            }
            startji[piv] <- backsolve(R,forwardsolve(t(R),Xty[piv]))
            startji[!is.finite(startji)] <- 0
            start[jj[[k]]] <- startji
            
          } 
        } else {
          jj <- attr(x,"lpi")
          start <- rep(0,ncol(x))
          for (k in 1:length(jj)) { ## loop over the linear predictors
            yt1 <- mult * ( logP[ , k + 1] - logP[ , 1] )
            x1 <- x[,jj[[k]],drop=FALSE]
            e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
            if (use.unscaled) {
              qrx <- qr(rbind(x1,e1))
              x1 <- rbind(x1,e1)
              startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
              startji[!is.finite(startji)] <- 0
            } else startji <- penReg(x1,e1,yt1)
            start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
          } ## lp loop
        }
        
        # Produce matrix of linear predictors
        nu <- sapply(1:length(jj), 
                     function(ii) {
                       if (is.list(x)) {
                         Xbd(x$Xd,start,k=x$kd,ks=x$ks,ts=x$ts,dt=x$dt,v=x$v,qc=x$qc,
                             drop=x$drop,lt=x$lpid[[ii]])
                       } else {
                         x[, jj[[ii]], drop = FALSE] %*% start[jj[[ii]]]
                       }
                     })
        
        nu1 <- cbind(0, nu)
        nuCen <- nu1 - rowMaxs(nu1)
        a <- exp(nuCen) / rowSums(exp(nuCen))
        
        - sum(log(rowSums(a * P))) + sum( (E%*%start) ^ 2 )
        
      }
      
      # Optimize, then initialize using the optimal mult
      tmpOpt <- optimize(f = objFun, interval = c(0.01, 100))
      mult <- tmpOpt$minimum
      
      if (is.list(x)) { ## discrete case
        start <- rep(0,max(unlist(jj)))
        for (k in 1:length(jj)) { ## loop over the linear predictors
          yt1 <- mult * ( logP[ , k + 1] - logP[ , 1] )
          R <- suppressWarnings(
            chol(XWXd(x$Xd,w=rep(1,length(y)),k=x$kd,ks=x$ks,ts=x$ts,dt=x$dt,v=x$v,
                      qc=x$qc,nthreads=1,drop=x$drop,lt=x$lpid[[k]]) + 
                   crossprod(E[,jj[[k]]]),
                 pivot=TRUE))
          Xty <- XWyd(x$Xd,rep(1,length(y)),yt1,x$kd,x$ks,x$ts,x$dt,x$v,
                      x$qc,x$drop,lt=x$lpid[[k]])
          piv <- attr(R,"pivot")
          rrank <- attr(R,"rank")
          startji <- rep(0,ncol(R))
          if (rrank<ncol(R)) {
            R <- R[1:rrank,1:rrank]
            piv <- piv[1:rrank]
          }
          startji[piv] <- backsolve(R,forwardsolve(t(R),Xty[piv]))
          startji[!is.finite(startji)] <- 0
          start[jj[[k]]] <- startji
          
        } 
      } else {
        jj <- attr(x,"lpi")
        start <- rep(0,ncol(x))
        for (k in 1:length(jj)) { ## loop over the linear predictors
          yt1 <- mult * ( logP[ , k + 1] - logP[ , 1] )
          x1 <- x[,jj[[k]],drop=FALSE]
          e1 <- E[,jj[[k]],drop=FALSE] ## square root of total penalty
          if (use.unscaled) {
            qrx <- qr(rbind(x1,e1))
            x1 <- rbind(x1,e1)
            startji <- qr.coef(qr(x1),c(yt1,rep(0,nrow(E))))
            startji[!is.finite(startji)] <- 0
          } else startji <- penReg(x1,e1,yt1)
          start[jj[[k]]] <- startji ## copy coefficients back into overall start coef vector
        } ## lp loop
      }
    }
  })

[capezza]

After a first check, differences between discrete==TRUE and discrete==FALSE should not be relevant. Maybe this works! There is only the problem when you have only two experts and one linear predictor (I have opened an issue in your mgcv github repository)