R/cub0q.R

Defines functions cub0q

Documented in cub0q

#' @title Main function for CUB models with covariates for the feeling component
#' @description Function to estimate and validate a CUB model for given ordinal responses, with covariates for
#'  explaining the feeling component.
#' @aliases cub0q
#' @usage cub0q(m, ordinal, W, maxiter, toler)
#' @param m Number of ordinal categories
#' @param ordinal Vector of ordinal responses
#' @param W Matrix of selected covariates for explaining the feeling component, not including intercept
#' @param maxiter Maximum number of iterations allowed for running the optimization algorithm 
#' @param toler Fixed error tolerance for final estimates 
#' @import stats graphics
#' @return An object of the class "CUB"
#' @references
#'  Piccolo D. and D'Elia A. (2008), A new approach for modelling consumers' preferences,
#'   \emph{Food Quality and Preference}, \bold{18}, 247--259 \cr
#' @references 
#' Iannario M. and Piccolo D. (2010), A new statistical model for the analysis of customer
#'  satisfaction, #' \emph{Quality Technology and Quantity management}, \bold{7}(2) 149--168 \cr
#' Iannario M. and Piccolo D. (2012), CUB models: Statistical methods and empirical evidence, in:
#' Kenett R. S. and Salini S. (eds.), \emph{Modern Analysis of Customer Surveys: with applications using R},
#'  J. Wiley and Sons, Chichester, 231--258. 
#' @keywords internal 


cub0q<-function(m,ordinal,W,maxiter,toler){
  tt0<-proc.time()
  n<-length(ordinal)
  W<-as.matrix(W)
  if (ncol(W)==1){
    W<-as.numeric(W)
  }
  q<-NCOL(W)
  aver<-mean(ordinal)
  WW<-cbind(1,W)                   
  ##############################################################
  freq<-tabulate(ordinal,nbins=m); inipaicsi<-inibest(m,freq);  paijj<-inipaicsi[1]; 
  gamajj<-inibestgama(m,ordinal,W)
  ##############################################################
  loglikjj<-loglikcub0q(m,ordinal,W,paijj,gamajj)
  # ********************************************************************
  # ************* E-M algorithm for CUB(0,q) ***************************
  # ********************************************************************
  nniter<-1
  while(nniter<=maxiter){

    loglikold<-loglikjj
    vettn<-bitgama(m,ordinal,W,gamajj)
    ttau<-1/(1+(1-paijj)/(m*paijj*vettn)) 
    ################################# maximize w.r.t. gama  ########
    ordd<-ordinal;covar<-WW;
    gama<-gamajj
    optimgama<-optim(gama,effe01,esterno01=cbind(ttau,ordinal,WW),m=m) 
    ################################################################
    gamajj<-optimgama$par
    paijj<-sum(ttau)/n                    #updated pai estimate
 
    
    loglikjj<-loglikcub0q(m,ordinal,W,paijj,gamajj)## needed for nlm version
    # print(c(nniter,paijj,gamajj,loglikjj)); #OPTIONAL PRINTING OF ITERATIONS
    testll<-abs(loglikjj-loglikold)
    
   
    if(testll<=toler) break else {loglikold<-loglikjj}
    nniter<-nniter+1
  }
  pai<-paijj;  gama<-gamajj;  loglik<-loglikjj;
  ####################################################################
  AICCUB0q<- -2*loglik+2*(q+2)
  BICCUB0q<- -2*loglik+log(n)*(q+2)
  ####################################################################
  # Compute asymptotic standard errors of ML estimates
  ####################################################################
  varmat<-varcovcub0q(m,ordinal,W,pai,gama)
  nomi<-c("pai    ",paste("gamma",0:(length(gama)-1),sep="_"))
  stime<-c(pai,gama)
  nparam<-length(stime)
  if (isTRUE(varmat==matrix(NA,nrow=nparam,ncol=nparam))==TRUE){
    ddd<-cormat<-matrix(NA,nrow=nparam,ncol=nparam)
    ICOMP<-trvarmat<-NA
    errstd<-wald<-pval<-rep(NA,nparam)
  } else {
    ddd<-diag(sqrt(1/diag(varmat)))
    cormat<-(ddd%*%varmat)%*%ddd  
    trvarmat<-sum(diag(varmat))
    ICOMP<- -2*loglik + nparam*log(trvarmat/nparam) - log(det(varmat)) ## added
    errstd<-sqrt(diag(varmat));  wald<-stime/errstd;
    pval<-2*(1-pnorm(abs(wald)))
  }
  rownames(cormat)<-nomi;colnames(cormat)<-nomi; 
  durata<-proc.time()-tt0;durata<-durata[1];
  
  results<-list('estimates'=stime,'ordinal'=ordinal,'time'=durata,
                'loglik'=loglik,'niter'=nniter,'varmat'=varmat,
                'BIC'=BICCUB0q)
  # class(results)<-"cub"
  return(results)
}

Try the CUB package in your browser

Any scripts or data that you put into this service are public.

CUB documentation built on March 31, 2020, 5:14 p.m.