R/lca.model.R

In LCAextend: Latent Class Analysis (LCA) with Familial Dependence in Extended Pedigrees

lca.model <-
function(ped,probs,param,optim.param,fit=TRUE,optim.probs.indic=c(TRUE,TRUE,TRUE,TRUE),tol=0.001,x=NULL,var.list=NULL,famdep=TRUE,modify.init=NULL)
{
    fam <- ped[,1]
    id <- ped[,2]
    dad <- ped[,3]
    mom <- ped[,4]
    sex <- ped[,5]
    status <- ped[,6]

	y <- as.matrix(ped[,-(1:6)])

    K <- length(probs$p)

	if(famdep==TRUE&(is.null(probs$p)|is.null(probs$p0)|is.null(probs$p.trans))) stop("One component of probs is missing for the familial dependence case.")
	
	optim.param <- attrib.dens(optim.param)
	if(is.null(attr(optim.param,"type"))) stop("Function optim.param must have an attribute of name \'type\'.")
	else
	{
		if(attr(optim.param,"type")=="norm") dens <- dens.norm
		if(attr(optim.param,"type")=="ordi") 
		{
			dens <- dens.prod.ordi
			for (j in 1:ncol(y))
			{
				if (!is.factor(y[,j]) & !is.integer(y[,j])) stop("Symptom ",j," must be of type integer or factor with multinomial distribution.")
				my = min(y[,j],na.rm=T)
				if (my<1) 
				{
				  # On doit modifier ped aussi car il est passe a la fonction e.step
				  ped[,6+j]	= y[,j] = as.integer(y[,j] + 1 - my)
					warning("Values of symptom ",j," shifted by ",1-my," to set the minimum to 1.")
				}
			}
		}
			
	}
	
	if(famdep==TRUE&is.null(probs$p.found))
	{
		#if all founds are affected
		if(all(status[dad==0]==2)) probs$p.found <- 1
		else
		{
			#if all founds are unaffected
			if(all(status[dad==0]==1)) probs$p.found <- 0
			else
			#if at least two different statuses for founders
			stop("component p.found of probs is needed")
		}
	}
	
	if(famdep==TRUE&is.null(probs$p.child))
	{
		#if all children are affected
		if(all(status[dad>0]==2)) probs$p.child <- 1
		else
		{
			#if all children are unaffected
			if(all(status[dad>0]==1)) probs$p.child <- 0
			else
			#if at least two different statuses for children
			stop("component p.child of probs is needed")
		}
	}

	if(famdep==FALSE&(is.null(probs$p)|is.null(probs$p0)|is.null(probs$p.aff))) stop("one component of probs is missing, case indpendent individuals")

	children <- which(dad>0)
	connects <- intersect(children,union(which(id%in%dad),which(id%in%mom)))
	if(length(connects)>0&famdep==TRUE&is.null(probs$p0connect)) stop("p0connect must be given when analysing extended pedigrees")

	optim.param <- attrib.dens(optim.param)
	if(is.null(attr(optim.param,"type"))) stop("function optim.param must have an attribute of name \'type\'")
	else
	{
		if(attr(optim.param,"type")=="ordi") dens <- dens.prod.ordi
		if(attr(optim.param,"type")=="norm") dens <- dens.norm
	}

    if(any(dad!=0&mom==0)|any(dad==0&mom!=0)) stop("there is individuals with a single parent\n")
    if(all(dad==0)&all(mom==0)) stop("it needs at least one child for computing p.trans\n")
    if(length(id)<ncol(y)) warning("not enough individuals to have invertible variance-covariance matrix\n")

	if(any(status[!id%in%c(dad,mom)]!=2)) stop("All children without descendants must have a symptom status = 2")

    peel <- list()
    for(f in unique(fam))
    {
        id.fam <- id[fam==f]
        dad.fam <- dad[fam==f]
        mom.fam <- mom[fam==f]
        sex.fam <- sex[fam==f]
        status.fam <- status[fam==f]

        couple <- cbind(unique(dad.fam[dad.fam>0]),unique(mom.fam[mom.fam>0]))
        couple <- rbind(couple,cbind(couple[,2],couple[,1]))
        # Changement par rapport aux versions 1.1 et ant?rieur pour s'adapter au module kinship2
        depth <- kindepth(id.fam,dad.fam,mom.fam,align=TRUE)
        generation <- max(depth)
        out <- NULL
        for(generat in 1:generation) out[generat] <- sum(depth==generat-1)/2
        peel.connect <- matrix(0,nrow=generation,ncol=max(out))
        peel.connect[generation,1] <- id.fam[depth==0][1]
        if(generation>1) for(generat in 1:(generation-1))
        {
            connect.generat <- intersect(id.fam[depth==generat][(id.fam[depth==generat])&(dad.fam[depth==generat]>0)],union(dad.fam,mom.fam))
            peel.connect[generation-generat,1:length(connect.generat)] <- connect.generat
        }
        peel[[which(unique(fam)==f)]] <- list("generation"=generation,"peel.connect"=peel.connect,"couple"=couple)
    }
    cat("CHECK OF ALL PEDIGREES DONE\n")

    if(any(probs$p<0)|abs(sum(probs$p)-1)>.Machine$double.eps) stop("There is a problem with the initial value of founder class probability p.\n")
    if(famdep) if(any(probs$p.trans<0)|any(abs(apply(array(probs$p.trans[,,1:K],dim=rep(K,times=3)),2:3,sum)-1)>.Machine$double.eps)|
       any(abs(apply(matrix(probs$p.trans[,1:K,K+1],nrow=K,ncol=K),2,sum)-1)>.Machine$double.eps))
    stop("There is a problem with the initial value of transition probability p.trans.\n")

	y.aff <- as.matrix(y[status==2,])
	y.x.aff = y.aff
    if(!is.null(x)) y.x.aff <- cbind(y.x.aff,as.matrix(x[status==2,]))

    fyc <- matrix(1,nrow=length(id),ncol=K+1)
    fyc[status==2,1:K] <- t(apply(y.x.aff,1,dens,param,var.list))
    if(K>1) prob0 <- apply(fyc[status==2,1:K],2,function(vec) all(vec<.Machine$double.eps))
    else prob0 <- fyc[status==2,1]<.Machine$double.eps
    if(fit&any(prob0))
    {
        if(!is.null(modify.init)) while(any(prob0))
        {
            param <- modify.init(prob0,param)
            fyc[status==2,1:K] <- t(apply(y.x.aff,1,dens,param,var.list))
            if(K>1) prob0 <- apply(fyc[status==2,1:K],2,function(vec) all(vec<.Machine$double.eps))
            else prob0 <- fyc[status==2,1]<.Machine$double.eps
        }
        else stop("initial parameter values produce zero probabilities for at least one subject.\n")
    }

    iter <- n.warning <- 0
    ll <- -Inf
    continue <- 1
    while(continue)
    {
        iter <- iter+1
        cat("%%%%%%%%%%%%%%%%%%%%%%%%%%% iteration",iter,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n")
        res.weight <- e.step(ped,probs,param,dens,peel,x,var.list,famdep)
        ll.new <- res.weight$ll
        if(!fit) continue <- 0
        else
        {
            probs <- optim.probs(ped,probs,optim.probs.indic,res.weight,famdep)

            weight <- as.matrix(res.weight$w[,1,1:K])

            param <- optim.param(y.aff,status,weight,param,x=x,var.list=var.list)

            if(ll.new<ll)
            {
                cat("log-likelihood is decreasing!!!!!\n")
                n.warning <- n.warning+1
            }
            else if(ll.new-ll<tol|n.warning>=20) continue <- 0
        }
        ll <- ll.new
    }
    res <- list("probs"=probs,"param"=param,"weight"=res.weight$w,"ll"=ll)
    res
}