R/startvalues.R

In lava: Latent Variable Models

mgstart <- function(x,p) {
  if (is.list(p)) {
    npar <- with(x, npar+npar.mean+length(unlist(expar)))
    pos <- modelPar(x,seq(npar))$p
    start <- matrix(NA,ncol=npar,nrow=x$ngroup)
    startl <- lapply(pos, function(x) rep(NA,length(x)))
    for (i in seq_len(x$ngroup)) {
      p0 <- p[[i]]
      pos0 <- pos[[i]]
      if (!is.null(names(p0))) {
        ii <- parpos(Model(x)[[i]],names(p0))
        startl[[i]][ii] <- p0[names(ii)]
      } else {
        ii <- seq(min(length(p0),length(startl[[i]])))
        startl[[i]][ii] <- p0[ii]
      }
      start[i,pos0] <- startl[[i]]
    }
    start0 <- apply(start,2,function(x) mean(x,na.rm=TRUE))
  }
  return(start0)
}

##' @export
starter.multigroup <- function(x, starterfun=startvalues2, meanstructure=TRUE,messages=0,...) {
  ## Initial values:
  W <- c() ## Weight-vector
  s <- list()
  for (i in seq_len(x$ngroup)) {
         mydata <- x$data[[i]][,manifest(x$lvm[[i]]),drop=FALSE]
    W <- c(W, nrow(mydata))
    if (nrow(mydata)<3) {
         ii <- index(x$lvm[[i]])
      nn <- ifelse(meanstructure, ii$npar+ii$npar.mean, ii$npar)
      s0 <- rep(1,nn)
    }
    else {
      S <- x$samplestat[[i]]$S
      mu <- if (meanstructure) x$samplestat[[i]]$mu else NULL;
      ##      S <- cov(mydata); mu <- if (meanstructure) colMeans(mydata) else NULL;
      s0 <- do.call(starterfun, list(x$lvm[[i]], S=S, mu=mu,messages=messages))
    }
    s <- c(s, list(s0))
  }

  Wtotal <- sum(W); W <- W/Wtotal

  pg <- vector("list", x$npar); for (i in seq_len(length(pg))) pg[[i]] <- rep(0,x$ngroup)
  meang <- vector("list", x$npar.mean); for (i in seq_len(length(meang))) meang[[i]] <- rep(0,x$ngroup)

  for (i in seq_len(x$ngroup)) {
         pp <- modelPar(x$lvm[[i]],s[[i]])
    pos <- sapply(x$parlist[[i]], function(y) char2num(substr(y,2,nchar(y))))
    for (j in seq_len(length(pos)))
         pg[[ pos[j] ]][i] <-  pp$p[j]

    pos <- sapply(x$meanlist[[i]], function(y) char2num(substr(y,2,nchar(y))))
    ptype <- sapply(x$meanlist[[i]], function(y) substr(y,1,1)=="m")
    if (!(any(ptype)))
         pos <- NULL
    else
         pos <- pos[ptype]
    if (length(pos)>0)
         for (j in seq_len(length(pos))) {
         meang[[ pos[j] ]][i] <-  pp$meanpar[j]
    }
  }

  ## Weighted average
  wp <- unlist(lapply(pg, function(y) {
    ppos <- !is.na(y)
    myweight <- W[ppos]/sum(W[ppos])
    sum(y[ppos]*myweight)
  }))
  wmean <- unlist(lapply(meang, function(y) {
    ppos <- !is.na(y)
    myweight <- W[ppos]/sum(W[ppos])
    sum(y[ppos]*myweight)
  }))
  res <- c(wmean,wp)
  res[!is.finite(res) | is.nan(res) | is.na(res) | is.complex(res)] <- .5
  return(as.numeric(res))
}

