Nothing
###{{{ missingModel
missingModel <- function(model,data,var=endogenous(model),fix=FALSE,type=2,keep=NULL,weights=NULL,data2=NULL,cluster=NULL,...) {
if (!inherits(model,"lvm")) stop("Needs a lvm-object")
if (type==3) {
var <- manifest(model)
}
data.mis <- is.na(data[,var,drop=FALSE])
colnames(data.mis) <- var
patterns <- unique(data.mis,MARGIN=1)
mis.type <- apply(data.mis,1,
function(x) which(apply(patterns,1,function(y) identical(x,y))))
pattern.allmis <- which(apply(patterns,1,all)) ## Remove entry with all missing
data2.list.provided <- is.list(data2) & !is.data.frame(data2)
weights.list.provided <- is.list(weights) & !is.data.frame(weights)
models <- datasets <- datasets2 <- weightsdata <- clusters <- c()
mymodel <- baptize(model)
pattern.compl <- 0
count <- 0
A <- index(model)$A
topendo <- endogenous(model,top=TRUE)
exo <- exogenous(model)
exclude <- c()
warned <- FALSE
for (i in setdiff(seq_len(nrow(patterns)),pattern.allmis)) {
exoremove <- c()
count <- count+1
mypattern <- patterns[i,]
m0 <- mymodel;
if (any(mypattern)) {
latent(m0) <- colnames(data.mis)[mypattern]
if (type>1) {
mytop <- intersect(topendo,colnames(data.mis)[mypattern])
if (!is.null(mytop)) {
rmvar(m0) <- mytop
for (xx in exo) {
## If exogenous variable only have effect on missing variables,
## then remove it from the model
if (all(c(rownames(A)[A[xx,]==1])%in%mytop) &&
!(xx%in%m0$par)
##&& !(xx%in%names(index(m0))$parval)
) {
exoremove <- c(exoremove,xx)
rmvar(m0) <- xx
}
}
}
}
} else
pattern.compl <- count
## d0 <- data[mis.type==i,manifest(m0),drop=FALSE];
d0 <- data[which(mis.type==i),c(manifest(m0),keep),drop=FALSE];
w0 <- NULL
if (!is.null(weights))
if (!weights.list.provided) {
if (NCOL(weights)==1) {
w0 <- as.matrix(weights) [which(mis.type==i)];
} else {
w0.var <- intersect(manifest(m0),colnames(weights))
w0 <- weights[which(mis.type==i),w0.var,drop=FALSE];
}
}
if (!data2.list.provided) {
w02.var <- intersect(manifest(m0),colnames(data2))
w02 <- data2[which(mis.type==i),w02.var,drop=FALSE];
}
clust0 <- cluster[which(mis.type==i)]
ex0 <- exogenous(m0) <- setdiff(exo,exoremove)
xmis <- which(apply(d0[,ex0,drop=FALSE],1,function(x) any(is.na(x))))
if (length(xmis)>0) {
misx <- ex0[apply(d0[xmis,ex0,drop=FALSE],2,function(x) any(is.na(x)))]
if (!warned)
warning("Missing exogenous variables: ", paste(misx,collapse=","),
". Removing rows...")
warned <- TRUE
d0 <- d0[-xmis,,drop=FALSE]
if (NCOL(weights)==1) {
w0 <- w0[-xmis]
} else {
w0 <- w0[-xmis,,drop=FALSE]
}
clust0 <- clust0[-xmis]
w02 <- w02[-xmis,,drop=FALSE]
}
if (length(misx <- intersect(ex0,latent(m0)))>0) {
warning("Missing exogenous variables:", paste(misx,collapse=","),
"! Remove manually!.")
}
## else
{
if( sum(unlist(index(m0)[c("npar","npar.mean")]))>0 ) {
models <- c(models, list(m0))
datasets <- c(datasets, list(d0))
if (!weights.list.provided)
weightsdata <- c(weightsdata, list(w0))
if (!data2.list.provided)
datasets2 <- c(datasets2, list(w02))
clusters <- c(clusters, list(clust0))
} else {
exclude <- c(exclude,count)
}
}
}
rmset <- c()
for (i in seq_len(length(datasets))) {
if (nrow(datasets[[i]])==0) rmset <- c(rmset,i)
}
if (length(rmset)>0) {
models[[rmset]] <- NULL
datasets[[rmset]] <- NULL
weightsdata[[rmset]] <- NULL
datasets2[[rmset]] <- NULL
clusters[[rmset]] <- NULL
patterns <- patterns[-rmset,,drop=FALSE]
}
if (data2.list.provided) datasets2 <- data2
if (weights.list.provided) weightsdata <- weights
Patterns <- patterns
if (length(exclude)>0)
Patterns <- Patterns[-exclude,]
pattern.allcomp<- which(apply(Patterns,1,function(x) all(!x))) ## Complete cases
res <- list(models=models, datasets=datasets,
weights=weightsdata,
data2=datasets2,
clusters=clusters,
patterns=Patterns,
pattern.compl=pattern.compl,
pattern.allmis=pattern.allmis,
pattern.allcomp=pattern.allcomp,
mis.type=mis.type)
return(res)
}
###}}}
###{{{ estimate.MAR.lvm
##' @export
estimate.MAR <- function(x,data,which=endogenous(x),fix,type=2,startcc=FALSE,control=list(),messages=lava.options()$messages,weights,data2,cluster,onlymodel=FALSE,estimator="gaussian",hessian=TRUE,keep=NULL,...) {
cl <- match.call()
Debug("estimate.MAR")
redvar <- intersect(intersect(parlabels(x),latent(x)),colnames(data))
if (length(redvar)>0 & (messages>0))
warning(paste("Remove latent variable colnames from dataset",redvar))
xfix <- setdiff(colnames(data)[(colnames(data)%in%parlabels(x,exo=TRUE))],latent(x))
if (missing(fix))
fix <- ifelse(length(xfix)>0,FALSE,TRUE)
S <- diag(nrow=length(manifest(x)));
mu <- rep(0,nrow(S));
K <- length(exogenous(x))
vnames <- index(x)$manifest
names(mu) <- rownames(S) <- colnames(S) <- vnames
if (K>0) {
xx <- subset(Model(x),exogenous(x))
exogenous(xx) <- NULL
covfix(xx, vars(xx)) <- NA
xx <- covariance(xx,exogenous(x),exogenous(x))
datax <- data[,exogenous(x),drop=FALSE]
exo.idx <- match(exogenous(x),manifest(x))
mu0 <- colMeans(datax,na.rm=TRUE)
cov0 <- cov(datax,use="pairwise.complete.obs")*(nrow(datax)-1)/nrow(datax)
cov0upper <- cov0[upper.tri(cov0,diag=TRUE)]
exogenous(xx) <- NULL
coefpos <- matrices(xx,seq_len(K*(K-1)/2+K))$P
ii <- coefpos[upper.tri(coefpos,diag=TRUE)]
start <- c(mu0, cov0upper[order(ii)])
S[exo.idx,exo.idx] <- cov0
mu[exo.idx] <- mu0
## message("\n")
}
x0 <- x
x <- fixsome(x, measurement.fix=fix, exo.fix=TRUE, S=S, mu=mu, n=1)
if (messages>1)
message("Identifying missing patterns...")
val <- missingModel(x,data,var=which,type=type,keep=c(keep,xfix),weights=weights,data2=data2,cluster=cluster,...)
if (messages>1)
message("\n")
if (nrow(val$patterns)==1) {
res <- estimate(x, data=data, fix=fix, weights=weights, data2=data2,
estimator=estimator, messages=messages, control=control, ...)
return(res)
}
if (startcc & is.null(control$start)) {
if (messages>1)
message("Obtaining starting value...")
start0 <- rep(1,sum(unlist(index(x)[c("npar","npar.mean")])))
mystart <- tryCatch(
(estimate(x,data=na.omit(data),messages=0,
weights=weights,data2=data2,estimator=estimator,quick=TRUE,...
)),
error=function(e) rep(1,sum(unlist(index(x)[c("npar","npar.mean")])))
)
control$start <- mystart
if (messages>1)
message("\n")
}
if (is.null(control$meanstructure))
control$meanstructure <- TRUE
mg0 <- with(val, suppressWarnings(multigroup(models,datasets,fix=FALSE,exo.fix=FALSE,missing=FALSE)))
if (!is.null(names(control$start))) {
parorder1 <- attributes(parpos(mg0,p=names(control$start)))$name
paridx <- match(parorder1,names(control$start))
## newpos <- paridx[which(!is.na(paridx))]
start0 <- control$start
start0[which(!is.na(paridx))] <- control$start[na.omit(paridx)]
names(start0)[which(!is.na(paridx))] <- names(control$start[na.omit(paridx)])
control$start <- start0
}
if (onlymodel) return(list(mg=mg0,val=val,weights=val$weights,data2=val$data2,cluster=val$clusters))
if (all(unlist(lapply(val$weights,is.null)))) val$weights <- NULL
if (all(unlist(lapply(val$data2,is.null)))) val$data2 <- NULL
if (all(unlist(lapply(val$clusters,is.null)))) val$clusters <- NULL
e.mis <- estimate(mg0,control=control,messages=messages,
weights=val$weights,data2=val$data2,
cluster=val$clusters,estimator=estimator,...)
cc <- coef(e.mis,type=1)
mynames <- c()
if (e.mis$model$npar.mean>0)
mynames <- c(mynames,paste0("m",seq_len(e.mis$model$npar.mean)))
if (e.mis$model$npar>0)
mynames <- c(mynames,paste0("p",seq_len(e.mis$model$npar)))
rownames(cc) <- mynames
mycc <- val$pattern.allcomp ## Position of complete-case model
nmis <- with(val, as.numeric(table(mis.type)[pattern.allmis])) ## Number of completely missing observations
if (length(nmis)>0 & length(mycc)>0) ## Any individuals with all missing?
if (val$pattern.allmis<mycc)
mycc <- mycc-1
if (length(xfix)>0) {
nrow <- length(vars(x))
xpos <- lapply(xfix,function(y) which(regfix(x)$labels==y))
colpos <- lapply(xpos, function(y) ceiling(y/nrow))
rowpos <- lapply(xpos, function(y) (y-1)%%nrow+1)
for (i in seq_along(xfix))
regfix(x, from=vars(x)[rowpos[[i]]],to=vars(x)[colpos[[i]]]) <-
rep(colMeans(data[,xfix[i],drop=FALSE],na.rm=TRUE),length(rowpos[[i]]))
x <- updatelvm(x,zeroones=TRUE,deriv=TRUE)
}
ord <- c()
ordlist <- list()
for (i in seq_len(nrow(val$patterns))) {
ordlist <- c(ordlist, list(which(val$mis.type==i)))
ord <- c(ord, ordlist[[i]])
}
res <- with(val, list(coef=cc,
patterns=patterns, table=table(mis.type),
mis.type=mis.type,
order=ord,
orderlist=ordlist,
nmis=nmis,
allmis=pattern.allmis,
cc=mycc,
ncc=as.numeric(table(mis.type)[pattern.allcomp]),
multigroup=e.mis$model,
estimate=e.mis,
model=x,
model0=x0,
vcov=e.mis$vcov, opt=e.mis$opt,
control=control,
data=list(model.frame=data),
estimator=estimator,
call=cl
))
class(res) <- c("lvm.missing","lvmfit")
if (inherits(e.mis,"lvmfit.randomslope"))
class(res) <- c(class(res),"lvmfit.randomslope")
if (hessian & is.null(cluster)) {
if (messages>1)
message("Calculating asymptotic variance...\n")
res$vcov <- solve(information(res$estimate,type="hessian"))
cc[] <- coef(e.mis,type=1,vcov=res$vcov)
res$coef <- cc
}
return(res)
}
###}}} estimate.MAR.lvm
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.