R/plot.LMlatent.R

In LMest: Generalized Latent Markov Models

plot.LMlatent<-function(x, what = c("modSel", "CondProb", "transitions","marginal"),verbose=interactive(),...)
  {
  object <- x
  what <- match.arg(what,several.ok = TRUE)
  n <- object$n; TT <- object$TT; k <- object$k
  d <- dim(object$Psi)
  nc <- d[1]
  if(length(d)==2) r <- d[2]
  else  r <- d[3]

  if(k>1){
    PMarg <- array(0,c(n,k,TT))
    PMarg[,,1] <- as.matrix(object$Piv)
    for(i in 1:n) for(t in 2:TT) PMarg[i,,t]= t(object$PI[,,i,t])%*%PMarg[i,,t-1]
    Pmarg <-apply(PMarg,c(2,3),mean)
  }
  if(interactive() & length(what) > 1){
    title <- "LM model plots:"
    # present menu waiting user choice
    choice <- menu(what, graphics = FALSE, title = title)
    while(choice != 0){
      if(what[choice] == "modSel"){
        par(mar=c(5, 4, 4, 2) + 0.1,mfrow=c(1,1))
        if(is.null(object$Aic)){
          object$Aic <- object$aic
          names(object$Aic) <- paste("k",k,sep = "=")
          object$Bic <- object$bic
          names(object$Bic) <- paste("k",k,sep = "=")
        }
        ylim1 =  min(c(object$Aic,object$Bic))
        ylim2 = max(c(object$Aic,object$Bic))*1.05
        kv <-as.numeric(substr(names(object$Bic),3,3))
        matplot(kv,cbind(object$Bic,object$Aic),type="b",pch=1:2,lty=1,col=c(2,4),xaxt="n",xlab="Number of states",ylab="",ylim = c(ylim1,ylim2))
        #lines(object$Aic,type="b",pch=2,col=4)
        axis(side=1,at=kv)
        legend("topright",legend=c("BIC","AIC"),col=c(2,4),lty=1,bty="n")
      }
      if(what[choice] == "CondProb"){
        title = "Conditional response probabilities"
        if(r==1){
          if(k==1){
            par(mfrow=c(1,1),mar=c(4,4,2,0)+0.1)
            pi.class <- matrix(NA,nrow=r,ncol=nc)
            pi.class[1,] <- object$Psi
            ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
            tt =as.table(ds.plot$value)
            dimnames(tt)[[1]]=0:(nc-1)
            plot(tt, type = "h", col = "blue", lwd = 10,col.axis="blue",bty="n",main=title,ylab="Prob.Cond.",xlab="Item categories",ylim=c(0,1))
          }else{
            par(mfrow=c(k,1),mar=c(4,4,2,0)+0.1)
            for (u in 1:k) {
              title = paste("Class ",u,": Averaged initial probabilites = ",round(mean(object$Piv[,u]),3),sep="")
              pi.class <- matrix(NA,nrow=r,ncol=nc)
              pi.class[1,] <- object$Psi[,u,]
              ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
              tt =as.table(ds.plot$value)
              dimnames(tt)[[1]]=0:(nc-1)
              plot(tt, type = "h", col = "blue", lwd = 10,col.axis="blue",bty="n",main=title,ylab="Prob.Cond.",xlab="Item categories",ylim=c(0,1))
            }
          }
        }else{
          if(k==1){
            par(mfrow=c(1,1),mar=c(5,4,4,2)+0.1)
            pi.class <- matrix(NA,nrow=r,ncol=nc)
            for (j in 1:r) pi.class[j,] <- object$Psi[,j]
            dimnames(pi.class) <- list(item=1:r,category=0:(nc-1))
            ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
            vis <- scatterplot3d(ds.plot,type="h",lwd=7,pch=" ",x.ticklabs=1:r,y.ticklabs=colnames(pi.class),z.ticklabs = "",
                                 xlab="Items",zlab="Cond. prob.",main=title,lab=c(r-1,nc-1),zlim=c(0,1)
                                 ,mar=c(3,3,2,3), highlight.3d=FALSE,col.grid="lightblue",col.axis="lightblue",cex.main=1.5,angle=75,box=FALSE,color=4)
          }else{
            par(mfrow=c(k,1),mar=c(5,4,4,2)+0.1)
            for (u in 1:k) {
              title = paste("Class ",u,": Averaged initial probabilites = ",round(mean(object$Piv[,u]),3),sep="")
              pi.class <- matrix(NA,nrow=r,ncol=nc)
              for (j in 1:r) pi.class[j,] <- object$Psi[,u,j]
              dimnames(pi.class) <- list(item=1:r,category=0:(nc-1))
              ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
              vis <- scatterplot3d(ds.plot,type="h",lwd=7,pch=" ",x.ticklabs=1:r,y.ticklabs=colnames(pi.class),z.ticklabs = "",
                                   xlab="Items",zlab="Cond. prob.",main=title,lab=c(r-1,nc-1),zlim=c(0,1)
                                   ,mar=c(3,3,2,3), highlight.3d=FALSE,col.grid="lightblue",col.axis="lightblue",cex.main=1.5,angle=75,box=FALSE,color=4)
            }
          }
        }
      }
      if(what[choice] == "transitions"){
        if(k==1){
          stop("Transition probabilities not available when k=1")
        }

        par(mar=c(2,1,5,1),mfrow=c(1,1))
        PM<-round(apply(object$PI[,,,2:TT],c(1,2),mean),3)
        PM = round(diag(1/rowSums(PM))%*%PM,3)
        plotmat(t(PM),relsize=0.7,box.col="lightblue",lwd = 1,
                box.lwd = 1,self.cex = 0.8,
                cex.txt = 0.8, box.size = 0.1,box.prop = 0.5,main="Averaged transition probabilities")
      }
      if(what[choice] == "marginal") {
        if(k==1){
          stop("Marginal distribution of the latent states not available when k=1")
        }
        par(mar=c(5, 4, 4, 2) + 0.1,mfrow=c(1,1))
        plot(1:TT,Pmarg[1,],ylim=c(0,1),xaxt="n",xlab="Time",ylab="Estimated marginal distribution",type="l",lwd=1)
        for(u in 2:k) lines(Pmarg[u,],col=u,lwd=1)
        leg = paste("state",1:k)
        axis(side=1,at=1:TT)
        legend(x="topright",leg,col=1:k,lwd=1,bty="n")
      }

      # re-present menu waiting user choice
      choice <- menu(what, graphics = FALSE, title = title)
    }
  }
  else{
    if(any(what == "modSel")){
    par(mar=c(5, 4, 4, 2) + 0.1,mfrow=c(1,1))
      if(is.null(object$Aic)){
        object$Aic <- object$aic
        names(object$Aic) <- paste("k",k,sep = "=")
        object$Bic <- object$bic
        names(object$Bic) <- paste("k",k,sep = "=")
      }
    ylim1 =  min(c(object$Aic,object$Bic))
    ylim2 = max(c(object$Aic,object$Bic))*1.05
    kv <-as.numeric(substr(names(object$Bic),3,3))
    matplot(kv,cbind(object$Bic,object$Aic),type="b",pch=1:2,lty=1,col=c(2,4),xaxt="n",xlab="Number of states",ylab="",ylim = c(ylim1,ylim2))
    #lines(object$Aic,type="b",pch=2,col=4)
    axis(side=1,at=kv)
    legend("topright",legend=c("BIC","AIC"),col=c(2,4),lty=1,bty="n")
  }
    if(any(what == "CondProb")){
      title = "Conditional response probabilities"
      if(r==1){
        if(k==1){
          par(mfrow=c(1,1),mar=c(4,4,2,0)+0.1)
          pi.class <- matrix(NA,nrow=r,ncol=nc)
          pi.class[1,] <- object$Psi
          ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
          tt =as.table(ds.plot$value)
          dimnames(tt)[[1]]=0:(nc-1)
          plot(tt, type = "h", col = "blue", lwd = 10,col.axis="blue",bty="n",main=title,ylab="Prob.Cond.",xlab="Item categories",ylim=c(0,1))
        }else{
          par(mfrow=c(k,1),mar=c(4,4,2,0)+0.1)
          for (u in 1:k) {
            title = paste("Class ",u,": Averaged initial probabilites = ",round(mean(object$Piv[,u]),3),sep="")
            pi.class <- matrix(NA,nrow=r,ncol=nc)
            pi.class[1,] <- object$Psi[,u,]
            ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
            tt =as.table(ds.plot$value)
            dimnames(tt)[[1]]=0:(nc-1)
            plot(tt, type = "h", col = "blue", lwd = 10,col.axis="blue",bty="n",main=title,ylab="Prob.Cond.",xlab="Item categories",ylim=c(0,1))
          }
        }
      }else{
        if(k==1){
          par(mfrow=c(1,1),mar=c(5,4,4,2)+0.1)
          pi.class <- matrix(NA,nrow=r,ncol=nc)
          for (j in 1:r) pi.class[j,] <- object$Psi[,j]
          dimnames(pi.class) <- list(item=1:r,category=0:(nc-1))
          ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
          vis <- scatterplot3d(ds.plot,type="h",lwd=7,pch=" ",x.ticklabs=1:r,y.ticklabs=colnames(pi.class),z.ticklabs = "",
                               xlab="Items",zlab="Cond. prob.",main=title,lab=c(r-1,nc-1),zlim=c(0,1)
                               ,mar=c(3,3,2,3), highlight.3d=FALSE,col.grid="lightblue",col.axis="lightblue",cex.main=1.5,angle=75,box=FALSE,color=4)
        }else{
          par(mfrow=c(k,1),mar=c(5,4,4,2)+0.1)
          for (u in 1:k) {
            title = paste("Class ",u,": Averaged initial probabilites = ",round(mean(object$Piv[,u]),3),sep="")
            pi.class <- matrix(NA,nrow=r,ncol=nc)
            for (j in 1:r) pi.class[j,] <- object$Psi[,u,j]
            dimnames(pi.class) <- list(item=1:r,category=0:(nc-1))
            ds.plot <- data.frame(Items=as.vector(row(pi.class)),Categories=as.vector(col(pi.class)),value=as.vector(pi.class))
            vis <- scatterplot3d(ds.plot,type="h",lwd=7,pch=" ",x.ticklabs=1:r,y.ticklabs=colnames(pi.class),z.ticklabs = "",
                                 xlab="Items",zlab="Cond. prob.",main=title,lab=c(r-1,nc-1),zlim=c(0,1)
                                 ,mar=c(3,3,2,3), highlight.3d=FALSE,col.grid="lightblue",col.axis="lightblue",cex.main=1.5,angle=75,box=FALSE,color=4)
          }
        }
      }

    }
    if(any(what == "transitions")) {
      if(k==1){
        stop("Transition probabilities not available when k=1")
      }
      par(mar=c(2,1,5,1),mfrow=c(1,1))
      TT <- dim(object$PI)[4]
      PM<-round(apply(object$PI[,,,2:TT],c(1,2),mean),3)
      PM = round(diag(1/rowSums(PM))%*%PM,3)
      plotmat(t(PM),relsize=0.7,box.col="lightblue",lwd = 1,
              box.lwd = 1,self.cex = 0.8,
              cex.txt = 0.8, box.size = 0.1,box.prop = 0.5,main="Averaged transition probabilities")
    }
  }
  if(any(what == "marginal")) {
    if(k==1){
      stop("Marginal distribution of the latent states not available when k=1")
    }
    par(mar=c(5, 4, 4, 2) + 0.1,mfrow=c(1,1))
    plot(1:TT,Pmarg[1,],ylim=c(0,1),xaxt="n",xlab="Time",ylab="Estimated marginal distribution",type="l",lwd=1)
    for(u in 2:k) lines(Pmarg[u,],col=u,lwd=1)
    axis(side=1,at=1:TT)
    leg = paste("state",1:k)
    legend(x="topright",leg,col=1:k,lwd=1,bty="n")
  }
}