R/ccrs.R

In ccrs: Correct and Cluster Response Style Biased Data

#' Correcting and Clustering response style biased data
#'
#' @export
#' @description Applies CCRS to \code{ccrsdata.list}.
#' @usage ccrs(ccrsdata.list,K=K,lam=lam, tandem.initial=FALSE,
#'             tol = 1e-5, maxit = 50, trace = 1, nstart = 3, parallel=F,verbose=T)
#' @param ccrsdata.list A list generated by \code{create.ccrsdata}.
#' @param K An integer indicating the number of content-based clusters used for CCRS estimation.
#' @param lam A numeric value indicating \code{lambda} used for CCRS estimation.
#' @param tandem.initial A logical value indicating whether the 1st initial value is generated by CCRS tandem initialization. See Section 3.3 in the paper for the detail.
#' @param tol A numeric value indicating the absolute convergence tolerance
#' @param maxit An integer indicating the maximum number of iterations
#' @param trace An non-negative integer. If positive, tracing information on the progress of the optimization is produced. Higher values produce more tracing information.
#' @param nstart An integer indicating the number of random initial values.
#' @param parallel A logical value indicating parallelization over starts is used.
#' @param verbose A logical value indicaitng if the progress is printed during the iteration (only when \code{parallel==FALSE}).
#' @return Returns a list with the following elements.
#' \item{\code{G}}{A K by m matrix of content-based cluster centroid.}
#' \item{\code{cls.cont.vec}}{A vector of integers (from 1:K) indicating the content-based cluster to which each respondent is allocated.}
#' \item{\code{opt.obval}}{An optimal value of objective function.}
#' \item{\code{crs.list}}{A list of class \code{crs}, same as the one generated by \link{correct.rs}.}
#' @seealso \code{\link{correct.rs}}
#' @references Takagishi, M., Velden, M. van de & Yadohisa, H. (2019). Clustering preference data in the presence of response style bias, to appear in British Journal of Mathematical and Statistical Psychology.
#' @importFrom lsbclust indarr
#' @importFrom parallel parLapply makeCluster detectCores clusterExport stopCluster
#' @importFrom stats kmeans runif
#' @export
#' @examples
#' ###data setting
#' n <- 30 ; m <- 10 ; H.true <- 2 ; K.true <- 2 ; q <- 5
#' datagene <- generate.rsdata(n=n,m=m,K.true=K.true,H.true=H.true,q=q,clustered.rs = TRUE)
#' ###obtain n x m data matrix
#' X <- datagene$X
#' ccrsdata.list <- create.ccrsdata(X,q=q)
#' ###CCRS
#' lam <- 0.8 ; K <- 2
#' ccrs.list <- ccrs(ccrsdata.list,K=K,lam=lam)
#' ###check content-based clustering result
#' ccrs.list$cls.cont.vec
#' ###check correction result
#' plot(ccrs.list$crs.list)


