normgam.fit: Normal-gamma Maximum Likelihood Estimator

In NormalGamma: Normal-gamma convolution model

Description

Computes the Maximum Likelihood Estimator for the normal-gamma distribution, either from a normal-gamma distributed sample or from two samples respectively normal-gamma and normally distributed.

Usage

normgam.fit(X, N = NULL, par.init = NULL, lower = NULL, upper = NULL,
            control = NULL, verbose = FALSE)

Arguments

X	vector of normal-gamma distributed values.
N	vector of normal distributed values.
par.init	vector of initial values for parameters (optional). par.init[1] and par.init[2] are the mean and standard deviation of the normal distribution, and par.init[3] and par.init[4] are the shape and scale parameters of the gamma distribution. See details for default initial values.
lower, upper	Bounds on the variables for maximization (optional).
control	list of control parameters (see details).
verbose	logical; if TRUE initial values of the parameters are printed.

Details

Likelihood maximization is run by the R function optimx.

By default, maximization is run with the following control parameters: the maximum number of iterations is equal to 1000 and the vector of scaling values for the parameters is (par0[1], par0[2], par0[3]*par0[4], sqrt(par0[3])*par0[4])/10 where par0 is the vector of default initial parameters. In case of unsuccessful convergence, maximization is run with optimx default control parameters. A list of control parameters can also be chosen by the user (see optimx).

If par.init == NULL, the initial parameters are computed in two ways depending if N is NULL or not. If N != NULL, the initial parameters are computed following the method of the moments (see Plancade S., Rozenholc Y. and Lund E., BMC Bioinfo 2012). If N == NULL, the initial parameters (par0[1], par0[2]) of the normal distribution are computed following the RMA procedure of Xie Y., Wang X. and Story M. (2009) for the normal-exponential convolution model, and the initial parameters of the gamma distribution, computed following the method of the moments, are (par0[4]=(sd(X)^2-par0[2])/(mean(X)-par0[1]), par0[3]=(mean(X)-par0[1]/par0[4]). Note that the RMA procedure for initial parameter computation when N == NULL stems from an heuristic adapted to microarray data. For parameters with different magnitude, user should specify initial parameters.

Value

par	vector of estimated parameters; par[1] and par[2] are the mean and standard deviation of the normal distribution and par[3] and par[4] are the shape and scale parameters of the gamma distribution
lik	value of the normal-gamma log-likelihood corresponding to par.
conv	integer code: 0 indicates successful convergence. This parameter has the value of the output parameter conv from the procedure optimx used for likelihood maximization (see optimx for details).

Author(s)

Plancade S. and Rozenholc Y.

References

Plancade S., Rozenholc Y. and Lund E. "Generalization of the normal-exponential model : exploration of a more accurate parametrisation for the signal distribution on Illumina BeadArrays", BMC Bioinfo 2012, 13(329).

Xie Y. Wang X. and Story M. (2009), Statistical methods of background correction for Illumina BeadArray data, Bioinformatics, 25(6), 751–757.

See Also

dnormgam computes the density of the normal-gamma distribution and normgam.signal implements the background correction using the normal-gamma model. Intensities provides an example of Illumina microarray data.

Examples

# Example 1: simulated data

## Not run: 

n = 1000
par = c(60,5,0.15,400)
X = rnorm(n, mean=par[1], sd=par[2]) + rgamma(n, shape=par[3], scale=par[4]) 
N = rnorm(100, mean=par[1], sd=par[2])

par1 = normgam.fit(X, N)$par
par2 = normgam.fit(X)$par


F1 = dnormgam(par1, plot=FALSE)
F2 = dnormgam(par2, plot=FALSE)

par(mfrow=c(2,1))

H = histogram(X, type='irregular', verbose=FALSE, plot=FALSE)

plot(H, xlim=c(0,500))
lines(F1$xout, F1$dout,col='red')

plot(H, xlim=c(0,500))
lines(F2$xout, F2$dout,col='blue')

## End(Not run)

# Example 2: Illumina data

## Not run: 

data(RegNegIntensities_Example)
 
X = Intensities$Regular
N = Intensities$Negative

par1 = normgam.fit(X, N)$par
par2 = normgam.fit(X)$par

F1 = dnormgam(par1, plot=FALSE)
F2 = dnormgam(par2, plot=FALSE)

par(mfrow=c(2,1))

H = histogram(X, type='irregular', verbose=FALSE, plot=FALSE)

plot(H, xlim=c(0,500))
lines(F1$xout, F1$dout, col='red')

plot(H, xlim=c(0,500))
lines(F2$xout, F2$dout, col='blue')

## End(Not run)

NormalGamma documentation built on May 2, 2019, 10:25 a.m.