1###{{{ multigroup
2
3##' @export
4multigroup <- function(models, datasets, fix, exo.fix=TRUE, keep=NULL, missing=FALSE, ...) {
5    nm <- length(models)
6    if (nm!=length(datasets)) stop("Supply dataset for each model")
7    if (nm<2) stop("Two or more groups neeeded")
8    mynames <- names(models)
9
10    ## Check for random slopes
11    xfix <- list()
12    for (i in seq_len(nm)) {
13        x0 <- models[[i]]
14        data0 <- datasets[[i]]
15        xfix0 <- colnames(data0)[(colnames(data0)%in%parlabels(x0,exo=TRUE))]
16        xfix <- c(xfix, list(xfix0))
17    }
18    if (missing(fix)) {
19        fix <- !any(unlist(lapply(xfix, function(x) length(x)>0)))
20    }
21
22    for (i in seq_len(nm)) {
23        x0 <- models[[i]]
24        data0 <- datasets[[i]]
25        if (length(exogenous(x0)>0)) {
26            catx <- categorical2dummy(x0,data0)
27            models[[i]] <- catx$x; datasets[[i]] <- catx$data
28        }
29        if (!lava.options()$exogenous) exogenous(models[[i]]) <- NULL
30    }
31
32    models.orig <- NULL
33    ######################
34    ### MLE with MAR mechanism
35    ######################
36    if (missing) {
37
38        reservedpars <- c()
39        mynpar <- c()
40        for (i in seq_len(nm)) {
41            ## Fix some parameters (predictors,latent variables,...)
42
43            d0 <- datasets[[i]][1,,drop=FALSE]; d0[,] <- 1
44            if (fix)
45                models[[i]] <- fixsome(models[[i]], exo.fix=exo.fix, measurement.fix=fix, data=d0)
46            ## Find named/labelled parameters
47            rpar <- unique(parlabels(models[[i]]))
48            reservedpars <- c(reservedpars, rpar)
49            mynpar <- c(mynpar, with(index(models[[1]]), npar+npar.mean+npar.ex))
50        }; reservedpars <- unique(reservedpars)
51        nonamepar <- sum(mynpar)
52        ## Find unique parameter-names for all parameters
53        newpars <- c()
54        i <- 0
55        pos <- 1
56        while(pos<=nonamepar) {
57            i <- i+1
58            newname <- paste0("par",i)
59            if (!(newname%in%reservedpars)) {
60                newpars <- c(newpars,newname)
61                pos <- pos+1
62            }
63        }
64
65        pos <- 0
66        models0 <- list()
67        datasets0 <- list()
68        complidx <- c()
69        nmodels <- 0
70        modelclass <- c()
71        nmis <- c()
72        for (i in seq_len(nm)) {
73            myvars <- unlist(intersect(colnames(datasets[[i]]),c(vars(models[[i]]),xfix[[i]],keep)))
74            mydata <- datasets[[i]][,myvars]
75            if (any(is.na(mydata))) {
76                if (i>1) pos <- pos+mynpar[i-1]
77                models[[i]] <- baptize(models[[i]],newpars[pos+seq_len(mynpar[i])] ,overwrite=FALSE)
78                val <- missingModel(models[[i]],mydata,fix=FALSE,keep=keep,...)
79                nmodels <- c(nmodels,length(val$models))
80                complidx <- c(complidx,val$pattern.allcomp+nmodels[i]+1)
81                nmis0 <- rowSums(val$patterns);
82                allmis <- which(nmis0==ncol(val$patterns))
83                if (length(allmis)>0) nmis0 <- nmis0[-allmis]
84                nmis <- c(nmis,nmis0)
85                datasets0 <- c(datasets0, val$datasets)
86                models0 <- c(models0, val$models)
87                modelclass <- c(modelclass,rep(i,length(val$models)))
88            } else {
89                datasets0 <- c(datasets0, list(mydata))
90                models0 <- c(models0, list(models[[i]]))
91                modelclass <- c(modelclass,i)
92                nmis <- c(nmis,0)
93            }
94        }
95
96        models.orig <- models
97
98        suppressWarnings(
99            val <- multigroup(models0,datasets0,fix=FALSE,missing=FALSE,exo.fix=TRUE,...)
100        )
101        val$models.orig <- models.orig; val$missing <- TRUE
102        val$complete <- complidx-1
103        val$mnames <- mynames
104        attributes(val)$modelclass <- modelclass
105        attributes(val)$nmis <- nmis
106        return(val)
107    }
108
109
110    ######################
111    ### Usual analysis:
112    ######################
113    warned <- FALSE
114    for (i in seq_len(nm)) {
115        if (inherits(datasets[[i]],c("data.frame","matrix"))) {
116            myvars <- intersect(colnames(datasets[[i]]),c(vars(models[[i]]),xfix[[i]],keep))
117            if (any(is.na(datasets[[i]][,myvars]))) {
118                if (!warned)
119                    warning(paste0("Missing data encountered. Going for complete-case analysis"))
120                warned  <- TRUE
121                datasets[[i]] <- na.omit(datasets[[i]][,myvars,drop=FALSE])
122            }
123        }
124    }
125
126    exo <- exogenous(models)
127    means <- lvms <- As <- Ps <- ps <- exs <- datas <- samplestat <- list()
128    for (i in seq_len(nm)) {
129
130        if (!is.null(exogenous(models[[i]]))) {
131            if (any(is.na(exogenous(models[[i]])))) {
132                exogenous(models[[i]]) <- exo
133            }
134        }
135
136        mydata <- datasets[[i]]
137        mymodel <- fixsome(models[[i]], data=mydata, measurement.fix=fix, exo.fix=exo.fix)
138        mymodel <- updatelvm(mymodel,zeroones=TRUE,deriv=TRUE)
139
140        P <- index(mymodel)$P1; P[P==0] <- NA
141        P[!is.na(P) & !is.na(mymodel$covpar)] <- mymodel$covpar[!is.na(P) & !is.na(mymodel$covpar)]
142
143        A <- index(mymodel)$M1; A[A==0] <- NA
144        A[!is.na(A) & !is.na(mymodel$par)] <- mymodel$par[!is.na(A) & !is.na(mymodel$par)]
145
146        mu <- unlist(mymodel$mean)[which(index(mymodel)$v1==1)]
147        #ex <- names(mymodel$expar)[which(index(mymodel)$e1==1)]
148        ex <- mymodel$exfix
149        if (length(ex)>0) {
150            if (any(is.na(ex))) ex[is.na(ex)] <- mymodel$expar[is.na(ex)]
151            ex <- ex[which(index(mymodel)$e1==1)]
152        }
153
154        p <- pars(mymodel, A, P, e=ex)
155        p[p=="1"] <- NA
156
157        means <- c(means, list(mu))
158        lvms <- c(lvms, list(mymodel))
159        datas <- c(datas, list(mydata))
160        samplestat <- c(samplestat, list(procdata.lvm(models[[i]],data=mydata)))
161        As <- c(As, list(A))
162        Ps <- c(Ps, list(P))
163        ps <- c(ps, list(p))
164        exs <- c(exs, list(ex))
165    };
166
167    ######
168    pp <- unlist(ps)
169    parname <- unique(pp[!is.na(pp)])
170    pidx <- is.na(char2num(parname))
171    parname <- unique(unlist(pp[!is.na(pp)]));
172    nfree <- sum(is.na(pp)) + length(parname)
173
174    if (nfree>0) {
175        pp0 <- lapply(ps, is.na)
176        usedname <- cbind(parname, rep(NA,length(parname)))
177        counter <- 1
178        pres <- pres0 <- pp0
179        for (i in seq_len(length(pp0))) {
180            if (length(pp0[[i]]>0))
181                for (j in seq_len(length(pp0[[i]]))) {
182                    pidx <- match(ps[[i]][j],parname)
183                    if (pp0[[i]][j]) {
184                        pres[[i]][j] <- paste0("p",counter)
185                        pres0[[i]][j] <- counter
186                        counter <- counter+1
187                    } else if (!is.na(pidx)) {
188                        if (!is.na(usedname[pidx,2])) {
189                            pres[[i]][j] <- usedname[pidx,2]
190                            pres0[[i]][j] <- char2num(substr(pres[[i]][j],2,nchar(pres[[i]][j])))
191                        } else {
192                            val <- paste0("p",counter)
193                            pres[[i]][j] <- val
194                            pres0[[i]][j] <- counter
195                            usedname[pidx,2] <- val
196                            counter <- counter+1
197                        }
198                    } else {
199                        pres[[i]][j] <- NA
200                    }
201                }
202        }
203        mypar <- paste0("p",seq_len(nfree))
204        myparPos <- pres0
205        myparpos <- pres
206        myparlist <- lapply(pres, function(x) x[!is.na(x)])
207    } else {
208        myparPos <- NULL
209        mypar <- NULL
210        myparpos <- NULL
211        myparlist <- NULL
212    }
213
214    ### Mean parameter
215
216    mm <- unlist(means)
217    meanparname <- unique(mm[!is.na(mm)])
218    midx <- is.na(char2num(meanparname));
219    meanparname <- meanparname[midx]
220    any.mean <- sum(is.na(mm)) + length(meanparname)
221    nfree.mean <- sum(is.na(mm)) + length(setdiff(meanparname,parname))
222    ## mean.fixed <- na.omit(match(parname,mm))
223    mean.omit <- lapply(means,function(x) na.omit(match(parname,x)))
224
225    if (any.mean>0) {
226        mm0 <- lapply(means, is.na)
227        usedname <- cbind(meanparname, rep(NA,length(meanparname)))
228        counter <- 1
229        res0 <- res <- mm0
230        for (i in seq_len(length(mm0))) {
231            if (length(mm0[[i]])>0)
232                for (j in seq_len(length(mm0[[i]]))) {
233                    midx <- match(means[[i]][j],meanparname)
234                    if (mm0[[i]][j]) {
235                        res[[i]][j] <- paste0("m",counter)
236                        res0[[i]][j] <- counter
237                        counter <- counter+1
238                    } else if (!is.na(midx)) {
239                        pidx <- match(meanparname[midx],pp)
240                        if (!is.na(pidx)) {
241                            res[[i]][j] <- unlist(myparlist)[pidx]
242                            res0[[i]][j] <- char2num(substr(res[[i]][j],2,nchar(res[[i]][j]))) +
243                                nfree.mean
244                        } else {
245                            if (!is.na(usedname[midx,2])) {
246                                res[[i]][j] <- usedname[midx,2]
247                                res0[[i]][j] <- char2num(substr(res[[i]][j],2,nchar(res[[i]][j])))
248                            } else {
249                                val <- paste0("m",counter)
250                                res[[i]][j] <- val
251                                res0[[i]][j] <- counter
252                                usedname[midx,2] <- val
253                                counter <- counter+1
254                            }
255                        }
256                    } else {
257                        res[[i]][j] <- NA
258                    }
259                }
260        }
261        mymeanPos <- res0
262        mymeanpos <- res
263        mymeanlist <- lapply(res, function(x) x[!is.na(x)])
264        mymean <- unique(unlist(mymeanlist))
265    } else {
266        mymeanPos <- NULL
267        mymean <- NULL
268        mymeanpos <- NULL
269        mymeanlist <- NULL
270    }
271
272    ### Extra parameters
273
274    m0 <- p0 <- c()
275    coefs <- coefsm <- mm0 <- mm <- pp0 <- pp <- c()
276    for (i in seq_len(length(myparPos))) {
277        mi <- mymeanPos[[i]]
278        pi <- myparPos[[i]]
279        p1 <- setdiff(pi,p0)
280        p0 <- c(p0,p1)
281        ##    pp0 <- c(pp0,list(match(p1,pi)+nfree.mean))
282        pp0 <- c(pp0,list(match(p1,pi)))
283        if (length(mean.omit[[i]])>0) mi <- mi[-mean.omit[[i]]]
284        m1 <- setdiff(mi,m0)
285        m0 <- c(m0,m1)
286        mm0 <- c(mm0,list(match(m1,mi)))
287        pp <- c(pp,list(c(m1,p1+nfree.mean)))
288        if (length(p1)>0)
289            coefs <- c(coefs,paste(coef(lvms[[i]],fix=FALSE,mean=FALSE)[pp0[[i]]],i,sep="@"))
290            #coefs <- c(coefs,paste(i,coef(lvms[[i]],fix=FALSE,mean=FALSE)[pp0[[i]]],sep="@"))
291        if (length(m1)>0) {
292            coefsm0 <- paste(coef(lvms[[i]],fix=FALSE,mean=TRUE)[mm0[[i]]],i,sep="@")
293            #coefsm0 <- paste(i,coef(lvms[[i]],fix=FALSE,mean=TRUE)[mm0[[i]]],sep="@")
294            coefsm <- c(coefsm,coefsm0)
295        }
296    }
297    coefs <- c(coefsm,coefs)
298
299    res <- list(npar=nfree, npar.mean=nfree.mean,
300                ngroup=length(lvms), names=mynames,
301                lvm=lvms, data=datas, samplestat=samplestat,
302                A=As, P=Ps, expar=exs,
303                meanpar=names(mu), name=coefs, coef=pp, coef.idx=pp0,
304                par=mypar, parlist=myparlist,  parpos=myparpos,
305                mean=mymean, meanlist=mymeanlist, meanpos=mymeanpos,
306                parposN=myparPos,
307                meanposN=mymeanPos,
308                models.orig=models.orig, missing=missing
309                )
310    class(res) <- "multigroup"
311    checkmultigroup(res)
312    return(res)
313}
314
315###}}}
316
317###{{{ checkmultigroup
318
319checkmultigroup <- function(x) {
320    ## Check validity:
321    for (i in seq_len(x$ngroup)) {
322        if (nrow(x$data[[i]])<2) {
323            warning("With only one observation in the group, all parameters should be inherited from another a group!")
324        }
325    }
326}
327
328###}}} checkmultigroup
329