startmean <- function(x,p,mu) {
  if (is.null(mu))
         return(p)
  meanpar <- numeric(index(x)$npar.mean)
  mymeans <- vars(x)[index(x)$v1==1]
  midx <- na.omit(match(names(mu),mymeans))
  meanpar[midx] <- mu[midx]
  AP <- matrices(x,p,meanpar)
  nu <- numeric(length(vars(x)))
  nu[vars(x)%in%manifest(x)] <- mu
  meanstart <- ((diag(nrow=nrow(AP$A))-t(AP$A))%*%nu)[index(x)$v1==1]
  names(meanstart) <- vars(x)[index(x)$v1==1]
  return( c(meanstart, p) )
}

## Estimate sub-models (measurement models)
##' @export
`startvalues2` <- function(x, S, mu=NULL, debug=FALSE, messages=0,...) {
  if (messages>0) cat("Obtaining start values...\n")
  S <- reorderdata.lvm(x,S)
  ss <- startvalues(x,S)
  Debug(list("ss=",ss),debug);
  g <- measurement(x,messages=0)
  keep <- c()
  if (length(g)>1) {
    for (i in seq_len(length(g))) {
      if (length(endogenous(g[[i]]))>2)
        keep <- c(keep,i)
    }
    g <- g[keep]
  }
  if (length(g)<2)
    return(startmean(x,ss,mu=mu))
  ## if (messages>0) cat("Fitting marginal measurement models...\n")
  op <- options(warn=-1)
  e <- lapply(g, function(y) {
    estimate(y, data=list(S=S[manifest(y),manifest(y),drop=FALSE], mu=mu[manifest(y)], n=100), control=list(meanstructure=FALSE, starterfun="startvalues", estimator="Simple", method="nlminb1"), optcontrol=list(), debug=FALSE, messages=0)
  })
  for (l in e) {
    ##    a <- coef(l$estimate)[,1]
    a <- coef(l)
    for (i in seq_len(length(a))) {
      pos <- match(names(a)[i],names(ss))
      if (!is.na(pos))
        ss[pos] <- a[i]
    }
  }
  options(op)
  startmean(x,ss,mu=mu)
}

##' @export
startvalues1 <- function(x,S,mu=NULL,tol=1e-6,delta=1e-6,...) {
  p0 <- startvalues(x,S,mu,...)
  p0[index(x)$npar.mean+variances(x)] <- 0.1
  p0[index(x)$npar.mean+offdiags(x)] <- 0
  p0
}

##' @export
startvalues0 <- function(x,S,mu=NULL,tol=1e-6,delta=1e-6,...) {
  p0 <- startvalues(x,S,mu,...)
  A <- t(index(x)$M) ## Adjacency matrix
  P0 <- A0 <- matrix(0,nrow(A),ncol(A))
  A0[,index(x)$eta.idx] <- A[,index(x)$eta.idx]
  diag(P0)[index(x)$endo.idx] <- diag(S)[index(x)$endo.obsidx]/3
  lu <- 0.9
  diag(P0)[index(x)$eta.idx] <- lu##mean(diag(S)[index(x)$endo.idx])/2
  pp <- pars(x,A=t(A0),P=P0,v=rep(1,length(index(x)$vars)))
  pp[pp==1] <- p0[pp==1]
  names(pp) <- coef(x, messages=0, fixed=FALSE, mean=TRUE)[seq_len(length(pp))]
  pp[!is.finite(pp) | is.nan(pp) | is.na(pp)] <- 0.01
  return(pp)
}

