gng.fit: Function for Fitting GNG model parameters

In DIME: Differential Identification using Mixture Ensemble

Function for Fitting GNG model parameters

Description

Function to estimate parameters for GNG model, mixture of exponential and k-normal. Parameters are estimated using EM algorithm.

Usage

gng.fit(data, avg = NULL, K = 2, weights = NULL, weights.cutoff = -1.345, 
  pi = NULL, mu = NULL, sigma = NULL, beta = NULL, tol = 1e-05, 
  max.iter = 2000, th = NULL)

Arguments

data	an R list of vector of normalized intensities (counts). Each element can correspond to particular chromosome. User can construct their own list containing only the chromosome(s) they want to analyze.
avg	optional vector of mean data (or log intensities). Only required when any one of huber weight (lower, upper or full) is selected.
K	optional number of normal component that will be fitted in GNG model.
weights	optional weights to be used for robust fitting. Can be a matrix the same length as data, or a character description of the huber weight method to be employed: "lower" - only value below weights.cutoff are weighted,\ "upper" - only value above weights.cutoff are weighted,\ "full" - both values above and below weights.cutoff are weighted,\ If selected, mean of data (avg) is required.
weights.cutoff	optional cutoff to be used with the Huber weighting scheme.
pi	optional vector containing initial estimates for proportion of the GNG mixture components. The first and last entries are for the estimates of negative and positive exponentials, respectively. The middle k entries are for normal components.
mu	optional vector containing initial estimates of the Gaussian means in GNG model.
sigma	optional vector containing initial estimates of the Gaussian standard deviation in GNG model. Must have K entries.
beta	optional vector containing initial estimates for the shape parameter in exponential components in GNG model. Must have 2 entries, one for negative exponential the other for positive exponential components.
tol	optional threshold for convergence for EM algorithm to estimate GNG parameters.
max.iter	optional maximum number of iterations for EM algorithm to estimate GNG parameters.
th	optional location parameter used to fit the negative and positive exponential model.

Value

A list of object:

name	the name of the model "GNG"
pi	a vector of estimated proportion of each components in the model
mu	a vector of estimated Gaussian means for k-normal components.
sigma	a vector of estimated Gaussian standard deviation for k-normal components.
beta	a vector of estimated exponential shape values.
th1	negative location parameter used to fit the negative exponential model.
th2	positive location parameter used to fit the positive exponential model.
K	the number of normal components in the corresponding mixture model.
loglike	the log likelihood for the fitted mixture model.
iter	the actual number of iterations run by the EM algorithm.
fdr	the local false discover rate estimated based on GNG model.
phi	a matrix of estimated GNG mixture component function.
AIC	Akaike Information Criteria.
BIC	Bayesian Information Criteria.

Author(s)

Cenny Taslim [e-mail], with contributions from Abbas Khalili [e-mail], Dustin Potter [e-mail], and Shili Lin [e-mail]

See Also

DIME, inudge.fit, nudge.fit

Examples

library(DIME)
# generate simulated datasets with underlying exponential-normal components
N1 <- 1500; N2 <- 500; K <- 4; rmu <- c(-2.25,1.50); rsigma <- c(1,1); 
rpi <- c(.05,.45,.45,.05); rbeta <- c(12,10);
set.seed(1234)
chr1 <- c(-rgamma(ceiling(rpi[1]*N1),shape = 1,scale = rbeta[1]), 
  rnorm(ceiling(rpi[2]*N1),rmu[1],rsigma[1]), 
  rnorm(ceiling(rpi[3]*N1),rmu[2],rsigma[2]), 
  rgamma(ceiling(rpi[4]*N1),shape = 1,scale = rbeta[2]));
chr2 <- c(-rgamma(ceiling(rpi[1]*N2),shape = 1,scale = rbeta[1]), 
  rnorm(ceiling(rpi[2]*N2),rmu[1],rsigma[1]), 
  rnorm(ceiling(rpi[3]*N2),rmu[2],rsigma[2]), 
  rgamma(ceiling(rpi[4]*N2),shape = 1,scale = rbeta[2])); 
chr3 <- c(-rgamma(ceiling(rpi[1]*N2),shape = 1,scale = rbeta[1]), 
  rnorm(ceiling(rpi[2]*N2),rmu[1],rsigma[1]), 
  rnorm(ceiling(rpi[3]*N2),rmu[2],rsigma[2]), 
  rgamma(ceiling(rpi[4]*N2),shape = 1,scale = rbeta[2]));
# analyzing only chromosome 1 and chromosome 3
data <- list(chr1,chr3);

# fit GNG model with 2 normal components
test <- gng.fit(data, K = 2);

# Getting the best fitted GNG model (parameters)
test$pi # estimated proportion of each component in GNG
test$mu # estimated mean of the normal component(s) GNG
# estimated standard deviation of the normal component(s) in GNG
test$sigma 
# estimated shape parameter of the exponential components in best model  
test$beta

DIME documentation built on May 9, 2022, 5:05 p.m.