R/decomp.R

In FastCUB: Fast Estimation of CUB Models via Louis' Identity

#' @title Louis' identity for the observed information matrix of the incomplete problem for CUB models
#' @description Compute the matrices involved in Louis' identity for the observed information matrix
#' @aliases decomp
#' @usage decomp(ttau,ordinal,m,param,ai,Y,W)
#' @param ttau Vector of posterior probabilities that each observation has been generated by the first mixture component (feeling)
#' @param ordinal Vector of ordinal responses
#' @param m Number of ordered categories
#' @param param  Vector of estimable parameters
#' @param ai  Auxiliary scalars
#' @param Y Matrix of selected covariates for explaining the uncertainty component
#' @param W Matrix of selected covariates for explaining the feeling component
#' @return A list of the matrices in Louis' identity for the observed information matrix
#' @import stats
#' @seealso \code{\link{fastCUB}}
#' @keywords internal 

########################################

########################################

decomp<-function(ttau,ordinal,m,param,ai,Y,W){
  
  np<-length(param)
  n<-length(ai) # in long form
  
  
  np<-length(param)
  n<-length(ordinal) # in long form
  
  
  if (is.null(Y)){
    p<-0
    paijj<-param[1]
    Y1<-rep(1/(paijj*(1-paijj)),n)
    YY<-as.matrix(Y1)
    
    bi<-paijj^2 + ttau*(1-2*paijj)
    
    if (is.null(W)){
      q<-0
      csijj<-param[2]
      W1<-rep(1/(csijj*(1-csijj)),n)
      WW<-as.matrix(W1)
      
    } else {
      q<-NCOL(W)
      W1<-rep(1,n)
      gamajj<-param[2:np]
      csijj<-logis(W,gamajj)
      WW<-cbind(W1,W)
      
      ui<-matrix(0,nrow=q+1,ncol=q+1)
      
      csicsi<-csijj*(1-csijj)
      
      Wcsi<-(m-1)*apply(WW,2,function(x) x*csicsi)
      
      
    }
    
    
    
  } else {
    
    p<-NCOL(Y)
    Y1<-rep(1,n)
    YY<-cbind(Y1,Y)
    betajj<-param[1:(p+1)]
    paijj<-logis(Y,betajj)   
    bi<-paijj*(1-paijj)
    
    #}
    
    if (is.null(W)){
      q<-0
      csijj<-param[p+2]
      W1<-rep(1/(csijj*(1-csijj)),n)
      WW<-as.matrix(W1)
      
      
    } else {
      q<-NCOL(W)
      W1<-rep(1,n)
      gamajj<-param[(p+2):np]
      csijj<-logis(W,gamajj)
      WW<-cbind(W1,W)
      
      
      csicsi<-csijj*(1-csijj)
      
      Wcsi<-(m-1)*apply(WW,2,function(x) x*csicsi)
    }
  } ### da togliere
  
  
  mat1<-mat2<-matrix(0,nrow=np,ncol=np)
  denom<-paijj*(1-paijj)
  
  
  vi<-WW
  
  for (s in 1:(q+1)){
    vi[,s]<- -ai*WW[,s]
  }
  
  
  # matrici dei prodotti scalari
  
  auxpai<-rep(0,p+1)
  for (t in 1:(p+1)){
    auxpai[t]<- t(YY[,t])%*%as.matrix((ttau-paijj))
    
  }
  
  auxcsi<-rep(0,q+1)
  for (t in 1:(q+1)){
    auxcsi[t]<- t(vi[,t])%*%(ttau)
    
  }
  
  for (t in 1:(p+1)){
    for (s in 1:(p+1)){
      mat2[t,s]<-auxpai[t]*auxpai[s]
    }
    for (j in 1:(q+1)){
      mat2[t,(p+1+j)]<-mat2[(p+1+j),t]<-auxpai[t]*auxcsi[j]
      for (l in 1:(q+1)){
        mat2[(p+1+l),(p+1+j)]<-auxcsi[j]*auxcsi[l]
        
      }
    }
  }
  mat2[(p+2):(p+q+2),1:(p+1)]<-t(mat2[1:(p+1),(p+2):(p+q+2)])
  
  #######################
  qq<-ttau*(1-ttau)
  Ytau<-apply(YY,2,function(x) x*qq)
  vitau<-apply(vi,2,function(x) x*qq)
  
  
  for (t in 1:(p+1)){
    for (s in 1:(p+1)){
      mat1[t,s]<-   t(Ytau[,t])%*%YY[,s]   +  mat2[t,s]
    }
    for (j in 1:(q+1)){
      mat1[t,(p+1+j)]<-mat1[(p+1+j),t]<- mat2[t,(p+1+j)] + t(Ytau[,t])%*%vi[,j] 
      for (l in 1:(q+1)){
        mat1[(p+1+l),(p+1+j)]<-  t(vitau[,j])%*%vi[,l]       + (t(vi[,l])%*%ttau)*(t(vi[,j])%*%ttau)
        
      }
    }
  }
  mat1[(p+2):(p+q+2),1:(p+1)]<-t(mat1[1:(p+1),(p+2):(p+q+2)])
  
  ########################################
  
  
  
  mat<-matrix(0,nrow=np,ncol=np)
  
  Yb<-apply(YY,2,function(x) x*bi)
  
  for (t in 1:(p+1)){
    for (s in 1:(p+1)){
      mat[t,s]<- t(Yb[,t])%*%YY[,s]
    }
  }
  for (t in 1:(q+1)){
    for (s in 1:(q+1)){
      
      if (is.null(W)){
        ui<- (ordinal-1)/(1-csijj)^2 + (m-ordinal)/(csijj^2)
      } else{
        ui<-as.matrix(Wcsi[,t]*WW[,s])
      }
      
      mat[p+1+t,p+1+s]<- t(ui)%*%ttau
    }
  }
  
  
  
  return(list('vcScore'=mat2,'vcScorec'=mat1,'Ic'=mat))
}