ccrs <- function(ccrsdata.list,K=K,lam=lam,tandem.initial=FALSE,
          tol = 1e-5, maxit = 50, trace = 1, nstart = 3,
        parallel=F,verbose=T){

  ####prepare data####
  X <- ccrsdata.list$X ; Fmat <- ccrsdata.list$Fmat
  Mmat.q1 <- ccrsdata.list$Mmat.q1  ; Mmat.q <- ccrsdata.list$Mmat.q
  n <- nrow(X) ; m <- ncol(X) ; q <- nrow(Mmat.q)
  Z <- indarr(X, maxcat = q)

  #maxit <- control$maxit ; tol <- control$tol ; nstart <- control$nstart ; trace <- control$trace

  #if(is.null(K)) stop("specify K!")

  ##other setting

  paraname=c("beta","cluster center","cluster allocation")
  updateorder=c(1,2,3)

  cat("K=",K,", lam=",lam,"\n")
  #cat("iteration: nstart=",nstart,", maxit=",maxit,"\n")



  ###prepare for iteration of init#####
  allstep <- (maxit+1)*length(paraname)
  OB.mat <- matrix(0,nstart,allstep)
  down.para.mat <- matrix(0,nstart,allstep)
  res.list <- rep(list(NA),nstart)
  optval.vec <- itevec <- rep(0,nstart)
  tot.list <- sapply(c(1:nstart),list)

  #######start iteration for initial values########
  #for(ti in 1:nstart){
  oneinit.func <- function(ti){

   # if(trace>1) cat("  ",ti,"th initial value\n")

    tandem.initial.ti <- ifelse(ti==1,tandem.initial,F)

    ######paralist
    Beta.para <- G.para <- U.para <- rep(list(NA),maxit)


    ####initial value for U,G (not for beta)
    if(tandem.initial.ti){#smoothing->kmeans

      crs.list <- correct.rs(ccrsdata.list)
      if(is.null(crs.list)){
        crs.list <- correct.rs(ccrsdata.list)
        #crs.list <- CCRS.lam1(Fmat=Fmat,Mmat.q1=Mmat.q1,Mmat.q=Mmat.q,X=X,
        #                        K=K,nstart.k=nstart.k)

      }

      #Beta <- crs.list$Beta
      #apply(crs.list$Beta,1,sum)
      Y.hat <- crs.list$Y.hat
      kres <- kmeans(Y.hat,K,nstart = 100)
      cls.init <- kres$cluster

      U.para[[1]] <- 1.0 * outer(cls.init, 1:K, "==")
      G.para[[1]] <- solve(t(U.para[[1]])%*%U.para[[1]])%*%t(U.para[[1]])%*%Y.hat

    }else if(!tandem.initial.ti){##end tandem.initial

      kinit.moto <- rep(c(1:K),n)[c(1:n)]
      cls.init <- kinit.moto[sample(n,n)]

      U.para[[1]] <- 1.0 * outer(cls.init, 1:K, "==")
      G.para[[1]] <- matrix(runif((m*K),0,(1)),K,m)
    }

    break.ite<-F# ; ite.ob<-1 ;
    OB.vec<-rep(0,allstep)
    down.para.vec<-rep(0,allstep)

    #############start iteration##############
    ite.ob<-1
    ite.b<-1 ; ite.u<-1 ; ite.g<-1
    ite<-1#trace<-T
    #iite<-2 ; ite<-iite ; ite.b<-iite ; ite.u<-iite; ite.g<-iite
    for(ite in 1:maxit){

      if(trace>1) cat("  ite=",ite,"\n")
      #browser()
      paran<-updateorder[1]
      for(paran in updateorder){

        if(trace>2) cat("    ite.b=",ite.b,",ite.u=",ite.u,",ite.g=",ite.g,"\n")

        if(paran==1){# beta update
          ####update quantification of categories###
          if(trace>2) cat(paste("update",paraname[paran]),"\n")
          ite.b<-ite+1

          ###update
          Beta.para[[ite.b]] <- update.beta(Fmat=Fmat,U=U.para[[ite.u]],G=G.para[[ite.g]],lam=lam,
                                              Z=Z,Mmat.q=Mmat.q,Mmat.q1=Mmat.q1)

          Y.hat <- transformRSdata(X,Beta=Beta.para[[ite.b]],Mmat.q=Mmat.q)$Y.hat

        }else if(paran==2){# cluster allpcation

          if(trace>2) cat(paste("update",paraname[paran]),"\n")
          #browser()
          #####update cluster center#####
          ite.g<-ite+1
          G.para[[ite.g]] <- update.G(U=U.para[[ite.u]],Y.hat=Y.hat)

        }else if(paran==3){# cluster allocation

          if(trace>2) cat(paste("update",paraname[paran]),"\n")

          #####update cluster center#####
          ite.u <- ite+1

          updateCluster.res <- update.U(G=G.para[[ite.g]],Y.hat=Y.hat)

          if(updateCluster.res$empty.cls){
            U.para[[ite.u]] <- U.para[[ite.u-1]]
            #G.para[[ite.g]]<-G.para[[ite.u]]
          }else{
            #G.para[[ite.g]]<-updateCluster.res$G
            U.para[[ite.u]] <- updateCluster.res$U
          }


        }

        para.list.now <- list(Beta=Beta.para[[ite.b]],U=U.para[[ite.u]],G=G.para[[ite.g]],
                            Fmat=Fmat,X=X,Mmat.q1=Mmat.q1, Mmat.q=Mmat.q, lam=lam)

        ###calculate ob
        OB.vec[ite.ob] <- OBJfunc(para.list.now)
        #OB.vec[ite.ob]<-ob.list$obval

        #obcheck.res<-objcheck.func(para.list=para.list.now,ite=ite,OB.vec=OB.vec#,trace=trace
        #                           ,paraname.p=paraname[paran],ite.ob=ite.ob,obcheck.start=1,tol=tol)

        #OB.vec<-obcheck.res$OB.vec
        #down.para.vec[ite.ob] <- obcheck.res$down.para.save
        #ite.ob<-obcheck.res$ite.ob

        #browser()
        # if(verbose) print(paste("obvalue",OB.vec[ite.ob-1]))
        if(verbose) cat("Start", formatC(ti, width = 5, digits = 0, mode = "integer"),
                        "Iter", formatC(ite, width = 5, digits = 0, mode = "integer"),
                        "Loss", formatC(OB.vec[ite.ob], width = 10, digits = 8, format = "f"),
                        "Diff", formatC(ifelse((ite==1 & paran==updateorder[1]),0,OB.vec[ite.ob]-OB.vec[ite.ob-1])
                                        , width = 10, digits = 12, format = "f"),
                        "\n")#only first cat, Dif is 0, otherwise print obval0-obval.

        break.ite <- F
        if(ite.ob>1){
          if(abs(OB.vec[ite.ob]-OB.vec[ite.ob-1]) < tol){
            OB.vec[(ite.ob+1):length(OB.vec)]<-OB.vec[ite.ob]
            break.ite <- T
            break
          }
        }
        ite.ob <- ite.ob+1


      }#####end paran######


      if(break.ite)  {
        if(trace>1) cat("  stop it!!!\n")
        break
      }


    }###################end iteration#########################

    #if(trace)
    #time.e <- proc.time()

    #if(trace>1) cat("  increased para:",paste(unique(down.para.vec),collapse="/"),"\n")
    if(trace>0) cat("  converged at",(ite-1),"th iteration.\n")


    #print(unique(down.para.vec))
    #plot(OB.vec,type="l",main=mainname.func(ti))
    opt.obval <- OB.vec[ite.ob-1]
    #if(ite.ob==1)opt.obval<-0

    Beta.p <- Beta.para[[ite.b]]
    G.p <- G.para[[ite.g]]
    #U.p <- U.para[[ite.u]]


    cls.cont.vec <- apply(U.para[[ite.u]],1,which.max)
    crs.list <- list(Beta=Beta.p,Y.hat=NULL,MB=NULL)

    return(list(G=G.p,cls.cont.vec=cls.cont.vec,opt.obval=opt.obval,crs.list=crs.list))
                #,OB.vec=OB.vec,down.para.vec=down.para.vec,step.conv=ite-1))#,ARIvec=ARIvec

  }#######################end initial value func#########################

  ## Apply algorithm over all starts, possibly in parallel
  if(parallel){

   # if (.Platform$OS.type == "windows") {
      cl <- makeCluster(detectCores())
      clusterExport(cl, varlist=c("tot.list","oneinit.func", #for paralell
                                  "tandem.initial","maxit","allstep","Fmat","lam","K",#common for any iteration
                                  "updateorder","paraname","crs.list","X","Xindi","Mmat.q","Mmat.q1",
                                  "n","m","OBJfunc","lsei","transformRSdata",
                                  "verbose",#"trace", #"plotcheck",
                                  #"objcheck.func"
                                  #"adjustedRandIndex", #clustering
                                  "update.beta","update.G","update.U" #method specific func.
      )
      ,envir=environment())#"calc.obj.func"

      #parallel::
      res.list <- parLapply(cl = cl, X = tot.list, fun = oneinit.func) #ed.m
      stopCluster(cl)

   # } else {
   #   res.list <- parallel::mclapply(X = tot.list, FUN = oneinit.func, mc.cores = detectCores(),
   #                                  mc.allow.recursive = FALSE)
   # }##end mac case
  } else {
    res.list <- lapply(X = tot.list, FUN = oneinit.func)
  }

  ti <- 1
  for(ti in 1:nstart){
    optval.vec[ti] <- res.list[[ti]]$opt.obval
  }


  ######decide minimum obvalue######
  res.list.min <- res.list[[which.min(optval.vec)]]

  #res.list.min$OB.mat <- OB.mat
  #res.list.min$optval.vec <- optval.vec
  #res.list.min$itevec <- itevec

  ######row and col names#####
  Beta.final <- res.list.min$crs.list$Beta
  rownames(Beta.final) <- paste("resp",c(1:n),sep="")
  colnames(Beta.final) <- c("intercept",paste("coef",c(1:3),sep=""))
  rownames(res.list.min$G) <- paste("cls",c(1:K),sep="")
  colnames(res.list.min$G) <- paste("item",c(1:m),sep="")

  ####create crs.list####
  trans.list <- transformRSdata(X,Beta=Beta.final,Mmat.q=Mmat.q)
  res.list.min$crs.list$Y.hat <- trans.list$Y.hat
  res.list.min$crs.list$MB <- trans.list$MB
  class(res.list.min$crs.list) <- "crs"

  #class(res.list.min) <- "ccrsres"
  res.list.min




}