Likelihood: Compute the log-likelihood of a parameter

In ElstonStewart: Elston-Stewart Algorithm

Compute the log-likelihood of a parameter

Description

Compute the log-likelihood of a parameter theta, given a list of pedigrees and a modele, using multiple cores and memoization

Usage

Likelihood(ped.set, modele, theta, 
                  n.cores = getOption("mc.cores", 2L), 
                  optim.alloc = TRUE, sum.likelihoods = TRUE,
                  PSOCK = Sys.info()["sysname"] != "Linux")

Arguments

ped.set	a list of pedigrees, each created with es.pedigree
modele	a modele
theta	a parameter for the modele
n.cores	number of cores on which to run the computation
optim.alloc	cf Details
sum.likelihoods	cf Value
PSOCK	Use PSOCK cluster instead of fork cluster (defaults to TRUE on non-Linux)

Details

The parameter theta will be given to the functions in modele to compute the likelihoods.

If n.cores > 1 a cluster is created and left open for futur use with the same parameters ped.set and modele. Open clusters can be closed with es.stopCluster().

If optim.alloc = TRUE, the function tries to optimize the distribution of the computation (if n.cores > 1) between the cluster nodes. This has no effect on the first run for given parameters ped.set and modele, but will reduce running time if the algorithm is ran several times with different values of theta and same parameters ped.set and modele. This is typically usefull for likelihood maximization.

Value

If sum.likelihoods = TRUE, the function returns a single value, the sum of the log-likelihoods computed for each pedigree of ped.set. Else, the function returns a vector containing these log-likelihoods.

See Also

Elston, es.pedigree, es.stopCluster

Examples

data(fams)

# this data frame contains various families
# getting their famid
fam.ids <- unique(fams$fam);

# creating a list of genotypes corresponding to all individuals in fams
# NA -> 0, 1 or 2
genotypes <- lapply( fams$genotype, function(x) if(is.na(x)) 0:2 else x )

# creating a list of es.pedigree objects
X <- vector("list", length(fam.ids))
for(i in seq_along(fam.ids))
{
  w <- which(fams$fam == fam.ids[i])
  X[[i]] <- es.pedigree( id = fams$id[w], father = fams$father[w],
      mother = fams$mother[w], sex = fams$sex[w], pheno = rep(0, length(w)), 
      geno = genotypes[w], famid = fam.ids[i] )
}

## Not run: # computing the likelihood for a single value p
Likelihood(X, modele.di, theta = list( p=0.5), n.cores=1 )

# computing the likelihood for a vector p (Elston-Stewart is ran only once!)
p <- seq(0,1,length=501)
L <- Likelihood(X, modele.di, theta = list( p=p ), n.cores=1 ) 
plot( p, exp(L), type="l")

# running an optimization algorithm
# Elston-Stewart is ran several times
# here we run the algorithm with 2 cores
L <- function(p) -Likelihood(X, modele.di, theta = list( p=p ), n.cores=2 ) 
optimize(L , c(0.35,0.45) )
## End(Not run)

ElstonStewart documentation built on June 20, 2022, 5:15 p.m.