Nothing
R/lmcovlatent.R
In LMest: Generalized Latent Markov Models
Defines functions
lmcovlatent
lmcovlatent <- function(S,X1=NULL,X2=NULL,yv=rep(1,nrow(S)),k,start=0,tol=10^-8,maxit=1000,
paramLatent="multilogit",Psi,Be,Ga,fort=TRUE,output=FALSE,out_se=FALSE,
fixPsi=FALSE,miss = FALSE, R = NULL, ntry = 0){
# Fit the LM model with individual covariates in the distribution of the latent process
#
# INPUT:
# S = array of available configurations (n x TT x r)
# X1 = matrix of covariates affecting the initial probabilities
# X2 = array of covariates affecting the transition probabilities
# Psi = conditional response probabilities
# Be = parameters on the initial probabilities (if start=2)
# Ga = parameters on the transition probabilities (if start=2)
# start = initialization (0 = deterministic, 1 = random, 2 = initial values in input)
# param = type of parametrization for the transition probabilities:
# multilogit = standard multinomial logit for every row of the transition matrix
# difflogit = multinomial logit based on the difference between two sets of parameters
# fort = fortran use (FALSE for not use fortran)
# output = to return additional output
# out_se = TRUE for computing the information and standard errors
# fixPsi = TRUE if Psi is given in input and is not updated anymore
# ---- Repeat estimation if necessary ----
if(ntry>0){
cat("* Deterministic inditialization *\n")
out = lmcovlatent(S,X1,X2,yv,k,start=0,tol,maxit,paramLatent,fort=fort,output=output,
out_se=out_se,miss=miss,R=R)
lkv_glob = out$lk
for(it0 in 1:(k*ntry)){
cat("\n* Random initialization (",it0,"/",k*ntry,") *\n",sep="")
outr = try(lmcovlatent(S,X1,X2,yv,k,start=1,tol,maxit,paramLatent,fort=fort,output=output,
out_se=out_se,miss=miss,R=R))
if(!inherits(outr,"try-error")){
lkv_glob = c(lkv_glob,outr$lk)
out = outr
}
}
out$lkv_glob = lkv_glob
return(out)
}
#---- Preliminaries ----
if(paramLatent=="difflogit"){
cat("\n* With difflogit is not possible to avoid the intercept for the transition probabilities*\n\n")
X2 = X2[,,-1,drop=FALSE]
}
check_der = FALSE # to check score and info
param <- paramLatent
if(fort!=TRUE) fort = FALSE
sS = dim(S)
ns = sS[1]
TT = sS[2]
n = sum(yv)
if(is.null(X1)){
X1 = matrix(1,ns,1)
colnames(X1) = "(Intercept)"
}
if(is.null(X2)){
X2 = array(1,c(ns,TT-1,1))
dimnames(X2) = list(NULL,NULL,"(Intercept)")
}
if(length(sS)==2) r = 1
else r = sS[3]
if(ns!=length(yv)) stop("dimensions mismatch between S and yv")
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] # dimension of X1
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,Xndis))
for(i in 1:Xndis){
if(nc1==0){
xdis = 1}
else{
if(is.matrix(Xdis)) xdis = Xdis[i,] else xdis = Xdis
}
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] # dimension of X2
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),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 = 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,]))
}
}
}
# for information matrix
if(out_se){
Am = vector("list",r)
for(j in 1:r) Am[[j]] = rbind(rep(0,b[j]),diag(b[j]))
}
# When there is just 1 latent class
if(k == 1){
Piv = rep(1,n); Pi = 1
P = matrix(0,mb+1,r)
for(t in 1:TT){
for(j in 1:r){
for(y in 0:b[j]){
ind = which(S[,t,j]==y)
P[y+1,j] = P[y+1,j]+sum(yv[ind])
}
}
}
Psi = P/(n*TT)
pm = rep(1,ns)
if (miss) for(t in 1:TT) for(j in 1:r) pm = pm*(Psi[S[,t,j]+1,j]*R[,t,j]+(1-R[,t,j]))
else for(t in 1:TT) for(j in 1:r) pm = pm*Psi[S[,t,j]+1,j]
lk = sum(yv*log(pm))
if(r==1) np = k*mb*r else np = k*sum(b)
aic = -2*lk+np*2
bic = -2*lk+np*log(n)
out = list(lk=lk,Piv=Piv,Pi=Pi,Psi=Psi,np=np, k = k,aic=aic,bic=bic,lkv=NULL,V=NULL, n = n, TT = TT, paramLatent=param )
class(out)="LMlatent"
return(out)
}
time = proc.time()
# Starting values: deterministic initialization
if(start == 0){
if(fixPsi==FALSE){
P = matrix(NA,mb+1,r)
for(j in 1:r) P[1:(b[j]+1),j] = 0
for(t in 1:TT) for(j in 1:r) for(y in 0:mb){
ind = which(S[,t,j]==y)
if(miss) P[y+1,j] = P[y+1,j]+sum(t(R[ind,t,j])%*%yv[ind])
else P[y+1,j] = P[y+1,j]+sum(yv[ind])
}
if(r==1) E = Co%*%log(Ma%*%P) else{
E = matrix(NA,mb,r)
for(j in 1:r){
Co = Matr[[j]]$Co; Ma = Matr[[j]]$Ma
E[1:b[j],j] = Co%*%log(Ma%*%P[1:(b[j]+1),j])
}
}
Psi = array(NA,c(mb+1,k,r)); Eta = array(NA,c(mb,k,r))
grid = seq(-k,k,2*k/(k-1))
for(c in 1:k){
for(j in 1:r){
etac = E[1:b[j],j]+grid[c]
Eta[1:b[j],c,j] = etac
Psi[1:(b[j]+1),c,j] = invglob(etac)
}
}
}
# parameters on initial probabilities
be = rep(0,nc1*(k-1))
out = prob_multilogit(XXdis,be,Xlab,fort)
Piv = out$P; Pivdis = out$Pdis
# parameters on transition probabilities
if(param=="multilogit"){
Ga = matrix(0,nc2*(k-1),k)
Ga[1+(0:(k-2))*nc2,] = -log(10)
PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
for(h in 1:k){
tmp = ZZdis[,,,h]
if(nc2==1) tmp = array(tmp,c(k,(k-1),Zndis))
out = prob_multilogit(tmp,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"){
Ga = matrix(0,k*(k-1)+(k-1)*nc2)
Ga[1:((h-1)*k)] = -log(10)
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))
}
}
# random initialization
if(start==1){
if(fixPsi==FALSE){
Psi = array(NA,c(mb+1,k,r))
for(j in 1:r){
Psi[1:(b[j]+1),,j] = matrix(runif((b[j]+1)*k),b[j]+1,k)
for(c in 1:k) Psi[1:(b[j]+1),c,j] = Psi[1:(b[j]+1),c,j]/sum(Psi[1:(b[j]+1),c,j])
}
}
# parameters on initial probabilities
be = c(rnorm(1),rep(0,nc1-1))
if(k>2) for(h in 2:(k-1)) be = c(be,rnorm(1),rep(0,nc1-1))
out = prob_multilogit(XXdis,be,Xlab,fort)
Piv = out$P; Pivdis = out$Pdis
# parameters on transition probabilities
if(param=="multilogit"){
Ga = matrix(0,(nc2)*(k-1),k)
Ga[1+(0:(k-2))*(nc2),] = -abs(rnorm((k-1)))
PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
for(h in 1:k){
tmp = ZZdis[,,,h]
if(nc2==1) tmp = array(tmp,c(k,(k-1),Zndis))
out = prob_multilogit(tmp,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"){
Ga = c(-abs(rnorm(k*(k-1))),rep(0,(k-1)*nc2))
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))
}
}
# initialization as input
if(start==2){
if(is.null(Be)) stop("initial value of the parameters on the initial probabilities (Be) must be given in input")
if(is.null(Ga)) stop("initial value of parameters on the transition probabilities (Ga) must be given in input")
if(is.null(Psi)) stop("initial value of the conditional response probabilities (Psi) must be given in input")
# parameters on initial probabilities
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,max(1,nc2)*(k-1),k)
PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
for(h in 1:k){
tmp = ZZdis[,,,h]
if(nc2==1) tmp = array(tmp,c(k,(k-1),Zndis))
out = prob_multilogit(tmp,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))
}
}
###### standard EM #####
out = lk_comp_latent(S,R,yv,Piv,PI,Psi,k,fort=fort)
lk = out$lk; Phi = out$Phi; L = out$L; pv = out$pv
it = 0; lko = lk-10^10; lkv = NULL
par = c(as.vector(Piv),as.vector(PI),as.vector(Psi))
if(any(is.na(par))) par = par[-which(is.na(par))]
paro = par
# Iterate until convergence
# display output
cat("------------|-------------|-------------|-------------|-------------|-------------|\n");
cat(" k | start | step | lk | lk-lko | discrepancy |\n");
cat("------------|-------------|-------------|-------------|-------------|-------------|\n");
cat(sprintf("%11g",c(k,start,0,lk)),"\n",sep=" | ")
while((lk-lko)/abs(lk)>tol & it<maxit){
Psi0 = Psi; Piv0 = Piv; PI0 = PI
it = it+1
# ---- E-step ----
# Compute V and U
out = prob_post_cov(S,yv,Psi,Piv,PI,Phi,L,pv,fort=fort)
U = out$U; V = out$V
# If required store parameters
# ---- M-step ----
# Update Psi
if(fixPsi==FALSE){
Y1 = array(NA,c(mb+1,k,r))
for(j in 1:r) Y1[1:(b[j]+1)] = 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(j in 1:r) for(c in 1:k){
tmp = Y1[1:(b[j]+1),c,j]
if(any(is.na(tmp))) tmp[is.na(tmp)] = 0
tmp = pmax(tmp/sum(tmp),10^-10)
Psi[1:(b[j]+1),c,j] = tmp/sum(tmp)
}
}
# Update piv
out = est_multilogit(V[,,1],XXdis,Xlab,be,Pivdis,fort=fort)
be = out$be; Pivdis = out$Pdi; Piv = out$P
# Update Pi
if(param=="multilogit"){
for(h in 1:k){
UU = NULL
for(t in 2:TT) UU = rbind(UU,t(U[h,,,t]))
tmp = array(ZZdis[,,,h],dim(ZZdis)[1:3])
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)
PIdis[,,h] = out$Pdis; PI[h,,,2:TT] = array(as.vector(t(out$P)),c(1,k,ns,TT-1)); Ga[,h] = out$be
}
}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)
PIdis = out$Pdis; Ga = out$be
Tmp = array(out$P,c(k,ns,TT-1,k))
PI[,,,2:TT] = aperm(Tmp,c(1,4,2,3))
}
# Compute log-likelihood
paro = par; par = c(as.vector(Piv),as.vector(PI),as.vector(Psi))
if(any(is.na(par))) par = par[-which(is.na(par))]
lko = lk;
out = lk_comp_latent(S,R,yv,Piv,PI,Psi,k,fort=fort)
lk = out$lk; Phi = out$Phi; L = out$L; pv = out$pv
# Display output
if(it/10 == floor(it/10)){
cat(sprintf("%11g",c(k,start,it,lk,lk-lko,max(abs(par-paro)))),"\n",sep=" | ")
}
lkv = c(lkv,lk)
}
if(it/10 > floor(it/10)) cat(sprintf("%11g",c(k,start,it,lk,lk-lko,max(abs(par-paro)))),"\n",sep=" | ")
#### compute infomation matrix ####
if(out_se){
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,nc1*(k-1)))
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,nc2*(k-1)*k))
temp0 = rep(1,k); Temp0 = diag(k)
for(h in 1:k){
ind = count+(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+(nc2*(k-1))
th = c(th,Ga[,h])
}
dlPI = array(dlPI,c(k,k,ns,TT,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
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==1) 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==1) 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)]
}
# Compute number of parameters
if(r==1) np = k*mb*r else np = k*sum(b)
np = np+(k-1)*nc1
if(param=="multilogit") np = np+(k-1)*nc2*k else if(param=="difflogit") np = np+(k-1)*(nc2+k)
aic = -2*lk+np*2
bic = -2*lk+np*log(n)
# local decoding
Ul = matrix(0,ns,TT)
for(i in 1:ns) for(t in 1:TT) Ul[i,t] = which.max(V[i,,t])
if(all(yv==1)) V1=V else V1 = V/yv
if(out_se){
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)
}
}
}
Be = matrix(be,nc1,k-1)
if(is.null(nameBe)) nameBe = c("(Intercept)",paste("X1",1:(nc1-1),sep=""))
dimnames(Be) = list(nameBe,logit=2:k)
if(out_se){seBe = matrix(sebe,nc1,k-1); dimnames(seBe) = list(nameBe,logit=2:k)}
if(param=="multilogit"){
if(is.null(nameGa)){
nameGa = c("(Intercept)", paste("X2",1:(nc2-1),sep=""))
}
if(k>2) {
Ga = array(as.vector(Ga),c(nc2,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,2)
dimnames(seGa) = list(nameGa,logit=1:k)
}else if(k>2){
seGa = array(as.vector(sega),c(nc2,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)
}
}
}
# adjust output
lk = as.vector(lk)
dimnames(Piv)=list(subject=1:ns,state=1:k)
dimnames(PI)=list(state=1:k,state=1:k,subject=1:ns,time=1:TT)
if(r==1) dimnames(Psi) = list(category=0:b,state=1:k,item=1) else dimnames(Psi)=list(category=0:mb,state=1:k,item=1:r)
out = list(lk=lk,Be=Be,Ga=Ga,Psi=Psi,Piv = Piv, PI = PI, np=np,k = k, aic=aic,bic=bic,
lkv=lkv, n = n, TT = TT,paramLatent=param,ns=ns,yv=yv)
if(out_se){
out$sePsi = sePsi
out$seBe = seBe
out$seGa = seGa
}
# final output
if(output){
if(k>1){
PMarg <- array(0,c(ns,k,TT))
PMarg[,,1] <- as.matrix(Piv)
for(i in 1:ns) for(t in 2:TT) PMarg[i,,t]= t(PI[,,i,t])%*%PMarg[i,,t-1]
Pmarg <-apply(PMarg,c(2,3),mean)
}else Pmarg<- NULL
out = c(out,list(V = V1, Ul = Ul, S = S, Pmarg=Pmarg))
}
cat("------------|-------------|-------------|-------------|-------------|-------------|\n");
class(out)="LMlatent"
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 Sept. 16, 2024, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lmcovlatent
lmcovlatent <- function(S,X1=NULL,X2=NULL,yv=rep(1,nrow(S)),k,start=0,tol=10^-8,maxit=1000,
paramLatent="multilogit",Psi,Be,Ga,fort=TRUE,output=FALSE,out_se=FALSE,
fixPsi=FALSE,miss = FALSE, R = NULL, ntry = 0){
# Fit the LM model with individual covariates in the distribution of the latent process
#
# INPUT:
# S = array of available configurations (n x TT x r)
# X1 = matrix of covariates affecting the initial probabilities
# X2 = array of covariates affecting the transition probabilities
# Psi = conditional response probabilities
# Be = parameters on the initial probabilities (if start=2)
# Ga = parameters on the transition probabilities (if start=2)
# start = initialization (0 = deterministic, 1 = random, 2 = initial values in input)
# param = type of parametrization for the transition probabilities:
# multilogit = standard multinomial logit for every row of the transition matrix
# difflogit = multinomial logit based on the difference between two sets of parameters
# fort = fortran use (FALSE for not use fortran)
# output = to return additional output
# out_se = TRUE for computing the information and standard errors
# fixPsi = TRUE if Psi is given in input and is not updated anymore
# ---- Repeat estimation if necessary ----
if(ntry>0){
cat("* Deterministic inditialization *\n")
out = lmcovlatent(S,X1,X2,yv,k,start=0,tol,maxit,paramLatent,fort=fort,output=output,
out_se=out_se,miss=miss,R=R)
lkv_glob = out$lk
for(it0 in 1:(k*ntry)){
cat("\n* Random initialization (",it0,"/",k*ntry,") *\n",sep="")
outr = try(lmcovlatent(S,X1,X2,yv,k,start=1,tol,maxit,paramLatent,fort=fort,output=output,
out_se=out_se,miss=miss,R=R))
if(!inherits(outr,"try-error")){
lkv_glob = c(lkv_glob,outr$lk)
out = outr
}
}
out$lkv_glob = lkv_glob
return(out)
}
#---- Preliminaries ----
if(paramLatent=="difflogit"){
cat("\n* With difflogit is not possible to avoid the intercept for the transition probabilities*\n\n")
X2 = X2[,,-1,drop=FALSE]
}
check_der = FALSE # to check score and info
param <- paramLatent
if(fort!=TRUE) fort = FALSE
sS = dim(S)
ns = sS[1]
TT = sS[2]
n = sum(yv)
if(is.null(X1)){
X1 = matrix(1,ns,1)
colnames(X1) = "(Intercept)"
}
if(is.null(X2)){
X2 = array(1,c(ns,TT-1,1))
dimnames(X2) = list(NULL,NULL,"(Intercept)")
}
if(length(sS)==2) r = 1
else r = sS[3]
if(ns!=length(yv)) stop("dimensions mismatch between S and yv")
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] # dimension of X1
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,Xndis))
for(i in 1:Xndis){
if(nc1==0){
xdis = 1}
else{
if(is.matrix(Xdis)) xdis = Xdis[i,] else xdis = Xdis
}
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] # dimension of X2
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),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 = 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,]))
}
}
}
# for information matrix
if(out_se){
Am = vector("list",r)
for(j in 1:r) Am[[j]] = rbind(rep(0,b[j]),diag(b[j]))
}
# When there is just 1 latent class
if(k == 1){
Piv = rep(1,n); Pi = 1
P = matrix(0,mb+1,r)
for(t in 1:TT){
for(j in 1:r){
for(y in 0:b[j]){
ind = which(S[,t,j]==y)
P[y+1,j] = P[y+1,j]+sum(yv[ind])
}
}
}
Psi = P/(n*TT)
pm = rep(1,ns)
if (miss) for(t in 1:TT) for(j in 1:r) pm = pm*(Psi[S[,t,j]+1,j]*R[,t,j]+(1-R[,t,j]))
else for(t in 1:TT) for(j in 1:r) pm = pm*Psi[S[,t,j]+1,j]
lk = sum(yv*log(pm))
if(r==1) np = k*mb*r else np = k*sum(b)
aic = -2*lk+np*2
bic = -2*lk+np*log(n)
out = list(lk=lk,Piv=Piv,Pi=Pi,Psi=Psi,np=np, k = k,aic=aic,bic=bic,lkv=NULL,V=NULL, n = n, TT = TT, paramLatent=param )
class(out)="LMlatent"
return(out)
}
time = proc.time()
# Starting values: deterministic initialization
if(start == 0){
if(fixPsi==FALSE){
P = matrix(NA,mb+1,r)
for(j in 1:r) P[1:(b[j]+1),j] = 0
for(t in 1:TT) for(j in 1:r) for(y in 0:mb){
ind = which(S[,t,j]==y)
if(miss) P[y+1,j] = P[y+1,j]+sum(t(R[ind,t,j])%*%yv[ind])
else P[y+1,j] = P[y+1,j]+sum(yv[ind])
}
if(r==1) E = Co%*%log(Ma%*%P) else{
E = matrix(NA,mb,r)
for(j in 1:r){
Co = Matr[[j]]$Co; Ma = Matr[[j]]$Ma
E[1:b[j],j] = Co%*%log(Ma%*%P[1:(b[j]+1),j])
}
}
Psi = array(NA,c(mb+1,k,r)); Eta = array(NA,c(mb,k,r))
grid = seq(-k,k,2*k/(k-1))
for(c in 1:k){
for(j in 1:r){
etac = E[1:b[j],j]+grid[c]
Eta[1:b[j],c,j] = etac
Psi[1:(b[j]+1),c,j] = invglob(etac)
}
}
}
# parameters on initial probabilities
be = rep(0,nc1*(k-1))
out = prob_multilogit(XXdis,be,Xlab,fort)
Piv = out$P; Pivdis = out$Pdis
# parameters on transition probabilities
if(param=="multilogit"){
Ga = matrix(0,nc2*(k-1),k)
Ga[1+(0:(k-2))*nc2,] = -log(10)
PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
for(h in 1:k){
tmp = ZZdis[,,,h]
if(nc2==1) tmp = array(tmp,c(k,(k-1),Zndis))
out = prob_multilogit(tmp,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"){
Ga = matrix(0,k*(k-1)+(k-1)*nc2)
Ga[1:((h-1)*k)] = -log(10)
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))
}
}
# random initialization
if(start==1){
if(fixPsi==FALSE){
Psi = array(NA,c(mb+1,k,r))
for(j in 1:r){
Psi[1:(b[j]+1),,j] = matrix(runif((b[j]+1)*k),b[j]+1,k)
for(c in 1:k) Psi[1:(b[j]+1),c,j] = Psi[1:(b[j]+1),c,j]/sum(Psi[1:(b[j]+1),c,j])
}
}
# parameters on initial probabilities
be = c(rnorm(1),rep(0,nc1-1))
if(k>2) for(h in 2:(k-1)) be = c(be,rnorm(1),rep(0,nc1-1))
out = prob_multilogit(XXdis,be,Xlab,fort)
Piv = out$P; Pivdis = out$Pdis
# parameters on transition probabilities
if(param=="multilogit"){
Ga = matrix(0,(nc2)*(k-1),k)
Ga[1+(0:(k-2))*(nc2),] = -abs(rnorm((k-1)))
PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
for(h in 1:k){
tmp = ZZdis[,,,h]
if(nc2==1) tmp = array(tmp,c(k,(k-1),Zndis))
out = prob_multilogit(tmp,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"){
Ga = c(-abs(rnorm(k*(k-1))),rep(0,(k-1)*nc2))
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))
}
}
# initialization as input
if(start==2){
if(is.null(Be)) stop("initial value of the parameters on the initial probabilities (Be) must be given in input")
if(is.null(Ga)) stop("initial value of parameters on the transition probabilities (Ga) must be given in input")
if(is.null(Psi)) stop("initial value of the conditional response probabilities (Psi) must be given in input")
# parameters on initial probabilities
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,max(1,nc2)*(k-1),k)
PIdis = array(0,c(Zndis,k,k)); PI = array(0,c(k,k,ns,TT))
for(h in 1:k){
tmp = ZZdis[,,,h]
if(nc2==1) tmp = array(tmp,c(k,(k-1),Zndis))
out = prob_multilogit(tmp,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))
}
}
###### standard EM #####
out = lk_comp_latent(S,R,yv,Piv,PI,Psi,k,fort=fort)
lk = out$lk; Phi = out$Phi; L = out$L; pv = out$pv
it = 0; lko = lk-10^10; lkv = NULL
par = c(as.vector(Piv),as.vector(PI),as.vector(Psi))
if(any(is.na(par))) par = par[-which(is.na(par))]
paro = par
# Iterate until convergence
# display output
cat("------------|-------------|-------------|-------------|-------------|-------------|\n");
cat(" k | start | step | lk | lk-lko | discrepancy |\n");
cat("------------|-------------|-------------|-------------|-------------|-------------|\n");
cat(sprintf("%11g",c(k,start,0,lk)),"\n",sep=" | ")
while((lk-lko)/abs(lk)>tol & it<maxit){
Psi0 = Psi; Piv0 = Piv; PI0 = PI
it = it+1
# ---- E-step ----
# Compute V and U
out = prob_post_cov(S,yv,Psi,Piv,PI,Phi,L,pv,fort=fort)
U = out$U; V = out$V
# If required store parameters
# ---- M-step ----
# Update Psi
if(fixPsi==FALSE){
Y1 = array(NA,c(mb+1,k,r))
for(j in 1:r) Y1[1:(b[j]+1)] = 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(j in 1:r) for(c in 1:k){
tmp = Y1[1:(b[j]+1),c,j]
if(any(is.na(tmp))) tmp[is.na(tmp)] = 0
tmp = pmax(tmp/sum(tmp),10^-10)
Psi[1:(b[j]+1),c,j] = tmp/sum(tmp)
}
}
# Update piv
out = est_multilogit(V[,,1],XXdis,Xlab,be,Pivdis,fort=fort)
be = out$be; Pivdis = out$Pdi; Piv = out$P
# Update Pi
if(param=="multilogit"){
for(h in 1:k){
UU = NULL
for(t in 2:TT) UU = rbind(UU,t(U[h,,,t]))
tmp = array(ZZdis[,,,h],dim(ZZdis)[1:3])
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)
PIdis[,,h] = out$Pdis; PI[h,,,2:TT] = array(as.vector(t(out$P)),c(1,k,ns,TT-1)); Ga[,h] = out$be
}
}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)
PIdis = out$Pdis; Ga = out$be
Tmp = array(out$P,c(k,ns,TT-1,k))
PI[,,,2:TT] = aperm(Tmp,c(1,4,2,3))
}
# Compute log-likelihood
paro = par; par = c(as.vector(Piv),as.vector(PI),as.vector(Psi))
if(any(is.na(par))) par = par[-which(is.na(par))]
lko = lk;
out = lk_comp_latent(S,R,yv,Piv,PI,Psi,k,fort=fort)
lk = out$lk; Phi = out$Phi; L = out$L; pv = out$pv
# Display output
if(it/10 == floor(it/10)){
cat(sprintf("%11g",c(k,start,it,lk,lk-lko,max(abs(par-paro)))),"\n",sep=" | ")
}
lkv = c(lkv,lk)
}
if(it/10 > floor(it/10)) cat(sprintf("%11g",c(k,start,it,lk,lk-lko,max(abs(par-paro)))),"\n",sep=" | ")
#### compute infomation matrix ####
if(out_se){
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,nc1*(k-1)))
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,nc2*(k-1)*k))
temp0 = rep(1,k); Temp0 = diag(k)
for(h in 1:k){
ind = count+(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+(nc2*(k-1))
th = c(th,Ga[,h])
}
dlPI = array(dlPI,c(k,k,ns,TT,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
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==1) 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==1) 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)]
}
# Compute number of parameters
if(r==1) np = k*mb*r else np = k*sum(b)
np = np+(k-1)*nc1
if(param=="multilogit") np = np+(k-1)*nc2*k else if(param=="difflogit") np = np+(k-1)*(nc2+k)
aic = -2*lk+np*2
bic = -2*lk+np*log(n)
# local decoding
Ul = matrix(0,ns,TT)
for(i in 1:ns) for(t in 1:TT) Ul[i,t] = which.max(V[i,,t])
if(all(yv==1)) V1=V else V1 = V/yv
if(out_se){
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)
}
}
}
Be = matrix(be,nc1,k-1)
if(is.null(nameBe)) nameBe = c("(Intercept)",paste("X1",1:(nc1-1),sep=""))
dimnames(Be) = list(nameBe,logit=2:k)
if(out_se){seBe = matrix(sebe,nc1,k-1); dimnames(seBe) = list(nameBe,logit=2:k)}
if(param=="multilogit"){
if(is.null(nameGa)){
nameGa = c("(Intercept)", paste("X2",1:(nc2-1),sep=""))
}
if(k>2) {
Ga = array(as.vector(Ga),c(nc2,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,2)
dimnames(seGa) = list(nameGa,logit=1:k)
}else if(k>2){
seGa = array(as.vector(sega),c(nc2,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)
}
}
}
# adjust output
lk = as.vector(lk)
dimnames(Piv)=list(subject=1:ns,state=1:k)
dimnames(PI)=list(state=1:k,state=1:k,subject=1:ns,time=1:TT)
if(r==1) dimnames(Psi) = list(category=0:b,state=1:k,item=1) else dimnames(Psi)=list(category=0:mb,state=1:k,item=1:r)
out = list(lk=lk,Be=Be,Ga=Ga,Psi=Psi,Piv = Piv, PI = PI, np=np,k = k, aic=aic,bic=bic,
lkv=lkv, n = n, TT = TT,paramLatent=param,ns=ns,yv=yv)
if(out_se){
out$sePsi = sePsi
out$seBe = seBe
out$seGa = seGa
}
# final output
if(output){
if(k>1){
PMarg <- array(0,c(ns,k,TT))
PMarg[,,1] <- as.matrix(Piv)
for(i in 1:ns) for(t in 2:TT) PMarg[i,,t]= t(PI[,,i,t])%*%PMarg[i,,t-1]
Pmarg <-apply(PMarg,c(2,3),mean)
}else Pmarg<- NULL
out = c(out,list(V = V1, Ul = Ul, S = S, Pmarg=Pmarg))
}
cat("------------|-------------|-------------|-------------|-------------|-------------|\n");
class(out)="LMlatent"
return(out)
}
Try the LMest package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.