##' @export
startvalues <- ## McDonald & Hartmann, 1992
function(x, S, mu=NULL, debug=FALSE, messages=lava.options()$messages, tol=1e-6, delta=1e-6,...) {
  ## As proposed by McDonald & Hartmann, 1992.
  ## Implementation based on John Fox's implementation in the 'sem' R-package
  S <- reorderdata.lvm(x,S)
  if (nrow(S)!=length(manifest(x))) stop("Number of observed variables in data and models does not agree.")
  J <- index(x)$J ## Manifest selection
  P0 <- index(x)$P0 ## covariance 'adjacency'
  A <- t(index(x)$M) ## Adjacency matrix
  n <- nrow(S) ## Number of manifest variables
  m <- nrow(A) ## Number of variables
  A0 <- t(index(x)$M0) ## Adjacency matrix (without fixed parameters)
  obs.idx <- as.vector(J%*%(seq_len(m)));  latent.idx <- setdiff(seq_len(m), obs.idx)
  s <- sqrt(diag(S))
  suppressWarnings(R <- (cov2cor(S))) ## S/outer(s,s)
  C <- P0
  Debug(list("obs.idx", obs.idx), debug)
  C[obs.idx,obs.idx] <- R
  ## Estimates of covariance between latent and manifest variables
  Debug((C), debug)
  for (i in latent.idx) {
    inRelation <- A[obs.idx,i]==1
    for (j in seq_len(length(obs.idx))) {
      Debug((j), debug)
      C[obs.idx[j],i] <- C[i,obs.idx[j]] <- if (any(inRelation)) {
        numerator <- sum(R[j, which(inRelation)])
        denominator <- sqrt(sum(R[which(inRelation), which(inRelation)]))
        numerator/denominator ## as proposed by McDonald & Hartmann
      } else {
        runif(1, .3, .5) ## No arrows => small random covariance
      }
    }
  }
  ## Estimates of covariance between latent variables
  for (i in latent.idx) {
    for (j in latent.idx) {
      C[i,j] <- C[j,i] <-
        if (i==j) {
          1
        } else {
          inRelation.i <- A[obs.idx, i]==1
          inRelation.j <- A[obs.idx, j]==1
          if ((any(inRelation.i)) | (any(inRelation.j))) {
            numerator <- sum(R[which(inRelation.i), which(inRelation.j)])
            denominator <- sqrt( sum(R[which(inRelation.i), which(inRelation.i)])
                                * sum(R[which(inRelation.j), which(inRelation.j)]))
            numerator/(denominator+0.01) ## Avoid division by zero
          } else {
            runif(1, .3, .5)
          }
        }
    }
  }
  if (debug) {
    print("C="); print(C);
  }
  Ahat <- matrix(0,m,m)
  C[is.nan(C)] <- 0
  for (j in seq_len(m)) { ## OLS-estimates
    relation <- A[j,]==1
    if (!any(relation)) next
    Ahat[j, relation] <- tryCatch(Inverse(C[relation,relation] + diag(nrow=sum(relation))*delta,tol=1e-3) %*% C[relation,j], error=function(...) 0)
  }
  Ahat[obs.idx,] <- Ahat[obs.idx,]*matrix(s, n, m)
  Ahat[,obs.idx] <- Ahat[,obs.idx]/matrix(s, m, n, byrow=TRUE)
  Chat <- C
  Chat[obs.idx,] <- Chat[obs.idx,]*matrix(s,n,m)  ##
  Chat[,obs.idx] <- Chat[,obs.idx]*matrix(s,m,n,byrow=TRUE)  ##
  Phat <- (diag(m)-Ahat)%*%Chat%*%t(diag(m)-Ahat)
  ##diag(Phat) <- abs(diag(Phat))
  ## Guarantee PD-matrix:
  Phat[is.nan(Phat) | is.na(Phat)] <- 0
  diag(Phat)[diag(Phat)==0] <- 1
  eig <- eigen(Phat)
  L <- abs(eig$values); L[L<1e-3] <- 1e-3
  Phat <- eig$vectors%*%diag(L,ncol=ncol(eig$vectors))%*%t(eig$vectors)

  Debug(list("start=",start), debug)
  start <- pars(x, A=t(Ahat*A0), P=(Phat*P0))
  names(start) <- coef(x, messages=0, fixed=FALSE, mean=FALSE)[seq_len(length(start))]
  res <- startmean(x,start,mu)
  res[!is.finite(res) | is.nan(res) | is.na(res)] <- 1
  res
}