Nothing
###{{{ 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) {
parcount <- 0
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 <- paste("par",i,sep="")
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+1: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(paste("Missing data encountered. Going for complete-case analysis", sep=""))
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)]
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)])
opt <- options(warn=-1); pidx <- is.na(as.numeric(parname)); options(opt)
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] <- paste("p",counter,sep="")
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] <- as.numeric(substr(pres[[i]][j],2,nchar(pres[[i]][j])))
} else {
val <- paste("p",counter,sep="")
pres[[i]][j] <- val
pres0[[i]][j] <- counter
usedname[pidx,2] <- val
counter <- counter+1
}
} else {
pres[[i]][j] <- NA
}
}
}
mypar <- paste("p",seq_len(nfree),sep="")
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)])
opt <- options(warn=-1); midx <- is.na(as.numeric(meanparname)); options(opt)
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)))
nmean <- lapply(means,length)
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] <- paste("m",counter,sep="")
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] <- as.numeric(substr(res[[i]][j],2,nchar(res[[i]][j]))) +
nfree.mean
##nmean[[i]]
} else {
if (!is.na(usedname[midx,2])) {
res[[i]][j] <- usedname[midx,2]
res0[[i]][j] <- as.numeric(substr(res[[i]][j],2,nchar(res[[i]][j])))
} else {
val <- paste("m",counter,sep="")
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
N <- nfree+nfree.mean
m0 <- p0 <- c()
coefs <- coefsm <- mm0 <- mm <- pp0 <- pp <- c()
for (i in seq_len(length(myparPos))) {
mi <- mymeanPos[[i]]
nmi <- length(mi)
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(i,coef(lvms[[i]],fix=FALSE,mean=FALSE)[pp0[[i]]],sep="@"))
if (length(m1)>0) {
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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.