R/se.R

Defines functions se.LMlatentcont se.LMbasiccont se.LMlatent se.LMbasic se

Documented in se se.LMbasic se.LMbasiccont se.LMlatent se.LMlatentcont

se <- function(est,...) {
  UseMethod("se")
}

se.LMbasic <- function(est, ...){

  check_der = FALSE
  if(!is.null(est$sepiv))
  {
    out = list(sepiv = est$sepiv,sePi = est$sePi,sePsi = est$sePsi)
    return(out)
  }

  data <- est$data
  responsesFormula = attributes(est)$responsesFormula
  latentFormula = attributes(est)$latentFormula

  tv.which <- attributes(est)$whichtv
  id.which <- attributes(est)$whichid
  id <- attributes(est)$id
  tv <- attributes(est)$time
  k <- est$k
  mod <- est$modBasic

  data.new <- data[,-c(id.which,tv.which), drop = FALSE]

  if(is.null(responsesFormula))
  {
    Y <- data.new
    Xmanifest <- NULL
    Xinitial <- NULL
    Xtrans <- NULL
  }else{
    temp <- getResponses(data = data.new,formula = responsesFormula)
    Y <- temp$Y
    Xmanifest <- temp$X
    Xinitial <- NULL
    Xtrans <- NULL
  }
  if(min(Y,na.rm=T)>0){
    Y <- apply(Y,2,function(x) x - min(x, na.rm = TRUE))
  }

  tmp <- long2matrices.internal(Y = Y, id = id, time = tv, yv = NULL,
                                Xinitial = Xinitial, Xmanifest = Xmanifest, Xtrans = Xtrans)
  #model <- tmp$model
  S <- tmp$Y
  yv = tmp$freq

  miss = any(is.na(S))
  if(miss){
    R = 1 * (!is.na(S))
    Y[is.na(S)] = 0
  }else{
    R = NULL
  }

  piv = est$piv
  Pi = est$Pi
  Psi = est$Psi

  n = sum(yv)
  sS = dim(S)
  ns = sS[1]
  TT = sS[2]

  if(length(sS)==2){
    r = 1
    if(is.matrix(S)) S = array(S,c(dim(S),1))
  }else
  {
    r = sS[3]
  }

  Sv = matrix(S,ns*TT,r)

  bv = apply(Sv,2,max)
  b = max(bv)
  m = vector("list",r)

  for(j in 1:r){
    m[[j]]$A = cbind(-rep(1,bv[j]),diag(bv[j]))
    m[[j]]$Am = rbind(rep(0,bv[j]),diag(bv[j]))
  }

  B = cbind(-rep(1,k-1),diag(k-1))
  Bm = rbind(rep(0,k-1),diag(k-1))
  C = array(0,c(k-1,k,k))
  Cm = array(0,c(k,k-1,k))
  for(u in 1:k){
    C[,,u] = rbind(cbind(diag(u-1),-rep(1,u-1),matrix(0,u-1,k-u)),
                   cbind(matrix(0,k-u,u-1),-rep(1,k-u),diag(k-u)))
    Cm[,,u] = rbind(cbind(diag(u-1),matrix(0,u-1,k-u)),
                    rep(0,k-1),
                    cbind(matrix(0,k-u,u-1),diag(k-u)))
  }



  # Compute information matrix if required
  th = NULL
  for(u in 1:k) for(j in 1:r) th = c(th,m[[j]]$A%*%log(Psi[1:(bv[j]+1),u,j]))
  th = c(th,B%*%log(piv))
  if(mod==0) for(t in 2:TT) for(u in 1:k) th = c(th,C[,,u]%*%log(Pi[u,,t]))
  if(mod==1) for(u in 1:k) th = c(th,C[,,u]%*%log(Pi[u,,2]))
  if(mod>1) {
    for(u in 1:k) th = c(th,C[,,u]%*%log(Pi[u,,2]))
    for(u in 1:k) th = c(th,C[,,u]%*%log(Pi[u,,mod+1]))
  }
  lth = length(th)
  out = recursions(S,R,yv,Psi,piv,Pi,k,lth,m,Bm,Cm,bv,mod)
  F1 = out$F1; F2 = out$F2; F1d = out$F1d; F2d = out$F2d

  sc = NULL
  Y = array(NA,c(b+1,TT,k,r))
  for(j in 1:r) Y[1:(bv[j]+1),,,j] = 0
  for(j in 1:r) for(t in 1:TT) for(jb in 0:bv[j]){
    ind = which(S[,t,j]==jb)
    if(length(ind)==1){
      if(miss) Y[jb+1,t,,j] = F1[,t,ind]*R[ind,t,j]*yv[ind]
      else Y[jb+1,t,,j] = F1[,t,ind]*yv[ind]
    }
    if(length(ind)>1){
      if(miss) Y[jb+1,t,,j] = F1[,t,ind]%*%(R[ind,t,j]*yv[ind])
      else Y[jb+1,t,,j] = F1[,t,ind]%*%yv[ind]
    }
  }
  Y1 = apply(Y,c(1,3,4),sum)
  for(u in 1:k) for(j in 1:r){
    indj = 1:(bv[j]+1)
    sc = c(sc,t(m[[j]]$Am)%*%(Y1[indj,u,j]-sum(Y1[indj,u,j])*Psi[indj,u,j]))
  }
  Y2 = Y1

  for(u in 1:k) for(j in 1:r){
    indj = 1:(bv[j]+1)
    Juj = sum(Y1[indj,u,j])*t(m[[j]]$Am)%*%(diag(Psi[indj,u,j])-Psi[indj,u,j]%o%Psi[indj,u,j])%*%m[[j]]$Am
    if(u==1 && j==1) J = Juj
    else J = blkdiag(J,Juj)
  }
  bv1 = F1[,1,]%*%yv
  sc = c(sc,t(Bm)%*%(bv1-sum(bv1)*piv))
  J = blkdiag(J,n*t(Bm)%*%(diag(piv)-piv%o%piv)%*%Bm)
  if(mod==0){
    for(t in 2:TT){
      Ut = 0
      for(i in 1:ns) Ut = Ut+yv[i]*F2[,,t,i]
      for(u in 1:k){
        sc = c(sc,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,t]))
        J = blkdiag(J,sum(Ut[u,])*t(Cm[,,u])%*%(diag(Pi[u,,t])-Pi[u,,t]%o%Pi[u,,t])%*%Cm[,,u])
      }
    }
  }
  if(mod==1){
    Ut = 0
    for(i in 1:ns) for(t in 2:TT) Ut = Ut+yv[i]*F2[,,t,i]
    for(u in 1:k){
      sc = c(sc,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,2]))
      J = blkdiag(J,sum(Ut[u,])*t(Cm[,,u])%*%(diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2])%*%Cm[,,u])
    }
  }
  if(mod>1){
    Ut=0
    for(i in 1:ns) for(t in 2:mod) Ut = Ut+yv[i]*F2[,,t,i]
    for(u in 1:k){
      sc = c(sc,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,2]))
      J =  blkdiag(J,sum(Ut[u,])*t(Cm[,,u])%*%(diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2])%*%Cm[,,u])
    }
    Ut=0
    for(i in 1:ns) for(t in (mod+1):TT) Ut = Ut+yv[i]*F2[,,t,i]
    for(u in 1:k){
      sc = c(sc,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,mod+1]))
      J = blkdiag(J,sum(Ut[u,])*t(Cm[,,u])%*%(diag(Pi[u,,mod+1])-Pi[u,,mod+1]%o%Pi[u,,mod+1])%*%Cm[,,u])
    }

  }
  J = as.matrix(J)
  Jd = NULL
  for(pa in 1:lth){
    scj = NULL
    Y = array(NA,c(b+1,TT,k,r))
    for(j in 1:r) Y[1:(bv[j]+1),,,j] = 0
    for(j in 1:r) for(t in 1:TT) for(jb in 0:b){
      ind = which(S[,t,j]==jb)
      if(length(ind)==1){
        if(miss) Y[jb+1,t,,j] = F1d[,t,ind,pa]*R[ind,t,j]*yv[ind]
        else Y[jb+1,t,,j] = F1d[,t,ind,pa]*yv[ind]
      }
      if(length(ind)>1){
        if(miss) Y[jb+1,t,,j] = F1d[,t,ind,pa]%*%(R[ind,t,j]*yv[ind])
        else Y[jb+1,t,,j] = F1d[,t,ind,pa]%*%yv[ind]
      }
    }
    Y1 = apply(Y,c(1,3,4),sum)
    for(u in 1:k) for(j in 1:r){
      indj = 1:(bv[j]+1)
      scj = c(scj,t(m[[j]]$Am)%*%(Y1[indj,u,j]-sum(Y1[indj,u,j])*Psi[indj,u,j]))
    }
    bv1 = F1d[,1,,pa]%*%yv
    scj = c(scj,t(Bm)%*%(bv1-sum(bv1)*piv))
    if(mod==0) {
      for(t in 2:TT){
        Ut = 0
        for(i in 1:ns) Ut = Ut+yv[i]*F2d[,,t,i,pa]
        for(u in 1:k) scj = c(scj,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,t]))
      }
    }
    if(mod==1){
      Ut = 0
      for(i in 1:ns) for(t in 2:TT) Ut = Ut+yv[i]*F2d[,,t,i,pa]
      for(u in 1:k) scj = c(scj,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,2]))
    }
    if(mod>1){
      Ut = 0
      for(i in 1:ns) for(t in 2:mod) Ut = Ut+yv[i]*F2d[,,t,i,pa]
      for(u in 1:k) scj = c(scj,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,2]))
      Ut = 0
      for(i in 1:ns) for(t in (mod+1):TT) Ut = Ut+yv[i]*F2d[,,t,i,pa]
      for(u in 1:k) scj = c(scj,t(Cm[,,u])%*%(Ut[u,]-sum(Ut[u,])*Pi[u,,mod+1]))
    }

    Jd = cbind(Jd,scj)
  }
  J = J-Jd
  Va = ginv(J)
  for(u in 1:k) for(j in 1:r){
    indj = 1:(bv[j]+1)
    Om = diag(Psi[indj,u,j])-Psi[indj,u,j]%o%Psi[indj,u,j]
    if(u==1) {
      if(j==1) M = Om%*%m[[j]]$Am
      else M = blkdiag(M,Om%*%m[[j]]$Am)
    } else M =  blkdiag(M,Om%*%m[[j]]$Am)
  }
  Om = diag(piv)-tcrossprod(piv,piv)
  M =  blkdiag(M,Om%*%Bm)
  if(mod==0){
    for(t in 2:TT) for(u in 1:k){
      Om = diag(Pi[u,,t])-Pi[u,,t]%o%Pi[u,,t]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
  }
  if(mod==1){
    for(u in 1:k){
      Om = diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
  }
  if(mod>1){
    for(u in 1:k){
      Om = diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
    for(u in 1:k){
      Om = diag(Pi[u,,mod+1])-Pi[u,,mod+1]%o%Pi[u,,mod+1]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
  }
  M = as.matrix(M)
  Va = M%*%Va%*%t(M)
  dVa = diag(Va)
  if(any(dVa<0)) warning("Negative elements in the estimated variance-covariance matrix for the parameters estimates")
  se = sqrt(abs(dVa))
  # Divide parameters
  nPsi = sum(bv+1)*k
  sePsi = se[1:nPsi]
  sepiv = se[nPsi+(1:k)]
  if(mod==0) sePi = se[nPsi+k+(1:(k*k*(TT-1)))]
  if(mod==1) sePi = se[nPsi+k+(1:(k*k))]
  if(mod>1) sePi = se[nPsi+k+(1:(k*k*2))]
  #to check derivatives
  # if(check_der){
  #   J0 = J
  #   th0 = th-10^-5/2
  #   out = lk_obs(th0,m,Bm,Cm,bv,k,S,R=R,yv,TT,r,mod)
  #   lk0 = out$lk; sc0 = out$sc
  #   lth = length(th)
  #   scn = rep(0,lth)
  #   J = matrix(0,lth,lth)
  #   for(j in 1:lth){
  #     thj = th0; thj[j] = thj[j]+10^-5
  #     out = lk_obs(thj,m,Bm,Cm,bv,k,S,R=R,yv,TT,r,mod)
  #     scn[j] = (out$lk-lk0)/10^-5
  #     J[,j] = (out$sc-sc0)/10^-5
  #   }
  #   J = -(J+t(J))/2
  #   print(c(lk,lk0))
  #   print(round(cbind(sc,scn,sc0),5))
  #   print(round(cbind(diag(J),diag(J0),diag(J)-diag(J0)),4))
  #   browser()
  # }


  sePsi0 = sePsi
  sePsi = array(NA,c(b+1,k,r))
  ind = 0
  for(u in 1:k) for(j in 1:r){
    indj = 1:(bv[j]+1)
    ind = max(ind)+indj
    sePsi[indj,u,j] = sePsi0[ind]
  }
  sePsi = array(sePsi,c(b+1,k,r))
  sePi0 = sePi
  sePi = array(0,c(k,k,TT))
  if(mod>1){
    sePi0 = array(sePi0,c(k,k,2))
    sePi0 = aperm(sePi0,c(2,1,3))
    sePi[,,2:mod] = sePi0[,,1]
    sePi[,,(mod+1):TT] = sePi0[,,2]
  } else {
    sePi[,,2:TT] = sePi0
    sePi = aperm(sePi,c(2,1,3))
  }
  dimnames(sePsi) = list(category=0:b,state=1:k,item=1:r)
  dimnames(sePi) = list(state=1:k,state=1:k,time=1:TT)

  out = list(sepiv = sepiv, sePi = sePi,sePsi = sePsi)
  return(out)
}
se.LMlatent <- function(est,...){
  out_se = TRUE
  fort = TRUE
  check_der = FALSE

  if(!is.null(est$seBe))
  {
    out = list(sePsi = est$sePsi, seBe = est$seBe,seGa = est$seGa)
    return(out)
  }

  data <- est$data
  responsesFormula = attributes(est)$responsesFormula
  latentFormula = attributes(est)$latentFormula

  tv.which <- attributes(est)$whichtv
  id.which <- attributes(est)$whichid
  id <- attributes(est)$id
  tv <- attributes(est)$time
  k <- est$k
  param <- est$paramLatent

  data.new <- data[,-c(id.which,tv.which), drop = FALSE]

  if(is.null(responsesFormula))
  {
    Y <- data.new
    Xmanifest <- NULL
    Xinitial <- NULL
    Xtrans <- NULL
  }else{
    temp <-  getResponses(data = data.new,formula = responsesFormula)
    Y <- temp$Y
    Xmanifest <- temp$X
    Xinitial <- NULL
    Xtrans <- NULL
  }
  if(min(Y,na.rm=T)>0){
    Y <- apply(Y,2,function(x) x - min(x, na.rm = TRUE))
  }

  if(!is.null(latentFormula))
  {
    temp <-  getLatent(data = data.new,latent = latentFormula, responses = responsesFormula)
    Xinitial <- temp$Xinitial
    Xtrans <- temp$Xtrans
  }
  tmp <-  long2matrices.internal(Y = Y, id = id, time = tv, yv = NULL,
                                        Xinitial = Xinitial, Xmanifest = Xmanifest, Xtrans = Xtrans)
  model <- tmp$model
  X1 <- tmp$Xinitial
  X2 <- tmp$Xtrans
  S <- tmp$Y
  yv = tmp$freq

  miss = any(is.na(S))
  if(miss){
    R = 1 * (!is.na(S))
    Y[is.na(S)] = 0
  }else{
    R = NULL
  }

  n = sum(yv)
  sS = dim(S)
  ns = sS[1]
  TT = sS[2]
  nc1 = dim(X1)[2]
  nc2 = dim(X2)[3]

  if(length(sS)==2){
    r = 1
    if(is.matrix(S)) S = array(S,c(dim(S),1))
  }else
  {
    r = sS[3]
  }

  Sv = matrix(S,ns*TT,r)


  if(miss) Rv = matrix(R,ns*TT,r)

  if(r==1){
    if(is.matrix(S)) S = array(S,c(dim(S),1))
    b = max(S); mb = b; sb = b
    Co = cbind(-diag(b),diag(b))
    Ma = cbind(lower.tri(matrix(1,b,b), diag = TRUE),rep(0,b))
    Ma = rbind(Ma,1-Ma)
  }else{
    b = rep(0,r)
    for(j in 1:r) b[j] = max(S[,,j])
    mb = max(b); sb = sum(b)
    Matr = vector("list",r)
    for(j in 1:r){
      Matr[[j]]$Co = cbind(-diag(b[j]),diag(b[j]))
      Ma = cbind(lower.tri(matrix(1,b[j],b[j]), diag = TRUE),rep(0,b[j]))
      Matr[[j]]$Ma = rbind(Ma,1-Ma)
    }
  }
  th = NULL; sc = NULL
  J = NULL

  # Covariate structure and related matrices: initial probabilities
  if(k == 2) GBe = as.matrix(c(0,1)) else{
    GBe = diag(k); GBe = GBe[,-1]
  }
  if(is.null(X1)){
    nc1=0
    Xlab = rep(1,ns)
    nameBe = NULL
  }else{
    if(is.vector(X1)) X1 = matrix(X1,ns,1)
    nc1 = dim(X1)[2] # number of covariates on the initial probabilities
    if(ns!= dim(X1)[1]) stop("dimension mismatch between S and X1")

    nameBe = colnames(X1)
    out =  aggr_data(X1,fort=fort)
    Xdis = out$data_dis
    if(nc1==1) Xdis = matrix(Xdis,length(Xdis),1)
    Xlab = out$label
  }
  Xndis = max(Xlab)
  XXdis = array(0,c(k,(k-1)*(nc1+1),Xndis))
  for(i in 1:Xndis){
    if(nc1==0) xdis = 1 else xdis = c(1,Xdis[i,])
    XXdis[,,i] = GBe%*%(diag(k-1)%x%t(xdis))
  }


  # for the transition probabilities
  if(is.null(X2) | any(dimnames(X2)[2] == "(Intercept)")){
    if(param=="difflogit"){
      warning("with X2=NULL parametrization difflogit not considered")
      param="multilogit"
    }
    nc2 = 0
    Zlab = rep(1,ns*(TT-1))
    nameGa = NULL
    Zndis = max(Zlab)
  }else{
    if(TT==2) X2 = array(X2,c(ns,1,dim(X2)[2]))
    if(is.matrix(X2)) X2 = array(X2,c(ns,TT-1,1))
    nc2 = dim(X2)[3] # number of covariates on the transition probabilities
    if(ns!= dim(X2)[1]) stop("dimension mismatch between S and X2")

    nameGa = colnames(aperm(X2,c(1,3,2)))
    Z = NULL
    for(t in 1:(TT-1)) Z = rbind(Z,X2[,t,])
    if(nc2==1) Z = as.vector(X2)
    out =  aggr_data(Z,fort=fort); Zdis = out$data_dis; Zlab = out$label; Zndis = max(Zlab)
    if(nc2==1) Zdis=matrix(Zdis,length(Zdis),1)
  }
  if(param=="multilogit"){
    ZZdis = array(0,c(k,(k-1)*(nc2+1),Zndis,k))
    for(h in 1:k){
      if(k==2){
        if(h == 1) GGa = as.matrix(c(0,1)) else GGa = as.matrix(c(1,0))
      }else{
        GGa = diag(k); GGa = GGa[,-h]
      }
      for(i in 1:Zndis){
        if(nc2==0) zdis = 1 else zdis = c(1,Zdis[i,])
        ZZdis[,,i,h] = GGa%*%(diag(k-1)%x%t(zdis))
      }
    }
  }else if(param=="difflogit"){
    Zlab = (((Zlab-1)*k)%x%rep(1,k))+rep(1,ns*(TT-1))%x%(1:k)
    ZZdis = array(0,c(k,k*(k-1)+(k-1)*nc2,Zndis*k))
    j = 0
    for(i in 1:Zndis){
      for(h in 1:k){
        j = j+1
        if(k==2){
          if(h == 1) GGa = as.matrix(c(0,1)) else GGa = as.matrix(c(1,0))
        }else{
          GGa = diag(k); GGa = GGa[,-h]
        }
        u = matrix(0,1,k); u[1,h] = 1
        U = diag(k); U[,h] = U[,h]-1
        U = U[,-1]
        ZZdis[,,j] = cbind(u%x%GGa,U%x%t(Zdis[i,]))
      }
    }
  }
  Piv = est$Piv
  PI = est$PI
  Psi = est$Psi
  Be = est$Be
  Ga = est$Ga
  V = est$V

  Am = vector("list",r)
  for(j in 1:r) Am[[j]] = rbind(rep(0,b[j]),diag(b[j]))

  be = as.vector(Be)
  out =  prob_multilogit(XXdis,be,Xlab,fort)
  Piv = out$P; Pivdis = out$Pdis
  # parameters on transition probabilities
  if(param=="multilogit"){
    if(is.list(Ga)) stop("invalid mode (list) for Ga")
    Ga = matrix(Ga,(nc2+1)*(k-1),k)
    PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
    for(h in 1:k){
      out =  prob_multilogit(ZZdis[,,,h],Ga[,h],Zlab,fort)
      PIdis[,,h] = out$Pdis; PI[h,,,2:TT] = array(as.vector(t(out$P)),c(1,k,ns,TT-1))
    }
  }else if(param=="difflogit"){
    if(is.list(Ga)) Ga = c(as.vector(t(Ga[[1]])),as.vector(Ga[[2]]))
    if(length(Ga)!=k*(k-1)+(k-1)*nc2) stop("invalid dimensions for Ga")
    PI = array(0,c(k,k,ns,TT))
    out =  prob_multilogit(ZZdis,Ga,Zlab,fort)
    PIdis = out$Pdis; PI[,,,2:TT] = array(as.vector(t(out$P)),c(k,k,ns,TT-1))
    PI = aperm(PI,c(2,1,3,4))
  }

  dlPsi = array(NA,c(mb+1,k,r,k*sb))
  for(j in 1:r) dlPsi[1:(b[j]+1),,j,] = 0
  count = 0
  for(c in 1:k) for(j in 1:r){
    ind = count+(1:b[j])
    temp = pmax(Psi[1:(b[j]+1),c,j],10^-50)
    dlPsi[1:(b[j]+1),c,j,ind] = (diag(b[j]+1)-rep(1,b[j]+1)%o%temp)%*%Am[[j]]
    count = count+b[j]
    th = c(th,log(temp[-1]/temp[1]))
  }
  dlPiv = array(0,c(ns,k,(1+nc1)*(k-1)))
  be = as.vector(Be)
  out =  est_multilogit(V[,,1],XXdis,Xlab,be,Pivdis,fort=fort)
  Pivdis = out$Pdi
  for(j in 1:Xndis){
    temp = pmax(Pivdis[j,],10^-50)
    Temp = (diag(k)-rep(1,k)%o%temp)%*%XXdis[,,j]
    for(i in which(Xlab==j)) dlPiv[i,,] = Temp
  }
  th = c(th,be)

  count = 0
  if(param=="multilogit"){
    dlPI = array(0,c(k,k,ns*TT,(1+nc2)*(k-1)*k))
    temp0 = rep(1,k); Temp0 = diag(k)
    for(h in 1:k){
      ind = count+(1:((1+nc2)*(k-1)))
      for(j in 1:Zndis){
        temp = pmax(PIdis[j,,h],10^-50)
        Temp = (Temp0-temp0%o%temp)%*%ZZdis[,,j,h]
        for(i in which(Zlab==j)) dlPI[h,,ns+i,ind] = Temp
      }
      count = count+((1+nc2)*(k-1))
      th = c(th,Ga[,h])
    }
    dlPI = array(dlPI,c(k,k,ns,TT,(1+nc2)*(k-1)*k))
  }else if(param=="difflogit"){
    dlPI = array(0,c(k,k*ns*TT,(k+nc2)*(k-1)))
    temp0 = rep(1,k); Temp0 = diag(k)
    for(j in 1:(Zndis*k)){
      temp = pmax(PIdis[j,],10^-50)
      Temp = (Temp0-temp0%o%temp)%*%ZZdis[,,j]
      for(i in which(Zlab==j)) dlPI[,k*ns+i,] = Temp
    }
    th = c(th,Ga)
    dlPI = array(dlPI,c(k,k,ns,TT,(k+nc2)*(k-1)))
    dlPI = aperm(dlPI,c(2,1,3,4,5))
  }

  # Compute log-likelihood
  lk = est$lk
  lk2 = lk
  out =  lk_comp_latent(S,R,yv,Piv,PI,Psi,k,der=TRUE,fort=fort,dlPsi=dlPsi,dlPiv=dlPiv,dlPI=dlPI)
  sc = out$dlk; dlL = out$dlL; dlPhi = out$dlPhi; dlL2 = out$dlL2; dlpv = out$dlpv
  lk = out$lk; Phi = out$Phi; L = out$L; pv = out$pv
  sc2 = sc;
  it = 0; lko = lk-10^10; lkv = NULL; dev = NULL
  # backward recursion
  out =  prob_post_cov(S,yv,Psi,Piv,PI,Phi,L,pv,der=TRUE,fort=fort,dlPhi=dlPhi,dlPiv,
                              dlPI=dlPI,dlL=dlL,dlL2=dlL2,dlpv=dlpv)
  U = out$U; V = out$V; dlU = out$dlU; dlV = out$dlV
  # ---- M-step ----
  # score and info Psi
  sc = NULL
  Y1 = array(NA,c(mb+1,k,r))
  for(j in 1:r) Y1[1:(b[j]+1),,j] = 0
  Vv = matrix(aperm(V,c(1,3,2)),ns*TT,k)
  for(j in 1:r) for(jb in 0:b[j]) {
    ind = which(Sv[,j]==jb)
    if(length(ind)==1){
      if(miss) Y1[jb+1,,j] = Vv[ind,]*Rv[ind,j]
      else Y1[jb+1,,j] = Vv[ind,]
    }
    if(length(ind)>1){
      if(miss) Y1[jb+1,,j] = colSums(Vv[ind,]*Rv[ind,j])
      else Y1[jb+1,,j] = colSums(Vv[ind,])
    }

  }
  for(c in 1:k) for(j in 1:r){
    sc = c(sc,t(Am[[j]])%*%(Y1[1:(b[j]+1),c,j]-sum(Y1[1:(b[j]+1),c,j])*Psi[1:(b[j]+1),c,j]))
    tmp = Y1[1:(b[j]+1),c,j]
    tmp = pmax(tmp/sum(tmp),10^-10)
    Psi[1:(b[j]+1),c,j] = tmp/sum(tmp)

    temp = pmax(Psi[1:(b[j]+1),c,j],10^-50)
    Op = diag(temp)-temp%o%temp
    Temp  = sum(Y1[1:(b[j]+1),c,j])*t(Am[[j]])%*%Op%*%Am[[j]]
    if(j==1 & c==1) Fi = Temp else Fi =  blkdiag(Fi,Temp)
  }

  # score and info piv
  out =  est_multilogit(V[,,1],XXdis,Xlab,be,Pivdis,fort=fort,ex=TRUE)
  sc = c(sc,out$sc); Fi =  blkdiag(Fi,out$Fi)
  # score and info Pi
  if(param=="multilogit"){
    for(h in 1:k){
      UU = NULL
      for(t in 2:TT) UU = rbind(UU,t(U[h,,,t]))
      tmp = ZZdis[,,,h]
      if(nc2==0) tmp = array(tmp,c(k,(k-1),Zndis))
      tmp2 = PIdis[,,h]
      if(Zndis==1) tmp2 = matrix(tmp2,1,k)
      out =  est_multilogit(UU,tmp,Zlab,Ga[,h],tmp2,fort=fort,ex=TRUE)
      sc = c(sc,out$sc); Fi =  blkdiag(Fi,out$Fi)
    }
  }else if(param=="difflogit"){
    Tmp = aperm(U[,,,2:TT],c(1,3,4,2))
    Tmp = matrix(Tmp,ns*k*(TT-1),k)
    out =  est_multilogit(Tmp,ZZdis,Zlab,Ga,PIdis,fort=fort,ex=TRUE)
    sc = c(sc,out$sc); Fi =  blkdiag(Fi,out$Fi)
  }
  Fi = as.matrix(Fi)
  # compute correction matrix for the information
  nal = dim(dlPhi)[4]; nbe = dim(dlPiv)[3]; nga = dim(dlPI)[5]
  npar = nal+nbe+nga
  Cor = matrix(0,npar,npar)
  dY1 = array(NA,c(mb+1,k,r,npar))
  for(j in 1:r) dY1[1:(b[j]+1),,j,] = 0
  dV = array(V,c(ns,k,TT,npar))*dlV
  dVv = array(aperm(dV,c(1,3,2,4)),c(ns*TT,k,nal+nbe+nga))
  for(j in 1:r) for(jb in 0:b[j]) {
    ind = which(Sv[,j]==jb)
    if(length(ind)==1){
      if(miss) for(h in 1:(nal+nbe+nga)) dY1[jb+1,,j,h] = dVv[ind,,h]*Rv[ind,j]
      else for(h in 1:(nal+nbe+nga)) dY1[jb+1,,j,h] = dVv[ind,,h]
    }
    if(length(ind)>1){
      if(miss) for(h in 1:(nal+nbe+nga)) dY1[jb+1,,j,h] = colSums(dVv[ind,,h]*Rv[ind,j])
      else for(h in 1:(nal+nbe+nga)) dY1[jb+1,,j,h] = colSums(dVv[ind,,h])
    }
  }
  for(h in 1:(npar)){
    count = 0
    for(c in 1:k) for(j in 1:r){
      ind = count+(1:b[j])
      count = count+b[j]
      Cor[h,ind] = t(Am[[j]])%*%(dY1[1:(b[j]+1),c,j,h]-sum(dY1[1:(b[j]+1),c,j,h])*Psi[1:(b[j]+1),c,j])
    }
  }
  for(h in 1:(npar)){
    out =  est_multilogit(dV[,,1,h],XXdis,Xlab,be,Pivdis,fort=fort,ex=TRUE)
    Cor[h,nal+(1:nbe)] = out$sc
  }
  dU = array(U,c(k,k,ns,TT,npar))*dlU
  if(param=="multilogit"){
    rGa = dim(Ga)[1]
    for(h in 1:k){
      for(h1 in 1:npar){
        UU = NULL
        for(t in 2:TT) UU = rbind(UU,t(dU[h,,,t,h1]))
        tmp = ZZdis[,,,h]
        if(nc2==0) tmp = array(tmp,c(k,(k-1),Zndis))
        tmp2 = PIdis[,,h]
        if(Zndis==1) tmp2 = matrix(tmp2,1,k)

        out =  est_multilogit(UU,tmp,Zlab,Ga[,h],tmp2,fort=fort,ex=TRUE)
        ind = nal+nbe+(h-1)*rGa+(1:rGa)
        Cor[h1,ind] = out$sc
      }
    }
  }else if(param=="difflogit"){
    rGa = length(Ga)
    for(h1 in 1:npar){
      Tmp = aperm(dU[,,,2:TT,h1],c(1,3,4,2))
      Tmp = matrix(Tmp,ns*k*(TT-1),k)
      out =  est_multilogit(Tmp,ZZdis,Zlab,Ga,PIdis,fort=fort,ex=TRUE)
      ind = nal+nbe+(1:rGa)
      Cor[h1,ind] = out$sc
    }
  }
  # check score and information
  if(check_der){
    lk0 = lk
    out = lk_obs_latent(th,S,R,b,yv,Am,XXdis,Xlab,ZZdis,Zlab,param,fort)
    print(lk-out$lk)
    sc1 = out$sc
    lth = length(th)
    scn = rep(0,lth); Fn = matrix(0,lth,lth)
    for(h in 1:lth){
      thh = th; thh[h] = thh[h]+10^-6
      outh = lk_obs_latent(thh,S,R,b,yv,Am,XXdis,Xlab,ZZdis,Zlab,param,fort=fort)
      scn[h] = (outh$lk-lk)*10^6
      Fn[,h] = -(outh$sc-sc)*10^6
    }
    print(round(cbind(sc,sc1,scn,sc-scn),4))
    print(round(cbind(diag(Fi-Cor),diag(Fn),diag(Fi-Cor-Fn)),4))
    browser()
  }
  # Information matrix and standard errors
  Fi = Fi-Cor
  iFi = ginv(Fi)
  se = sqrt(diag(iFi))
  # Divide parameters
  sepsi = se[1:nal]
  sebe = se[nal+(1:nbe)]
  sega = se[nal+nbe+(1:nga)]


  if(r==1){
    psi = as.vector(aperm(Psi,c(1,3,2)))
    dPsi = diag(psi)%*%matrix(aperm(dlPsi,c(1,3,2,4)),c((b+1)*k,nal))
    sePsi = sqrt(diag(dPsi%*%iFi[1:nal,1:nal]%*%t(dPsi)))
    sePsi = aperm(array(sePsi,c(b+1,r,c)),c(1,3,2))
    dimnames(sePsi)=list(category=0:b,state=1:k)
  }else{
    sePsi = array(NA,dim(Psi))
    for(j in 1:r) sePsi[1:(b[j]+1)]=0
    ind = 0
    for(c in 1:k) for(j in 1:r){
      Tmp = iFi[ind+(1:b[j]),ind+(1:b[j])]
      ind = ind+b[j]
      psi = Psi[1:(b[j]+1),c,j]
      Tmp1 = matrix((diag(psi)-psi%o%psi)[,-1],b[j]+1,b[j])
      sePsi[1:(b[j]+1),c,j] = sqrt(diag(Tmp1%*%Tmp%*%t(Tmp1)))
      dimnames(sePsi)=list(category=0:mb,state=1:k,item=1:r)
    }
  }

  if(out_se) {seBe = matrix(sebe,nc1+1,k-1); dimnames(seBe) = list(nameBe,logit=2:k)}
  if(param=="multilogit"){
    if(is.null(nameGa)){
      if(nc2==0) nameGa = c("Intercept") else nameGa = c("intercept", paste("X2",1:nc2,sep=""))
    }else{
      nameGa = c("intercept",nameGa)
    }
    if(k>2) {
      Ga = array(as.vector(Ga),c(nc2+1,k-1,k))
      dimnames(Ga) = list(nameGa,logit=2:k,logit=1:k)
    }else if(k==2){

      dimnames(Ga) = 	list(nameGa,logit=1:k)
    }
    if(out_se){
      if(k==2){
        seGa = matrix(sega,nc2+1,2)
        dimnames(seGa) = list(nameGa,logit=1:k)
      }else if(k>2){
        seGa = array(as.vector(sega),c(nc2+1,k-1,k))
        dimnames(seGa) = list(nameGa,logit=2:k,logit=1:k)
      }
    }
  }else if(param=="difflogit"){
    Ga0 = Ga
    Ga = vector("list",2)
    seGa = vector("list",2)
    Ga[[1]] = t(matrix(Ga0[1:(k*(k-1))],k-1,k))
    Ga[[2]] = matrix(Ga0[(k*(k-1))+(1:((k-1)*nc2))],nc2,k-1)
    if(is.null(nameGa)){
      nameGa2 = paste("X2",1:nc2,sep="")
    }else{
      nameGa2 = nameGa
    }
    if (k==2) {
      dimnames(Ga[[1]]) = list(intercept=1:k,logit=k)
      dimnames(Ga[[2]])=list(nameGa2,logit=k)
    } else if (k>2){
      dimnames(Ga[[1]]) = list(intercept=1:k,logit=2:k)
      dimnames(Ga[[2]])=list(nameGa2,logit=2:k)
    }
    if(out_se){
      seGa[[1]] = t(matrix(sega[1:(k*(k-1))],k-1,k))
      seGa[[2]] = matrix(sega[(k*(k-1))+(1:((k-1)*nc2))],nc2,k-1)
      if(k==2){
        dimnames(seGa[[1]]) = list(intercept=1:k,logit=k)
        dimnames(seGa[[2]])=list(nameGa2,logit=k)
      }else if (k>2){
        dimnames(seGa[[1]]) = list(intercept=1:k,logit=2:k)
        dimnames(seGa[[2]])=list(nameGa2,logit=2:k)
      }
    }
  }

  out = list(sePsi = sePsi, seBe = seBe,seGa = seGa)
  return(out)
}
se.LMbasiccont <- function(est,...){



  if(!is.null(est$seMu))
  {
    out = list(sepiv = est$sepiv,sePi = est$sePi,seMu = est$seMu,seSi = est$seSi)
    return(out)
  }

  data <- est$data
  responsesFormula = attributes(est)$responsesFormula
  latentFormula = attributes(est)$latentFormula

  tv.which <- attributes(est)$whichtv
  id.which <- attributes(est)$whichid
  id <- attributes(est)$id
  tv <- attributes(est)$time
  k <- est$k
  mod <- est$modBasic

  data.new <- data[,-c(id.which,tv.which), drop = FALSE]

  if(is.null(responsesFormula))
  {
    Y <- data.new
    Xmanifest <- NULL
    Xinitial <- NULL
    Xtrans <- NULL
  }else{
    temp <-  getResponses(data = data.new,formula = responsesFormula)
    Y <- temp$Y
    Xmanifest <- temp$X
    Xinitial <- NULL
    Xtrans <- NULL
  }
  if(min(Y,na.rm=T)>0){
    Y <- apply(Y,2,function(x) x - min(x, na.rm = TRUE))
  }

  tmp <-  long2matrices.internal(Y = Y, id = id, time = tv, yv = NULL,
                                        Xinitial = Xinitial, Xmanifest = Xmanifest, Xtrans = Xtrans)
  #model <- tmp$model
  Y <- tmp$Y

  miss = any(is.na(Y))

  if(miss){
    Yv = cbind(1,Yv)
    pYv = prelim.mix(Yv,1)
    thhat = em.mix(prelim.mix(Yv,1))
    rngseed(1)
    Yv = as.matrix(imp.mix(pYv, da.mix(pYv,thhat,steps=100), Yv)[,-1])
    Y = array(Yv,c(n,TT,r))
    cat("Missing data in the dataset. imp.mix function (mix package) used for imputation.\n")
  }

  piv = est$piv
  Pi = est$Pi
  Mu = est$Mu
  Si = est$Si



  sY = dim(Y)
  n = sY[1]
  TT = sY[2]

  if(length(sY)==2){
    r = 1
    if(is.matrix(Y)) Y = array(Y,c(dim(Y),1))
  }else r = sY[3]

  Yv = matrix(Y,n*TT,r)

  th = NULL; sc = NULL; J = NULL

  B = cbind(-rep(1,k-1),diag(k-1))
  Bm = rbind(rep(0,k-1),diag(k-1))
  C = array(0,c(k-1,k,k))
  Cm = array(0,c(k,k-1,k))
  for(u in 1:k){
    C[,,u] = rbind(cbind(diag(u-1),-rep(1,u-1),matrix(0,u-1,k-u)),
                   cbind(matrix(0,k-u,u-1),-rep(1,k-u),diag(k-u)))
    Cm[,,u] = rbind(cbind(diag(u-1),matrix(0,u-1,k-u)),
                    rep(0,k-1),
                    cbind(matrix(0,k-u,u-1),diag(k-u)))
  }



  th = NULL
  th = c(th,as.vector(Mu))
  th = c(th,Si[upper.tri(Si,TRUE)])
  th = c(th,B%*%log(piv))
  if(mod==0) for(t in 2:TT) for(u in 1:k) th = c(th,C[,,u]%*%log(Pi[u,,t]))
  if(mod==1) for(u in 1:k) th = c(th,C[,,u]%*%log(Pi[u,,2]))

  th0 = th-10^-5/2
  out =  lk_obs_cont(th0,Bm,Cm,k,Y,TT,r,mod)
  lk0 = out$lk; sc0 = out$sc
  lth = length(th)
  scn = rep(0,lth)
  J = matrix(0,lth,lth)
  for(j in 1:lth){
    thj = th0; thj[j] = thj[j]+10^-5
    out =  lk_obs_cont(thj,Bm,Cm,k,Y,TT,r,mod)
    scn[j] = (out$lk-lk0)/10^-5
    J[,j] = (out$sc-sc0)/10^-5
  }
  J = -(J+t(J))/2
  Va = ginv(J)
  nMu = r*k
  nSi = r*(r+1)/2
  Va2 = Va[1:(nMu+nSi),1:(nMu+nSi)]
  se2 = sqrt(diag(Va2))

  Va = Va[-(1:(nMu+nSi)),-(1:(nMu+nSi))]
  Om = diag(piv)-tcrossprod(piv,piv)
  M = Om%*%Bm
  if(mod==0){
    for(t in 2:TT) for(u in 1:k){
      Om = diag(Pi[u,,t])-Pi[u,,t]%o%Pi[u,,t]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
  }
  if(mod==1){
    for(u in 1:k){
      Om = diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
  }
  if(mod>1){
    for(u in 1:k){
      Om = diag(Pi[u,,2])-Pi[u,,2]%o%Pi[u,,2]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
    for(u in 1:k){
      Om = diag(Pi[u,,mod+1])-Pi[u,,mod+1]%o%Pi[u,,mod+1]
      M =  blkdiag(M,Om%*%Cm[,,u])
    }
  }
  M = as.matrix(M)
  Va = M%*%Va%*%t(M)
  dVa = diag(Va)
  if(any(dVa<0)) warning("Negative elements in the estimated variance-covariance matrix for the parameters estimates")
  se = sqrt(abs(dVa))
  # Divide parameters
  se = c(se2,se)
  seMu = se[1:nMu]
  seSi = se[nMu+(1:nSi)]
  sepiv = se[nMu+nSi+(1:k)]

  if(mod==0) sePi = se[nMu+nSi+k+(1:(k*k*(TT-1)))]
  if(mod==1) sePi = se[nMu+nSi+k+(1:(k*k))]
  if(mod>1) sePi = se[nMu+nSi+k+(1:(k*k*2))]



  if(miss) out$Y = Y


  seMu = matrix(seMu,r,k)
  seSi2 = matrix(0,r,r)
  seSi2[upper.tri(seSi2,TRUE)]=seSi
  seSi2 = seSi2+t(seSi2-diag(diag(seSi2)))
  seSi = seSi2
  sePi0 = sePi
  sePi = array(0,c(k,k,TT))
  if(mod>1){
    sePi0 = array(sePi0,c(k,k,2))
    sePi0 = aperm(sePi0,c(2,1,3))
    sePi[,,2:mod] = sePi0[,,1]
    sePi[,,(mod+1):TT] = sePi0[,,2]
  } else {
    sePi[,,2:TT] = sePi0
    sePi = aperm(sePi,c(2,1,3))
  }
  dimnames(sePi) = list(state=1:k,state=1:k,time=1:TT)
  if(r==1) dimnames(seMu) = list(item=1,state=1:k) else dimnames(seMu)=list(item=1:r,state=1:k)

  out = list(sepiv = sepiv,sePi = sePi,seMu = seMu,seSi = seSi)

  return(out)


}
se.LMlatentcont <- function(est,...){
  out_se = TRUE
  #fort = TRUE
  #check_der = FALSE

  if(!is.null(est$seBe))
  {
    out = list(seMu = out$seMu,seSi = out$seSi, seBe = est$seBe,seGa = est$seGa)
    return(out)
  }

  data <- est$data
  responsesFormula = attributes(est)$responsesFormula
  latentFormula = attributes(est)$latentFormula

  tv.which <- attributes(est)$whichtv
  id.which <- attributes(est)$whichid
  id <- attributes(est)$id
  tv <- attributes(est)$time
  k <- est$k
  param <- est$paramLatent

  data.new <- data[,-c(id.which,tv.which), drop = FALSE]

  if(is.null(responsesFormula))
  {
    Y <- data.new
    Xmanifest <- NULL
    Xinitial <- NULL
    Xtrans <- NULL
  }else{
    temp <-  getResponses(data = data.new,formula = responsesFormula)
    Y <- temp$Y
    Xmanifest <- temp$X
    Xinitial <- NULL
    Xtrans <- NULL
  }
  if(min(Y,na.rm=T)>0){
    Y <- apply(Y,2,function(x) x - min(x, na.rm = TRUE))
  }

  if(!is.null(latentFormula))
  {
    temp <-  getLatent(data = data.new,latent = latentFormula, responses = responsesFormula)
    Xinitial <- temp$Xinitial
    Xtrans <- temp$Xtrans
  }
  tmp <-  long2matrices.internal(Y = Y, id = id, time = tv, yv = NULL,
                                        Xinitial = Xinitial, Xmanifest = Xmanifest, Xtrans = Xtrans)
  model <- tmp$model
  X1 <- tmp$Xinitial
  X2 <- tmp$Xtrans
  Y <- tmp$Y
  yv = tmp$freq


  miss = any(is.na(Y))

  if(miss){
    Yv = cbind(1,Yv)
    pYv = prelim.mix(Yv,1)
    thhat = em.mix(prelim.mix(Yv,1))
    rngseed(1)
    Yv = as.matrix(imp.mix(pYv, da.mix(pYv,thhat,steps=100), Yv)[,-1])
    Y = array(Yv,c(n,TT,r))
    cat("Missing data in the dataset. imp.mix function (mix package) used for imputation.\n")
  }

  Be = est$Be
  Ga = est$Ga
  Mu = est$Mu
  Si = est$Si



  sY = dim(Y)
  n = sY[1]
  TT = sY[2]

  if(length(sY)==2){
    r = 1
    if(is.matrix(Y)) Y = array(Y,c(dim(Y),1))
  }else r = sY[3]

  Yv = matrix(Y,n*TT,r)

  if(k == 2){
    GBe = as.matrix(c(0,1))
  }else{
    GBe = diag(k); GBe = GBe[,-1]
  }

  if(is.vector(X1)) X1 = matrix(X1,n,1)
  nc1 = dim(X1)[2] # number of covariates on the initial probabilities
  if(n!= dim(X1)[1]) stop("dimension mismatch between S and X1")
  nameBe = colnames(X1)
  out =  aggr_data(X1)
  Xdis = out$data_dis
  if(nc1==1) Xdis = matrix(Xdis,length(Xdis),1)
  Xlab = out$label

  Xndis = max(Xlab)
  XXdis = array(0,c(k,(k-1)*(nc1+1),Xndis))
  for(i in 1:Xndis){
    if(nc1==0) xdis = 1 else xdis = c(1,Xdis[i,])
    XXdis[,,i] = GBe%*%(diag(k-1)%x%t(xdis))
  }
  if(TT==2) X2 = array(X2,c(n,1,dim(X2)[2]))
  if(is.matrix(X2)) X2 = array(X2,c(n,TT-1,1))
  nc2 = dim(X2)[3] # number of covariates on the transition probabilities
  if(n!= dim(X2)[1]) stop("dimension mismatch between S and X2")
  nameGa = colnames(aperm(X2,c(1,3,2)))
  Z = NULL
  for(t in 1:(TT-1)) Z = rbind(Z,X2[,t,])
  if(nc2==1) Z = as.vector(X2)
  out =  aggr_data(Z); Zdis = out$data_dis; Zlab = out$label; Zndis = max(Zlab)
  if(nc2==1) Zdis=matrix(Zdis,length(Zdis),1)

  if(param=="multilogit"){
    ZZdis = array(0,c(k,(k-1)*(nc2+1),Zndis,k))
    for(h in 1:k){
      if(k==2){
        if(h == 1) GGa = as.matrix(c(0,1)) else GGa = as.matrix(c(1,0))
      }else{
        GGa = diag(k); GGa = GGa[,-h]
      }
      for(i in 1:Zndis){
        if(nc2==0) zdis = 1 else zdis = c(1,Zdis[i,])
        ZZdis[,,i,h] = GGa%*%(diag(k-1)%x%t(zdis))
      }
    }
  }else if(param=="difflogit"){
    Zlab = (((Zlab-1)*k)%x%rep(1,k))+rep(1,n*(TT-1))%x%(1:k)
    ZZdis = array(0,c(k,k*(k-1)+(k-1)*nc2,Zndis*k))
    j = 0
    for(i in 1:Zndis){
      for(h in 1:k){
        j = j+1
        if(k==2){
          if(h == 1) GGa = as.matrix(c(0,1)) else GGa = as.matrix(c(1,0))
        }else{
          GGa = diag(k); GGa = GGa[,-h]
        }
        u = matrix(0,1,k); u[1,h] = 1
        U = diag(k); U[,h] = U[,h]-1
        U = U[,-1]
        ZZdis[,,j] = cbind(u%x%GGa,U%x%t(Zdis[i,]))
      }
    }
  }

  be = as.vector(Be)
  out =  prob_multilogit(XXdis,be,Xlab)
  Piv = out$P; Pivdis = out$Pdis
  # parameters on transition probabilities
  if(param=="multilogit"){
    if(is.list(Ga)) stop("invalid mode (list) for Ga")
    Ga = matrix(Ga,(nc2+1)*(k-1),k)
    PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,n,TT))
    for(h in 1:k){
      out =  prob_multilogit(ZZdis[,,,h],Ga[,h],Zlab)
      PIdis[,,h] = out$Pdis; PI[h,,,2:TT] = array(as.vector(t(out$P)),c(1,k,n,TT-1))
    }
  }else if(param=="difflogit"){
    if(is.list(Ga)) Ga = c(as.vector(t(Ga[[1]])),as.vector(Ga[[2]]))
    if(length(Ga)!=k*(k-1)+(k-1)*nc2) stop("invalid dimensions for Ga")
    PI = array(0,c(k,k,n,TT))
    out =  prob_multilogit(ZZdis,Ga,Zlab)
    PIdis = out$Pdis; PI[,,,2:TT] = array(as.vector(t(out$P)),c(k,k,n,TT-1))
    PI = aperm(PI,c(2,1,3,4))
  }



  th = NULL
  th = c(th,as.vector(Mu))
  th = c(th,Si[upper.tri(Si,TRUE)])
  th = c(th, be)
  if(param=="multilogit"){
    for(h in 1:k) th = c(th, Ga[,h])
  }else if(param=="difflogit") th = c(th,Ga)

  out =  lk_obs_latent_cont(th,Y,yv,XXdis,Xlab,ZZdis,Zlab,param)
  lk0 = out$lk; sc0 = out$sc

  lth = length(th)
  scn = rep(0,lth); Fn = matrix(0,lth,lth)

  for(h in 1:lth){
    thh = th; thh[h] = thh[h]+10^-5
    outh =  lk_obs_latent_cont(thh,Y,yv,XXdis,Xlab,ZZdis,Zlab,param)
    scn[h] = (outh$lk-lk0)/10^-5
    Fn[,h] = (outh$sc-sc0)/10^-5
  }

  J = -(Fn+t(Fn))/2
  iJ = ginv(J)
  se = sqrt(diag(iJ))

  nMu = r*k
  nSi = r*(r+1)/2
  nbe = (1+nc1)*(k-1)
  if(param=="multilogit") nga=(1+nc2)*(k-1)*k else if(param=="difflogit") nga=(k+nc2)*(k-1)
  seMu = se[1:nMu]
  seSi = se[nMu+(1:nSi)]
  sebe = se[nMu+nSi+(1:nbe)]
  sega = se[nMu+nSi+nbe+(1:nga)]

  if (is.null(nameBe)){
    if(nc1==0) nameBe = c("Intercept") else nameBe = c("intercept",paste("X1",1:nc1,sep=""))
  }else{
    nameBe = c("intercept",nameBe)
  }
  seBe = matrix(sebe,nc1+1,k-1); dimnames(seBe) = list(nameBe,logit=2:k)


  if(param=="multilogit"){
    if(is.null(nameGa)){
      if(nc2==0) nameGa = c("Intercept") else nameGa = c("intercept", paste("X2",1:nc2,sep=""))
    }else{
      nameGa = c("intercept",nameGa)
    }
    if(k>2) {
      Ga = array(as.vector(Ga),c(nc2+1,k-1,k))
      dimnames(Ga) = list(nameGa,logit=2:k,logit=1:k)
    }else if(k==2){

      dimnames(Ga) = 	list(nameGa,logit=1:k)
    }
    if(out_se){
      if(k==2){
        seGa = matrix(sega,nc2+1,2)
        dimnames(seGa) = list(nameGa,logit=1:k)
      }else if(k>2){
        seGa = array(as.vector(sega),c(nc2+1,k-1,k))
        dimnames(seGa) = list(nameGa,logit=2:k,logit=1:k)
      }
    }
  }else if(param=="difflogit"){
    Ga0 = Ga
    Ga = vector("list",2)
    seGa = vector("list",2)
    Ga[[1]] = t(matrix(Ga0[1:(k*(k-1))],k-1,k))
    Ga[[2]] = matrix(Ga0[(k*(k-1))+(1:((k-1)*nc2))],nc2,k-1)
    if(is.null(nameGa)){
      nameGa2 = paste("X2",1:nc2,sep="")
    }else{
      nameGa2 = nameGa
    }
    if (k==2) {
      dimnames(Ga[[1]]) = list(intercept=1:k,logit=k)
      dimnames(Ga[[2]])=list(nameGa2,logit=k)
    } else if (k>2){
      dimnames(Ga[[1]]) = list(intercept=1:k,logit=2:k)
      dimnames(Ga[[2]])=list(nameGa2,logit=2:k)
    }
    if(out_se){
      seGa[[1]] = t(matrix(sega[1:(k*(k-1))],k-1,k))
      seGa[[2]] = matrix(sega[(k*(k-1))+(1:((k-1)*nc2))],nc2,k-1)
      if(k==2){
        dimnames(seGa[[1]]) = list(intercept=1:k,logit=k)
        dimnames(seGa[[2]])=list(nameGa2,logit=k)
      }else if (k>2){
        dimnames(seGa[[1]]) = list(intercept=1:k,logit=2:k)
        dimnames(seGa[[2]])=list(nameGa2,logit=2:k)
      }
    }
  }

  seMu = matrix(seMu,r,k)
  if(r==1) dimnames(seMu) = list(item=1,state=1:k) else dimnames(seMu)=list(item=1:r,state=1:k)
  seSi2 = matrix(0,r,r)
  seSi2[upper.tri(seSi2,TRUE)]=seSi
  seSi2 = seSi2+t(seSi2-diag(diag(seSi2)))
  seSi = seSi2
  dimnames(seSi)=list(item=1:r,item=1:r)

  out <- list(seMu = seMu,seSi = seSi,seBe = seBe,seGa = seGa)
  return(out)

}

Try the LMest package in your browser

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

LMest documentation built on Aug. 27, 2023, 5:06 p.m.