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