###{{{ multigroup
##' @export
multigroup <- function(models, datasets, fix, exo.fix=TRUE, keep=NULL, missing=FALSE, ...) {
nm <- length(models)
if (nm!=length(datasets)) stop("Supply dataset for each model")
if (nm<2) stop("Two or more groups neeeded")
mynames <- names(models)
## Check for random slopes
xfix <- list()
for (i in seq_len(nm)) {
x0 <- models[[i]]
data0 <- datasets[[i]]
xfix0 <- colnames(data0)[(colnames(data0)%in%parlabels(x0,exo=TRUE))]
xfix <- c(xfix, list(xfix0))
}
if (missing(fix)) {
fix <- !any(unlist(lapply(xfix, function(x) length(x)>0)))
}
for (i in seq_len(nm)) {
x0 <- models[[i]]
data0 <- datasets[[i]]
if (length(exogenous(x0)>0)) {
catx <- categorical2dummy(x0,data0)
models[[i]] <- catx$x; datasets[[i]] <- catx$data
}
if (!lava.options()$exogenous) exogenous(models[[i]]) <- NULL
}
models.orig <- NULL
######################
### MLE with MAR mechanism
######################
if (missing) {
reservedpars <- c()
mynpar <- c()
for (i in seq_len(nm)) {
## Fix some parameters (predictors,latent variables,...)
d0 <- datasets[[i]][1,,drop=FALSE]; d0[,] <- 1
if (fix)
models[[i]] <- fixsome(models[[i]], exo.fix=exo.fix, measurement.fix=fix, data=d0)
## Find named/labelled parameters
rpar <- unique(parlabels(models[[i]]))
reservedpars <- c(reservedpars, rpar)
mynpar <- c(mynpar, with(index(models[[1]]), npar+npar.mean+npar.ex))
}; reservedpars <- unique(reservedpars)
nonamepar <- sum(mynpar)
## Find unique parameter-names for all parameters
newpars <- c()
i <- 0
pos <- 1
while(pos<=nonamepar) {
i <- i+1
newname <- paste0("par",i)
if (!(newname%in%reservedpars)) {
newpars <- c(newpars,newname)
pos <- pos+1
}
}
pos <- 0
models0 <- list()
datasets0 <- list()
complidx <- c()
nmodels <- 0
modelclass <- c()
nmis <- c()
for (i in seq_len(nm)) {
myvars <- unlist(intersect(colnames(datasets[[i]]),c(vars(models[[i]]),xfix[[i]],keep)))
mydata <- datasets[[i]][,myvars]
if (any(is.na(mydata))) {
if (i>1) pos <- pos+mynpar[i-1]
models[[i]] <- baptize(models[[i]],newpars[pos+seq_len(mynpar[i])] ,overwrite=FALSE)
val <- missingModel(models[[i]],mydata,fix=FALSE,keep=keep,...)
nmodels <- c(nmodels,length(val$models))
complidx <- c(complidx,val$pattern.allcomp+nmodels[i]+1)
nmis0 <- rowSums(val$patterns);
allmis <- which(nmis0==ncol(val$patterns))
if (length(allmis)>0) nmis0 <- nmis0[-allmis]
nmis <- c(nmis,nmis0)
datasets0 <- c(datasets0, val$datasets)
models0 <- c(models0, val$models)
modelclass <- c(modelclass,rep(i,length(val$models)))
} else {
datasets0 <- c(datasets0, list(mydata))
models0 <- c(models0, list(models[[i]]))
modelclass <- c(modelclass,i)
nmis <- c(nmis,0)
}
}
models.orig <- models
suppressWarnings(
val <- multigroup(models0,datasets0,fix=FALSE,missing=FALSE,exo.fix=TRUE,...)
)
val$models.orig <- models.orig; val$missing <- TRUE
val$complete <- complidx-1
val$mnames <- mynames
attributes(val)$modelclass <- modelclass
attributes(val)$nmis <- nmis
return(val)
}
######################
### Usual analysis:
######################
warned <- FALSE
for (i in seq_len(nm)) {
if (inherits(datasets[[i]],c("data.frame","matrix"))) {
myvars <- intersect(colnames(datasets[[i]]),c(vars(models[[i]]),xfix[[i]],keep))
if (any(is.na(datasets[[i]][,myvars]))) {
if (!warned)
warning(paste0("Missing data encountered. Going for complete-case analysis"))
warned <- TRUE
datasets[[i]] <- na.omit(datasets[[i]][,myvars,drop=FALSE])
}
}
}
exo <- exogenous(models)
means <- lvms <- As <- Ps <- ps <- exs <- datas <- samplestat <- list()
for (i in seq_len(nm)) {
if (!is.null(exogenous(models[[i]]))) {
if (any(is.na(exogenous(models[[i]])))) {
exogenous(models[[i]]) <- exo
}
}
mydata <- datasets[[i]]
mymodel <- fixsome(models[[i]], data=mydata, measurement.fix=fix, exo.fix=exo.fix)
mymodel <- updatelvm(mymodel,zeroones=TRUE,deriv=TRUE)
P <- index(mymodel)$P1; P[P==0] <- NA
P[!is.na(P) & !is.na(mymodel$covpar)] <- mymodel$covpar[!is.na(P) & !is.na(mymodel$covpar)]
A <- index(mymodel)$M1; A[A==0] <- NA
A[!is.na(A) & !is.na(mymodel$par)] <- mymodel$par[!is.na(A) & !is.na(mymodel$par)]
mu <- unlist(mymodel$mean)[which(index(mymodel)$v1==1)]
#ex <- names(mymodel$expar)[which(index(mymodel)$e1==1)]
ex <- mymodel$exfix
if (length(ex)>0) {
if (any(is.na(ex))) ex[is.na(ex)] <- mymodel$expar[is.na(ex)]
ex <- ex[which(index(mymodel)$e1==1)]
}
p <- pars(mymodel, A, P, e=ex)
p[p=="1"] <- NA
means <- c(means, list(mu))
lvms <- c(lvms, list(mymodel))
datas <- c(datas, list(mydata))
samplestat <- c(samplestat, list(procdata.lvm(models[[i]],data=mydata)))
As <- c(As, list(A))
Ps <- c(Ps, list(P))
ps <- c(ps, list(p))
exs <- c(exs, list(ex))
};
######
pp <- unlist(ps)
parname <- unique(pp[!is.na(pp)])
pidx <- is.na(char2num(parname))
parname <- unique(unlist(pp[!is.na(pp)]));
nfree <- sum(is.na(pp)) + length(parname)
if (nfree>0) {
pp0 <- lapply(ps, is.na)
usedname <- cbind(parname, rep(NA,length(parname)))
counter <- 1
pres <- pres0 <- pp0
for (i in seq_len(length(pp0))) {
if (length(pp0[[i]]>0))
for (j in seq_len(length(pp0[[i]]))) {
pidx <- match(ps[[i]][j],parname)
if (pp0[[i]][j]) {
pres[[i]][j] <- paste0("p",counter)
pres0[[i]][j] <- counter
counter <- counter+1
} else if (!is.na(pidx)) {
if (!is.na(usedname[pidx,2])) {
pres[[i]][j] <- usedname[pidx,2]
pres0[[i]][j] <- char2num(substr(pres[[i]][j],2,nchar(pres[[i]][j])))
} else {
val <- paste0("p",counter)
pres[[i]][j] <- val
pres0[[i]][j] <- counter
usedname[pidx,2] <- val
counter <- counter+1
}
} else {
pres[[i]][j] <- NA
}
}
}
mypar <- paste0("p",seq_len(nfree))
myparPos <- pres0
myparpos <- pres
myparlist <- lapply(pres, function(x) x[!is.na(x)])
} else {
myparPos <- NULL
mypar <- NULL
myparpos <- NULL
myparlist <- NULL
}
### Mean parameter
mm <- unlist(means)
meanparname <- unique(mm[!is.na(mm)])
midx <- is.na(char2num(meanparname));
meanparname <- meanparname[midx]
any.mean <- sum(is.na(mm)) + length(meanparname)
nfree.mean <- sum(is.na(mm)) + length(setdiff(meanparname,parname))
## mean.fixed <- na.omit(match(parname,mm))
mean.omit <- lapply(means,function(x) na.omit(match(parname,x)))
if (any.mean>0) {
mm0 <- lapply(means, is.na)
usedname <- cbind(meanparname, rep(NA,length(meanparname)))
counter <- 1
res0 <- res <- mm0
for (i in seq_len(length(mm0))) {
if (length(mm0[[i]])>0)
for (j in seq_len(length(mm0[[i]]))) {
midx <- match(means[[i]][j],meanparname)
if (mm0[[i]][j]) {
res[[i]][j] <- paste0("m",counter)
res0[[i]][j] <- counter
counter <- counter+1
} else if (!is.na(midx)) {
pidx <- match(meanparname[midx],pp)
if (!is.na(pidx)) {
res[[i]][j] <- unlist(myparlist)[pidx]
res0[[i]][j] <- char2num(substr(res[[i]][j],2,nchar(res[[i]][j]))) +
nfree.mean
} else {
if (!is.na(usedname[midx,2])) {
res[[i]][j] <- usedname[midx,2]
res0[[i]][j] <- char2num(substr(res[[i]][j],2,nchar(res[[i]][j])))
} else {
val <- paste0("m",counter)
res[[i]][j] <- val
res0[[i]][j] <- counter
usedname[midx,2] <- val
counter <- counter+1
}
}
} else {
res[[i]][j] <- NA
}
}
}
mymeanPos <- res0
mymeanpos <- res
mymeanlist <- lapply(res, function(x) x[!is.na(x)])
mymean <- unique(unlist(mymeanlist))
} else {
mymeanPos <- NULL
mymean <- NULL
mymeanpos <- NULL
mymeanlist <- NULL
}
### Extra parameters
m0 <- p0 <- c()
coefs <- coefsm <- mm0 <- mm <- pp0 <- pp <- c()
for (i in seq_len(length(myparPos))) {
mi <- mymeanPos[[i]]
pi <- myparPos[[i]]
p1 <- setdiff(pi,p0)
p0 <- c(p0,p1)
## pp0 <- c(pp0,list(match(p1,pi)+nfree.mean))
pp0 <- c(pp0,list(match(p1,pi)))
if (length(mean.omit[[i]])>0) mi <- mi[-mean.omit[[i]]]
m1 <- setdiff(mi,m0)
m0 <- c(m0,m1)
mm0 <- c(mm0,list(match(m1,mi)))
pp <- c(pp,list(c(m1,p1+nfree.mean)))
if (length(p1)>0)
coefs <- c(coefs,paste(coef(lvms[[i]],fix=FALSE,mean=FALSE)[pp0[[i]]],i,sep="@"))
#coefs <- c(coefs,paste(i,coef(lvms[[i]],fix=FALSE,mean=FALSE)[pp0[[i]]],sep="@"))
if (length(m1)>0) {
coefsm0 <- paste(coef(lvms[[i]],fix=FALSE,mean=TRUE)[mm0[[i]]],i,sep="@")
#coefsm0 <- paste(i,coef(lvms[[i]],fix=FALSE,mean=TRUE)[mm0[[i]]],sep="@")
coefsm <- c(coefsm,coefsm0)
}
}
coefs <- c(coefsm,coefs)
res <- list(npar=nfree, npar.mean=nfree.mean,
ngroup=length(lvms), names=mynames,
lvm=lvms, data=datas, samplestat=samplestat,
A=As, P=Ps, expar=exs,
meanpar=names(mu), name=coefs, coef=pp, coef.idx=pp0,
par=mypar, parlist=myparlist, parpos=myparpos,
mean=mymean, meanlist=mymeanlist, meanpos=mymeanpos,
parposN=myparPos,
meanposN=mymeanPos,
models.orig=models.orig, missing=missing
)
class(res) <- "multigroup"
checkmultigroup(res)
return(res)
}
###}}}
###{{{ checkmultigroup
checkmultigroup <- function(x) {
## Check validity:
for (i in seq_len(x$ngroup)) {
if (nrow(x$data[[i]])<2) {
warning("With only one observation in the group, all parameters should be inherited from another a group!")
}
}
}
###}}} checkmultigroup
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.