Version:	2.0.0.1
Date:	2022-12-04
Title:	Tools for Analyzing Finite Mixture Models
Depends:	R (≥ 4.0.0)
Imports:	kernlab, MASS, plotly, scales, segmented, stats, survival
URL:	https://github.com/dsy109/mixtools
Description:	Analyzes finite mixture models for various parametric and semiparametric settings. This includes mixtures of parametric distributions (normal, multivariate normal, multinomial, gamma), various Reliability Mixture Models (RMMs), mixtures-of-regressions settings (linear regression, logistic regression, Poisson regression, linear regression with changepoints, predictor-dependent mixing proportions, random effects regressions, hierarchical mixtures-of-experts), and tools for selecting the number of components (bootstrapping the likelihood ratio test statistic, mixturegrams, and model selection criteria). Bayesian estimation of mixtures-of-linear-regressions models is available as well as a novel data depth method for obtaining credible bands. This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772 and the Chan Zuckerberg Initiative: Essential Open Source Software for Science (Grant No. 2020-255193).
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
NeedsCompilation:	yes
Packaged:	2025-03-08 07:39:45 UTC; hornik
Author:	Derek Young [aut, cre], Tatiana Benaglia [aut], Didier Chauveau [aut], David Hunter [aut], Kedai Cheng [aut], Ryan Elmore [ctb], Thomas Hettmansperger [ctb], Hoben Thomas [ctb], Fengjuan Xuan [ctb]
Maintainer:	Derek Young <derek.young@uky.edu>
Repository:	CRAN
Date/Publication:	2025-03-08 08:58:50 UTC

GNP and CO2 Data Set

Description

This data set gives the gross national product (GNP) per capita in 1996 for various countries as well as their estimated carbon dioxide (CO2) emission per capita for the same year.

Usage

data(CO2data)

Format

This data frame consists of 28 countries and the following columns:

GNP

The gross national product per capita in 1996.

CO2

The estimated carbon dioxide emission per capita in 1996.

country

An abbreviation pertaining to the country measured (e.g., "GRC" = Greece and "CH" = Switzerland).

References

Hurn, M., Justel, A. and Robert, C. P. (2003) Estimating Mixtures of Regressions, Journal of Computational and Graphical Statistics 12(1), 55–79.

Infant habituation data

Description

From Thomas et al (2011):

"Habituation is a standard method of studying infant behaviors. Indeed, much of what is known about infant memory and perception rests on habituation methods. Six-month infants (n = 51) were habituated to a checker-board pattern on two occasions, one week apart. On each occasion, the infant was presented with the checkerboard pattern and the length of time the infant viewed the pattern before disengaging was recorded; this denoted the end of a trial. After disengagement, another trial was presented. The procedure was implemented for eleven trials. The conventional index of habituation performance is the summed observed fixation to the checkerboard pattern over the eleven trials. Thus, an index of reliability focuses on how these fixation times, in seconds, on the two assessment occasions correlate: r = .29."

Usage

data(Habituationdata)

Format

A data frame with two variables, m1 and m2, and 51 cases. The two variables are the summed observations times for the two occasions described above.

Author(s)

Hoben Thomas

Source

Original source: Thomas et al. (2011). See references section.

References

Thomas, H., Lohaus, A., and Domsch, H. (2011), Extensions of Reliability Theory, in Nonparametric Statistics and Mixture Models: A Festschrift in Honor of Thomas Hettmansperger (Singapore: World Scientific), pp. 309-316.

Ethanol Fuel Data Set

Description

This data set gives the equivalence ratios and peak nitrogen oxide emissions in a study using pure ethanol as a spark-ignition engine fuel.

Usage

data(NOdata)

Format

This data frame consists of:

NO

The peak nitrogen oxide emission levels.

Equivalence

The equivalence ratios for the engine at compression ratios from 7.5 to 18.

Source

Brinkman, N. D. (1981) Ethanol Fuel – A Single-Cylinder Engine Study of Efficiency and Exhaust Emissions, S.A.E. Transactions, 68.

References

Hurn, M., Justel, A. and Robert, C. P. (2003) Estimating Mixtures of Regressions, Journal of Computational and Graphical Statistics 12(1), 55–79.

Reaction Time (RT) Data Set

Description

This data set involves normally developing children 9 years of age presented with two visual simuli on a computer monitor. The left image is the target stimuli and on the right is either an exact copy or a mirror image of the target stimuli. The child must press one key if it is a copy or another key if it is a mirror image. The data consists of the reaction times (RT) of the 197 children who provided correct responses for all 6 task trials.

Usage

data(RTdata)

Format

This data frame consists of 197 children (the rows) and their 6 responses (the columns) to the stimulus presented. The response (RT) is recorded in milliseconds.

References

Cruz-Medina, I. R., Hettmansperger, T. P. and Thomas, H. (2004) Semiparametric Mixture Models and Repeated Measures: The Multinomial Cut Point Model, Applied Statistics 53(3), 463–474.

Miller, C. A., Kail, R., Leonard, L. B. and Tomblin, J. B. (2001) Speed of Processing in Children with Specific Language Impairment, Journal of Speech, Language, and Hearing Research 44(2), 416–433.

See Also

RTdata2

Reaction Time (RT) Data Set (No. 2)

Description

This data set involves normally developing children 9 years of age presented visual simuli on a computer monitor. There are three different experimental conditions, according to the length of the delay after which the stimulus was displayed on the screen. Each subject experienced each condition eight times, and these 24 trials were given in random order. These data give the 82 children for whom there are complete measurements among over 200 total subjects.

Usage

data(RTdata2)

Format

This data frame consists of 82 children (the rows) and their 24 responses (the columns) to the stimulus presented. The response is recorded in milliseconds. The columns are not in the order in which the stimuli were presented to the children; rather, they are arranged into three blocks of eight columns each so that each eight-column block contains only trials from one of the three conditions.

References

Miller, C. A., Kail, R., Leonard, L. B. and Tomblin, J. B. (2001) Speed of Processing in Children with Specific Language Impairment, Journal of Speech, Language, and Hearing Research 44(2), 416–433.

See Also

RTdata

Simulated Data from 2-Component Mixture of Regressions with Random Effects

Description

This data set was generated from a 2-component mixture of regressions with random effects.

Usage

data(RanEffdata)

Format

This data set consists of a list with 100 25x3 matrices. The first column is the response variable, the second column is a column of 1's and the last column is the predictor variable.

See Also

regmixEM.mixed

Rod and Frame Task Data Set

Description

This data set involves assessing children longitudinally at 6 age points from ages 4 through 18 years for the rod and frame task. This task sits the child in a darkened room in front of a luminous square frame tilted at 28 degrees on its axis to the left or right. Centered inside the frame was a luminous rod also tilted 28 degrees to the left or right. The child's task was to adjust the rod to the vertical position and the absolute deviation from the vertical (in degrees) was the measured response.

Usage

data(RodFramedata)

Format

This data frame consists of 140 children (the rows). Column 1 is the subject number and column 2 is the sex (0=MALE and 1=FEMALE). Columns 3 through 26 give the 8 responses at each of the ages 4, 5, and 7. Columns 27 through 56 give the 10 responses at each of the ages 11, 14, and 18. A value of 99 denotes missing data.

Source

Thomas, H. and Dahlin, M. P. (2005) Individual Development and Latent Groups: Analytical Tools for Interpreting Heterogeneity, Developmental Review 25(2), 133–154.

Water-Level Task Data Set

Description

This data set arises from the water-level task proposed by the Swiss psychologist Jean Piaget to assess children's understanding of the physical world. This involves presenting a child with a rectangular vessel with a cap, affixed to a wall, that can be tilted (like the minute hand of a clock) to point in any direction. A separate disk with a water line indicated on it, which can similarly be spun so that the water line may assume any desired angle with the horizontal, is positioned so that by spinning this disk, the child subject may make the hypothetical surface of water inside the vessel assume any desired orientation. For each of eight different orientations of the vessel, corresponding to the clock angles at 1:00, 2:00, 4:00, 5:00, 7:00, 8:00, 10:00, and 11:00, the child subject is asked to position the water level as it would appear in reality if water were in the vessel. The measurement is the acute angle with the horizontal, in degrees, assumed by the water line after it is positioned by the child. A sign is attached to the measurement to indicate whether the line slopes up (positive) or down (negative) from left to right. Thus, each child has 8 repeated measurements, one for each vessel angle, and the range of possible values are from -90 to 90.

The setup of the experiment, along with a photograph of the testing apparatus, is given by Thomas and Jamison (1975). A more detailed analysis using a subset of 405 of the original 579 subjects is given by Thomas and Lohaus (1993); further analyses using the functions in mixtools are given by Benaglia et al (2008) and Levine et al (2011), among others.

There are two versions of the dataset included in mixtools. The full dataset, called WaterdataFull, has 579 individuals. The dataset called Waterdata is a subset of 405 individuals, comprising all children aged 11 years or more and omitting any individuals with any observations equal to 100, which in this context indicates a missing value (since all of the degree measurements should be in the range from -90 to +90, 100 is not a possible value).

Usage

data(Waterdata)

Format

These data frames consist of 405 or 579 rows, one row for each child. There are ten columns: The age (in years) and sex (where 1=male and 0=female) are given for each individual along with the degree of deviation from the horizontal for 8 specified clock-hour orientations (11, 4, 2, 7, 10, 5, 1, and 8 o'clock, in order).

Source

Benaglia, T., Chauveau, D., and Hunter, D.R. (2009), An EM-Like Algorithm for Semi- and Non-Parametric Estimation in Multivariate Mixtures, Journal of Computational and Graphical Statistics, 18: 505-526.

Levine, M., Hunter, D.R., and Chauveau, D. (2011), Maximum Smoothed Likelihood for Multivariate Mixtures, Biometrika, 98(2): 403-416.

Thomas, H. and Jamison, W. (1975), On the Acquisition of Understanding that Still Water is Horizontal, Merrill-Palmer Quarterly of Behavior and Development, 21(1): 31-44.

Thomas, H. and Lohaus, A. (1993), Modeling Growth and Individual Differences in Spatial Tasks, University of Chicago Press, Chicago, available on JSTOR.

Augmented Predictor Function

Description

Creates the augmented predictor matrix based on an appropriately defined changepoint structure.

Usage

aug.x(X, cp.locs, cp, delta = NULL)

Arguments

Details

This function is called by segregmixEM and the associated internal functions.

Value

aug.x returns a matrix of the original matrix X with the predictor adjusted for changepoints and (optional) discontinuities.

See Also

segregmixEM

Examples

x <- matrix(1:30, nrow = 10)
cp <- c(1, 3, 0)
cp.locs <- c(3, 12, 14, 16)
d <- rep(0, 4)
x1 <- aug.x(x, cp.locs, cp, delta = NULL)
x1
x2 <- aug.x(x, cp.locs, cp, delta = d)
x2

Performs Parametric Bootstrap for Sequentially Testing the Number of Components in Various Mixture Models

Description

Performs a parametric bootstrap by producing B bootstrap realizations of the likelihood ratio statistic for testing the null hypothesis of a k-component fit versus the alternative hypothesis of a (k+1)-component fit to various mixture models. This is performed for up to a specified number of maximum components, k. A p-value is calculated for each test and once the p-value is above a specified significance level, the testing terminates. An optional histogram showing the distribution of the likelihood ratio statistic along with the observed statistic can also be produced.

Usage

boot.comp(y, x = NULL, N = NULL, max.comp = 2, B = 100,
          sig = 0.05, arbmean = TRUE, arbvar = TRUE,
          mix.type = c("logisregmix", "multmix", "mvnormalmix",
          "normalmix", "poisregmix", "regmix", "regmix.mixed", 
          "repnormmix"), hist = TRUE, ...)

Arguments

Value

boot.comp returns a list with items:

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

See Also

logisregmixEM, multmixEM, mvnormalmixEM, normalmixEM, poisregmixEM, regmixEM, regmixEM.mixed, repnormmixEM

Examples

## Bootstrapping to test the number of components on the RTdata.

data(RTdata)
set.seed(100)
x <- as.matrix(RTdata[, 1:3])
y <- makemultdata(x, cuts = quantile(x, (1:9)/10))$y
a <- boot.comp(y = y, max.comp = 1, B = 5, mix.type = "multmix", 
               epsilon = 1e-3)
a$p.values

Performs Parametric Bootstrap for Standard Error Approximation

Description

Performs a parametric bootstrap by producing B bootstrap samples for the parameters in the specified mixture model.

Usage

boot.se(em.fit, B = 100, arbmean = TRUE, arbvar = TRUE, 
        N = NULL, ...)

Arguments

Value

boot.se returns a list with the bootstrap samples and standard errors for the mixture of interest.

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

Examples

## Bootstrapping standard errors for a regression mixture case.

data(NOdata)
attach(NOdata)
set.seed(100)
em.out <- regmixEM(Equivalence, NO, arbvar = FALSE)
out.bs <- boot.se(em.out, B = 10, arbvar = FALSE)
out.bs

Plot the Component CDF

Description

Plot the components' CDF via the posterior probabilities.

Usage

compCDF(data, weights, 
        x=seq(min(data, na.rm=TRUE), max(data, na.rm=TRUE), len=250), 
        comp=1:NCOL(weights), makeplot=TRUE, ...) 

Arguments

Details

When makeplot is TRUE, a line plot is produced of the CDFs evaluated at x. The plot is not a step function plot; the points (x, CDF(x)) are simply joined by line segments.

Value

A matrix with length(comp) rows and length(x) columns in which each row gives the CDF evaluated at each point of x.

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

Elmore, R. T., Hettmansperger, T. P. and Xuan, F. (2004) The Sign Statistic, One-Way Layouts and Mixture Models, Statistical Science 19(4), 579–587.

See Also

makemultdata, multmixmodel.sel, multmixEM.

Examples

## The sulfur content of the coal seams in Texas

set.seed(100)

A <- c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17)
B <- c(1.69, 0.64, .9, 1.41, 1.01, .84, 1.28, 1.59) 
C <- c(1.56, 1.22, 1.32, 1.39, 1.33, 1.54, 1.04, 2.25, 1.49) 
D <- c(1.3, .75, 1.26, .69, .62, .9, 1.2, .32) 
E <- c(.73, .8, .9, 1.24, .82, .72, .57, 1.18, .54, 1.3)

dis.coal <- makemultdata(A, B, C, D, E, 
                         cuts = median(c(A, B, C, D, E)))
temp <- multmixEM(dis.coal)

## Now plot the components' CDF via the posterior probabilities

compCDF(dis.coal$x, temp$posterior, xlab="Sulfur", ylab="", main="empirical CDFs")

Density Function for the Dirichlet Distribution

Description

Density function for the Dirichlet distribution.

Usage

ddirichlet(x, alpha)

Arguments

Details

This is usually not to be called by the user.

Normal kernel density estimate for nonparametric EM output

Description

Takes an object of class npEM and returns an object of class density giving the kernel density estimate for the selected component and, if applicable, the selected block.

Usage

## S3 method for class 'npEM'
density(x, u=NULL, component=1, block=1, scale=FALSE, ...)

Arguments

Details

The bandwidth is taken to be the same as that used to produce the npEM object, which is given by x$bandwidth.

Value

density.npEM returns a list of type "density". See density for details. In particular, the output of density.npEM may be used directly by functions such as plot or lines.

See Also

npEM, spEMsymloc, plot.npEM

Examples

 
## Look at histogram of Old Faithful waiting times
data(faithful)
Minutes <- faithful$waiting
hist(Minutes, freq=FALSE)

## Superimpose equal-variance normal mixture fit:
set.seed(100)
nm <- normalmixEM(Minutes, mu=c(50,80), sigma=5, arbvar=FALSE, fast=TRUE)
x <- seq(min(Minutes), max(Minutes), len=200)
for (j in 1:2) 
  lines(x, nm$lambda[j]*dnorm(x, mean=nm$mu[j], sd=nm$sigma), lwd=3, lty=2)
  
## Superimpose several semiparametric fits with different bandwidths:
bw <- c(1, 3, 5)
for (i in 1:3) {
  sp <- spEMsymloc(Minutes, c(50,80), bw=bw[i], eps=1e-3)
  for (j in 1:2) 
    lines(density(sp, component=j, scale=TRUE), col=1+i, lwd=2)    
}
legend("topleft", legend=paste("Bandwidth =",bw), fill=2:4)

Normal kernel density estimate for semiparametric EM output

Description

Takes an object of class spEM and returns an object of class density giving the kernel density estimate.

Usage

## S3 method for class 'spEM'
density(x, u=NULL, component=1, block=1, scale=FALSE, ...)

Arguments

Details

The bandwidth is taken to be the same as that used to produce the npEM object, which is given by x$bandwidth.

Value

density.spEM returns a list of type "density". See density for details. In particular, the output of density.spEM may be used directly by functions such as plot or lines.

See Also

spEM, spEMsymloc, plot.spEM

Examples

 
set.seed(100)
mu <- matrix(c(0, 15), 2, 3)
sigma <- matrix(c(1, 5), 2, 3)
x <- rmvnormmix(300, lambda = c(.4,.6), mu = mu, sigma = sigma)

d <- spEM(x, mu0 = 2, blockid = rep(1,3), constbw = TRUE) 
plot(d, xlim=c(-10, 40), ylim = c(0, .16), xlab = "", breaks = 30, 
     cex.lab=1.5, cex.axis=1.5) # plot.spEM calls density.spEM here

Elliptical and Spherical Depth

Description

Computation of spherical or elliptical depth.

Usage

 depth(pts, x, Cx = var(x))

Arguments

Value

depth returns a k-vector where each entry is the elliptical depth of a point in pts.

Note

depth is used in regcr.

References

Elmore, R. T., Hettmansperger, T. P. and Xuan, F. (2000) Spherical Data Depth and a Multivariate Median, Proceedings of Data Depth: Robust Multivariate Statistical Analysis, Computational Geometry and Applications.

See Also

regcr

Examples

  set.seed(100)
  x <- matrix(rnorm(200),nc = 2)
  depth(x[1:3, ], x)

The Multivariate Normal Density

Description

Density and log-density for the multivariate normal distribution with mean equal to mu and variance matrix equal to sigma.

Usage

dmvnorm(y, mu=NULL, sigma=NULL)
logdmvnorm(y, mu=NULL, sigma=NULL)

Arguments

Details

This code is written to be efficient, using the qr-decomposition of the covariance matrix (and using it only once, rather than recalculating it for both the determinant and the inverse of sigma).

Value

dmvnorm gives the densities, while logdmvnorm gives the logarithm of the densities.

See Also

qr, qr.solve, dnorm, rmvnorm

Draw Two-Dimensional Ellipse Based on Mean and Covariance

Description

Draw a two-dimensional ellipse that traces a bivariate normal density contour for a given mean vector, covariance matrix, and probability content.

Usage

ellipse(mu, sigma, alpha = .05, npoints = 250, newplot = FALSE,
        draw = TRUE, ...)

Arguments

Value

ellipse returns an npointsx2 matrix of the points forming the border of the ellipse.

References

Johnson, R. A. and Wichern, D. W. (2002) Applied Multivariate Statistical Analysis, Fifth Edition, Prentice Hall.

See Also

regcr

Examples

## Produce a 95% ellipse with the specified mean and covariance structure. 

mu <- c(1, 3)
sigma <- matrix(c(1, .3, .3, 1.5), 2, 2)

ellipse(mu, sigma, npoints = 200, newplot = TRUE)
 

EM algorithm for Reliability Mixture Models (RMM) with right Censoring

Description

Parametric EM algorithm for univariate finite mixture of exponentials distributions with randomly right censored data.

Usage

  expRMM_EM(x, d=NULL, lambda = NULL, rate = NULL, k = 2, 
		    complete = "tdz", epsilon = 1e-08, maxit = 1000, verb = FALSE) 

Arguments

Value

expRMM_EM returns a list of class "mixEM" with the following items:

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: plotexpRMM, summary.mixEM.

Other models and algorithms for censored lifetime data: weibullRMM_SEM, spRMM_SEM.

Examples

n <- 300 # sample size
m <- 2   # number of mixture components
lambda <- c(1/3,1-1/3); rate <- c(1,1/10) # mixture parameters
set.seed(1234)
x <- rexpmix(n, lambda, rate) # iid ~ exponential mixture
cs <- runif(n,0,max(x)) # Censoring (uniform) and incomplete data
t <- apply(cbind(x,cs),1,min) # observed or censored data
d <- 1*(x <= cs)              # censoring indicator

###### EM for RMM, exponential lifetimes
l0 <- rep(1/m,m); r0 <- c(1, 0.5) # "arbitrary" initial values
a <- expRMM_EM(t, d, lambda = l0, rate = r0, k = m)
summary(a)                 # EM estimates etc
plotexpRMM(a, lwd=2) # default plot of EM sequences
plot(a, which=2) # or equivalently, S3 method for "mixEM" object

EM Algorithm for Mixtures of Regressions with Flare

Description

Returns output for 2-component mixture of regressions with flaring using an EM algorithm with one step of Newton-Raphson requiring an adaptive barrier for maximization of the objective function. A mixture of regressions with flare occurs when there appears to be a common regression relationship for the data, but the error terms have a mixture structure of one normal component and one exponential component.

Usage

flaremixEM(y, x, lambda = NULL, beta = NULL, sigma = NULL, 
           alpha = NULL, nu = NULL, epsilon = 1e-04, 
           maxit = 10000, verb = FALSE, restart = 50)

Arguments

Value

flaremixEM returns a list of class mixEM with items:

See Also

regmixEM

Examples

## Simulation output.

set.seed(100)
j=1
while(j == 1){
    x1 <- runif(30, 0, 10)
    x2 <- runif(20, 10, 20)
    x3 <- runif(30, 20, 30)
    y1 <- 3+4*x1+rnorm(30, sd = 1)
    y2 <- 3+4*x2+rexp(20, rate = .05)
    y3 <- 3+4*x3+rnorm(30, sd = 1)
    x <- c(x1, x2, x3)
    y <- c(y1, y2, y3)
    nu <- (1:30)/2

    out <- try(flaremixEM(y, x, beta = c(3, 4), nu = nu,
               lambda = c(.75, .25), sigma = 1), silent = TRUE)
    if(any(class(out) == "try-error")){
        j <- 1
    } else j <- 2
}

out[4:7]
plot(x, y, pch = 19)
abline(out$beta)

EM Algorithm for Mixtures of Gamma Distributions

Description

Return EM algorithm output for mixtures of gamma distributions.

Usage

gammamixEM(x, lambda = NULL, alpha = NULL, beta = NULL, k = 2,
           mom.start = TRUE, fix.alpha = FALSE, epsilon = 1e-08, 
           maxit = 1000, maxrestarts = 20, verb = FALSE)

Arguments

Value

gammamixEM returns a list of class mixEM with items:

References

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977) Maximum Likelihood From Incomplete Data Via the EM Algorithm, Journal of the Royal Statistical Society, Series B, 39(1), 1–38.

Young, D. S., Chen, X., Hewage, D., and Nilo-Poyanco, R. (2019) Finite Mixture-of-Gamma Distributions: Estimation, Inference, and Model-Based Clustering, Advances in Data Analysis and Classification, 13(4), 1053–1082.

Examples

##Analyzing a 3-component mixture of gammas.

set.seed(100)
x <- c(rgamma(200, shape = 0.2, scale = 14), rgamma(200, 
     shape = 32, scale = 10), rgamma(200, shape = 5, scale = 6))
out <- gammamixEM(x, lambda = c(1, 1, 1)/3, verb = TRUE)
out[2:4]

EM Algorithm for Mixtures-of-Experts

Description

Returns EM algorithm output for a mixture-of-experts model. Currently, this code only handles a 2-component mixture-of-experts, but will be extended to the general k-component hierarchical mixture-of-experts.

Usage

hmeEM(y, x, lambda = NULL, beta = NULL, sigma = NULL, w = NULL,
      k = 2, addintercept = TRUE, epsilon = 1e-08, 
      maxit = 10000, verb = FALSE)

Arguments

Value

hmeEM returns a list of class mixEM with items:

References

Jacobs, R. A., Jordan, M. I., Nowlan, S. J. and Hinton, G. E. (1991) Adaptive Mixtures of Local Experts, Neural Computation 3(1), 79–87.

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

See Also

regmixEM

Examples

## EM output for NOdata.
 
data(NOdata)
attach(NOdata)
set.seed(100)
em.out <- regmixEM(Equivalence, NO)
hme.out <- hmeEM(Equivalence, NO, beta = em.out$beta)
hme.out[3:7]

Integrated Squared Error for a selected density from npEM output

Description

Computes the integrated squared error for a selected estimated density from npEM output (selected by specifying the component and block number), relative to a true pdf that must be specified by the user. The range for the numerical integration must be specified. This function also returns (by default) a plot of the true and estimated densities.

Usage

ise.npEM(npEMout, component=1, block=1, truepdf, lower=-Inf, 
         upper=Inf, plots = TRUE, ...)

Arguments

Details

This function calls the wkde (weighted kernel density estimate) function.

Value

Just as for the integrate function, a list of class "integrate" with components

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics, 18, 505-526.

• Benaglia, T., Chauveau, D., Hunter, D. R., and Young, D. (2009), mixtools: An R package for analyzing finite mixture models. Journal of Statistical Software, 32(6):1-29.

See Also

npEM, wkde, integrate

Examples

# Mixture with mv gaussian model
set.seed(100)
m <- 2 # no. of components
r <- 3 # no. of repeated measures (coordinates)
lambda <- c(0.4, 0.6)
# Note:  Need first 2 coordinates conditionally iid due to block structure
mu <- matrix(c(0, 0, 0, 3, 3, 5), m, r, byrow=TRUE) # means 
sigma <- matrix(rep(1, 6), m, r, byrow=TRUE) # stdevs
blockid = c(1,1,2) # block structure of coordinates
n <- 200
x <- rmvnormmix(n, lambda, mu, sigma) # simulated data

# fit the model with "arbitrary" initial centers
centers <- matrix(c(0, 0, 0, 4, 4, 4), 2, 3, byrow=TRUE) 
a <- npEM(x, centers, blockid, eps=1e-8, verb=FALSE)

# Calculate integrated squared error for j=2, b=1:
j <- 2 # component
b <- 1 # block
coords <- a$blockid == b
ise.npEM(a, j, b, dnorm, lower=0, upper=10, plots=TRUE,
         mean=mu[j,coords][1], sd=sigma[j, coords][1])


# The following (lengthy) example recreates the normal multivariate 
# mixture model simulation from Benaglia et al (2009).  
mu <- matrix(c(0, 0, 0, 3, 4, 5), m, r, byrow=TRUE) 
nbrep <- 5  # Benaglia et al use 300 replications

# matrix for storing sums of Integrated Squared Errors 
ISE <- matrix(0,m,r,dimnames=list(Components=1:m, Blocks=1:r)) 

nblabsw <- 0 # no. of label switches
for (mc in 1:nbrep) {
  print(paste("REPETITION", mc))
	x <- rmvnormmix(n,lambda,mu,sigma) # simulated data
  a <- npEM(x, centers, verb=FALSE) #default:
	if (a$lambda[1] > a$lambda[2]) nblabsw <- nblabsw + 1
	for (j in 1:m) {  # for each component
		for (k in 1:r) { # for each coordinate; not assuming iid!
      # dnorm with correct mean, sd is the true density:
      ISE[j,k] <- ISE[j,k] + ise.npEM(a, j, k, dnorm, lower=mu[j,k]-5, 
               upper=mu[j,k]+5, plots=FALSE, mean=mu[j,k], 
               sd=sigma[j,k])$value
    }
  }
	MISE <- ISE/nbrep # Mean ISE
	sqMISE <- sqrt(MISE) # root-mean-integrated-squared error
}
sqMISE

Local Estimation for Lambda in Mixtures of Regressions

Description

Return local estimates of the mixing proportions from each component of a mixture of regressions model using output from an EM algorithm.

Usage

lambda(z, x, xi, h = NULL, kernel = c("Gaussian", "Beta", 
       "Triangle", "Cosinus", "Optcosinus"), g = 0)

Arguments

Value

lambda returns local estimates of the mixing proportions for the inputted x vector.

Note

lambda is for use within regmixEM.loc.

See Also

regmixEM.loc

Perturbation of Mixing Proportions

Description

Perturbs a set of mixing proportions by first scaling the mixing proportions, then taking the logit of the scaled values, perturbing them, and inverting back to produce a set of new mixing proportions.

Usage

lambda.pert(lambda, pert)

Arguments

Details

This function is called by regmixMH.

Value

lambda.pert returns new lambda values perturbed by pert.

See Also

regmixMH

Examples

set.seed(100)
x <- c(0.5, 0.2, 0.3)
lambda.pert(x, rcauchy(3))

Log-Density for Multinomial Distribution

Description

Return the logarithm of the multinomial density function.

Usage

ldmult(y, theta)

Arguments

Details

This function is called by multmixEM.

Value

ldmult returns the logarithm of the multinomial density with parameter theta, evaluated at y.

See Also

multmixEM

Examples

y <- c(2, 2, 10)
theta <- c(0.2, 0.3, 0.5)
ldmult(y, theta)

EM Algorithm for Mixtures of Logistic Regressions

Description

Returns EM algorithm output for mixtures of logistic regressions with arbitrarily many components.

Usage

logisregmixEM(y, x, N = NULL, lambda = NULL, beta = NULL, k = 2,
              addintercept = TRUE, epsilon = 1e-08, 
              maxit = 10000, verb = FALSE)

Arguments

Value

logisregmixEM returns a list of class mixEM with items:

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

See Also

poisregmixEM

Examples

## EM output for data generated from a 2-component logistic regression model.

set.seed(100)
beta <- matrix(c(1, .5, 2, -.8), 2, 2)
x <- runif(50, 0, 10)
x1 <- cbind(1, x)
xbeta <- x1%*%beta
N <- ceiling(runif(50, 50, 75))
w <- rbinom(50, 1, .3)
y <- w*rbinom(50, size = N, prob = (1/(1+exp(-xbeta[, 1]))))+
              (1-w)*rbinom(50, size = N, prob = 
              (1/(1+exp(-xbeta[, 2]))))
out.1 <- logisregmixEM(y, x, N, verb = TRUE, epsilon = 1e-01)
out.1

## EM output for data generated from a 2-component binary logistic regression model.

beta <- matrix(c(-10, .1, 20, -.1), 2, 2)
x <- runif(500, 50, 250)
x1 <- cbind(1, x)
xbeta <- x1%*%beta
w <- rbinom(500, 1, .3)
y <- w*rbinom(500, size = 1, prob = (1/(1+exp(-xbeta[, 1]))))+
              (1-w)*rbinom(500, size = 1, prob = 
              (1/(1+exp(-xbeta[, 2]))))
out.2 <- logisregmixEM(y, x, beta = beta, lambda = c(.3, .7), 
                       verb = TRUE, epsilon = 1e-01)
out.2

Produce Cutpoint Multinomial Data

Description

Change data into a matrix of multinomial counts using the cutpoint method and generate EM algorithm starting values for a k-component mixture of multinomials.

Usage

makemultdata(..., cuts)

Arguments

Details

The (i, j)th entry of the matrix y (for j < p) is equal to the number of entries in the ith column of x that are less than or equal to cuts[j]. The (i, p)th entry is equal to the number of entries greater than cuts[j].

Value

makemultdata returns an object which is a list with components:

References

Elmore, R. T., Hettmansperger, T. P. and Xuan, F. (2004) The Sign Statistic, One-Way Layouts and Mixture Models, Statistical Science 19(4), 579–587.

See Also

compCDF, multmixmodel.sel, multmixEM

Examples

## Randomly generated data.

set.seed(100)
y <- matrix(rpois(70, 6), 10, 7)
cuts <- c(2, 5, 7)
out1 <- makemultdata(y, cuts = cuts)
out1

## The sulfur content of the coal seams in Texas.

A <- c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17)
B <- c(1.69, 0.64, .9, 1.41, 1.01, .84, 1.28, 1.59)
C <- c(1.56, 1.22, 1.32, 1.39, 1.33, 1.54, 1.04, 2.25, 1.49)
D <- c(1.3, .75, 1.26, .69, .62, .9, 1.2, .32)
E <- c(.73, .8, .9, 1.24, .82, .72, .57, 1.18, .54, 1.3)

out2 <- makemultdata(A, B, C, D, E, 
                     cuts = median(c(A, B, C, D, E)))
out2

## The reaction time data.

data(RTdata)
out3 <- makemultdata(RTdata, cuts = 
                     100*c(5, 10, 12, 14, 16, 20, 25, 30, 40, 50))
dim(out3$y)
out3$y[1:10,]

Calculates the Square Root of a Diagonalizable Matrix

Description

Returns the square root of a diagonalizable matrix.

Usage

matsqrt(x)

Arguments

Details

This function is called by regcr.

Value

matsqrt returns the square root of x.

See Also

regcr

Examples

a <- matrix(c(1, -0.2, -0.2, 1), 2, 2)
matsqrt(a)

Initializations for Various EM Algorithms in 'mixtools'

Description

Internal intialization functions for EM algorithms in the package mixtools.

Usage

flaremix.init(y, x, lambda = NULL, beta = NULL, sigma = NULL,
              alpha = NULL)
gammamix.init(x, lambda = NULL, alpha = NULL, beta = NULL, 
              k = 2)
logisregmix.init(y, x, N, lambda = NULL, beta = NULL, k = 2)
multmix.init(y, lambda = NULL, theta = NULL, k = 2)
mvnormalmix.init(x, lambda = NULL, mu = NULL, sigma = NULL, 
                 k = 2, arbmean = TRUE, arbvar = TRUE)
normalmix.init(x, lambda = NULL, mu = NULL, s = NULL, k = 2, 
               arbmean = TRUE, arbvar = TRUE)
poisregmix.init(y, x, lambda = NULL, beta = NULL, k = 2)
regmix.init(y, x, lambda = NULL, beta = NULL, s = NULL, k = 2, 
            addintercept = TRUE, arbmean = TRUE, arbvar=TRUE)
regmix.lambda.init(y, x, lambda = NULL, beta = NULL, s = NULL,
                   k = 2, addintercept = TRUE, arbmean = TRUE,
                   arbvar = TRUE)
regmix.mixed.init(y, x, w = NULL, sigma = NULL, 
                  arb.sigma = TRUE, alpha = NULL, lambda = NULL, 
                  mu = NULL, R = NULL, arb.R = TRUE, k = 2, 
                  mixed = FALSE, addintercept.fixed = FALSE,
                  addintercept.random = TRUE)
repnormmix.init(x, lambda = NULL, mu = NULL, s = NULL, k = 2, 
                arbmean = TRUE, arbvar = TRUE)
segregmix.init(y, x, lambda = NULL, beta = NULL, s = NULL, k = 2,
               seg.Z, psi, psi.locs = NULL)

Details

These are usually not to be called by the user. Definitions of the arguments appear in the respective EM algorithms.

Internal 'mixtools' Functions

Description

Internal kernel, semiparametric-related, and miscellaneous functions for the package mixtools.

Usage

dexpmixt(t, lam, rate)
HRkde(cpd, u = cpd[,1], kernelft = triang_wkde, 
      bw = rep(bw.nrd0(as.vector(cpd[,1])), length(cpd[,1])))
inv.logit(eta)
kern.B(x, xi, h, g = 0)
kern.C(x, xi, h)
kern.G(x, xi, h)
kern.O(x, xi, h)
kern.T(x, xi, h)
kfoldCV(h, x, nbsets = 2, w = rep(1, length(x)), 
        lower = mean(x) - 5*sd(x), upper = mean(x) + 5*sd(x))
KMintegrate(s) 
KMod(cpd, already.ordered = TRUE)
ldc(data, class, score)
logit(mu)
npMSL_old(x, mu0, blockid = 1:ncol(x),
          bw=bw.nrd0(as.vector(as.matrix(x))), samebw = TRUE,
          h=bw, eps=1e-8, maxiter=500, bwiter = maxiter,
          ngrid = 200, post = NULL, verb = TRUE)
plotseq(x, ...)
rlnormscalemix(n, lambda=1, meanlog=1, sdlog=1, scale=0.1)
splitsample(n, nbsets = 2)
triang_wkde(t, u=t, w=rep(1/length(t),length(t)), bw=rep(bw.nrd0(t), length(t)))
wbw.kCV(x, nbfold = 5, w = rep(1, length(x)), 
        hmin = 0.1*hmax, hmax = NULL)

Arguments

Details

These are usually not to be called by the user.

See Also

npMSL

Mixturegrams

Description

Construct a mixturegram for determining an apporpriate number of components.

Usage

mixturegram(data, pmbs, method = c("pca", "kpca", "lda"), all.n = FALSE,
			id.con = NULL, score = 1, iter.max = 50, nstart = 25, ...)

Arguments

Value

mixturegram returns a mixturegram where the profiles are plotted over component values of k = 1,...,K.

References

Young, D. S., Ke, C., and Zeng, X. (2018) The Mixturegram: A Visualization Tool for Assessing the Number of Components in Finite Mixture Models, Journal of Computational and Graphical Statistics, 27(3), 564–575.

See Also

boot.comp

Examples

##Data generated from a 2-component mixture of normals.


set.seed(100)
n <- 100
w <- rmultinom(n,1,c(.3,.7))
y <- sapply(1:n,function(i) w[1,i]*rnorm(1,-6,1) +
            w[2,i]*rnorm(1,0,1))

selection <- function(i,data,rep=30){
  out <- replicate(rep,normalmixEM(data,epsilon=1e-06,
  				   k=i,maxit=5000),simplify=FALSE)
  counts <- lapply(1:rep,function(j) 
                   table(apply(out[[j]]$posterior,1,
                   which.max)))
  counts.length <- sapply(counts, length)
  counts.min <- sapply(counts, min)
  counts.test <- (counts.length != i)|(counts.min < 5)
  if(sum(counts.test) > 0 & sum(counts.test) < rep) 
  	out <- out[!counts.test]
  l <- unlist(lapply(out, function(x) x$loglik))
  tmp <- out[[which.max(l)]]
}

all.out <- lapply(2:5, selection, data = y, rep = 2)
pmbs <- lapply(1:length(all.out), function(i) 
			   all.out[[i]]$post)
mixturegram(y, pmbs, method = "pca", all.n = FALSE,
		    id.con = NULL, score = 1, 
		    main = "Mixturegram (Well-Separated Data)")

EM Algorithm for Mixtures of Multinomials

Description

Return EM algorithm output for mixtures of multinomial distributions.

Usage

multmixEM(y, lambda = NULL, theta = NULL, k = 2,
          maxit = 10000, epsilon = 1e-08, verb = FALSE)

Arguments

Value

multmixEM returns a list of class mixEM with items:

References

• McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

• Elmore, R. T., Hettmansperger, T. P. and Xuan, F. (2004) The Sign Statistic, One-Way Layouts and Mixture Models, Statistical Science 19(4), 579–587.

See Also

compCDF, makemultdata, multmixmodel.sel

Examples

## The sulfur content of the coal seams in Texas

set.seed(100)
A <- c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17)
B <- c(1.69, 0.64, .9, 1.41, 1.01, .84, 1.28, 1.59) 
C <- c(1.56, 1.22, 1.32, 1.39, 1.33, 1.54, 1.04, 2.25, 1.49) 
D <- c(1.3, .75, 1.26, .69, .62, .9, 1.2, .32) 
E <- c(.73, .8, .9, 1.24, .82, .72, .57, 1.18, .54, 1.3)

dis.coal <- makemultdata(A, B, C, D, E, 
                         cuts = median(c(A, B, C, D, E)))
em.out <- multmixEM(dis.coal)
em.out[1:4]

Model Selection Mixtures of Multinomials

Description

Assess the number of components in a mixture of multinomials model using the Akaike's information criterion (AIC), Schwartz's Bayesian information criterion (BIC), Bozdogan's consistent AIC (CAIC), and Integrated Completed Likelihood (ICL).

Usage

multmixmodel.sel(y, comps = NULL, ...)

Arguments

Value

multmixmodel.sel returns a table summarizing the AIC, BIC, CAIC, ICL, and log-likelihood values along with the winner (the number with the lowest aforementioned values).

See Also

compCDF, makemultdata, multmixEM

Examples

##Data generated using the multinomial cutpoint method.

set.seed(100)
x <- matrix(rpois(70, 6), 10, 7) 
x.new <- makemultdata(x, cuts = 5)
multmixmodel.sel(x.new$y, comps = c(1,2), epsilon = 1e-03)

EM Algorithm for Mixtures of Multivariate Normals

Description

Return EM algorithm output for mixtures of multivariate normal distributions.

Usage

mvnormalmixEM(x, lambda = NULL, mu = NULL, sigma = NULL, k = 2,
              arbmean = TRUE, arbvar = TRUE, epsilon = 1e-08, 
              maxit = 10000, verb = FALSE)

Arguments

Value

normalmixEM returns a list of class mixEM with items:

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

See Also

normalmixEM

Examples

##Fitting randomly generated data with a 2-component location mixture of bivariate normals.

set.seed(100)
x.1 <- rmvnorm(40, c(0, 0))
x.2 <- rmvnorm(60, c(3, 4))
X.1 <- rbind(x.1, x.2)
mu <- list(c(0, 0), c(3, 4))

out.1 <- mvnormalmixEM(X.1, arbvar = FALSE, mu = mu,
                       epsilon = 1e-02)
out.1[2:5]

##Fitting randomly generated data with a 2-component scale mixture of bivariate normals.

x.3 <- rmvnorm(40, c(0, 0), sigma = 
               matrix(c(200, 1, 1, 150), 2, 2))
x.4 <- rmvnorm(60, c(0, 0))
X.2 <- rbind(x.3, x.4)
lambda <- c(0.40, 0.60)
sigma <- list(diag(1, 2), matrix(c(200, 1, 1, 150), 2, 2))
 
out.2 <- mvnormalmixEM(X.2, arbmean = FALSE,
                       sigma = sigma, lambda = lambda,
                       epsilon = 1e-02)
out.2[2:5]

EM-like Algorithm for Nonparametric Mixture Models with Conditionally Independent Multivariate Component Densities

Description

An extension of the original npEM algorithm, for mixtures of multivariate data where the coordinates of a row (case) in the data matrix are assumed to be made of independent but multivariate blocks (instead of just coordinates), conditional on the mixture component (subpopulation) from which they are drawn (Chauveau and Hoang 2015).

Usage

mvnpEM(x, mu0, blockid = 1:ncol(x), samebw = TRUE, 
       bwdefault = apply(x,2,bw.nrd0), init = NULL,
       eps = 1e-8, maxiter = 500, verb = TRUE)

Arguments

Value

mvnpEM returns a list of class mvnpEM with the following items:

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics, 18, 505-526.

• Benaglia, T., Chauveau, D. and Hunter, D.R. (2011), Bandwidth Selection in an EM-like algorithm for nonparametric multivariate mixtures. Nonparametric Statistics and Mixture Models: A Festschrift in Honor of Thomas P. Hettmansperger. World Scientific Publishing Co., pages 15-27.

• Chauveau, D., and Hoang, V. T. L. (2015), Nonparametric mixture models with conditionally independent multivariate component densities, Preprint under revision. https://hal.science/hal-01094837

See Also

plot.mvnpEM, npEM

Examples

# Example as in Chauveau and Hoang (2015) with 6 coordinates
## Not run: 
m=2; r=6; blockid <-c(1,1,2,2,3,3) # 3 bivariate blocks 
# generate some data x ...
a <- mvnpEM(x, mu0=2, blockid, samebw=F) # adaptive bandwidth
plot(a) # this S3 method produces 6 plots of univariate marginals
summary(a)
## End(Not run)

EM Algorithm for Mixtures of Univariate Normals

Description

Return EM algorithm output for mixtures of normal distributions.

Usage

normalmixEM(x, lambda = NULL, mu = NULL, sigma = NULL, k = 2, 
            mean.constr = NULL, sd.constr = NULL,
            epsilon = 1e-08, maxit = 1000, maxrestarts = 20, 
            verb = FALSE, fast = FALSE, ECM = FALSE,
            arbmean = TRUE, arbvar = TRUE) 

Arguments

Details

This is the standard EM algorithm for normal mixtures that maximizes the conditional expected complete-data log-likelihood at each M-step of the algorithm. If desired, the EM algorithm may be replaced by an ECM algorithm (see ECM argument) that alternates between maximizing with respect to the mu and lambda while holding sigma fixed, and maximizing with respect to sigma and lambda while holding mu fixed. In the case where arbmean is FALSE and arbvar is TRUE, there is no closed-form EM algorithm, so the ECM option is forced in this case.

Value

normalmixEM returns a list of class mixEM with items:

References

• McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

• Meng, X.-L. and Rubin, D. B. (1993) Maximum Likelihood Estimation Via the ECM Algorithm: A General Framework, Biometrika 80(2): 267-278.

• Benaglia, T., Chauveau, D., Hunter, D. R., and Young, D. mixtools: An R package for analyzing finite mixture models. Journal of Statistical Software, 32(6):1-29, 2009.

See Also

mvnormalmixEM, normalmixEM2comp, normalmixMMlc, spEMsymloc

Examples

##Analyzing the Old Faithful geyser data with a 2-component mixture of normals.

data(faithful)
attach(faithful)
set.seed(100)
system.time(out<-normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03))
out
system.time(out2<-normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03, fast=TRUE))
out2 # same thing but much faster

Fast EM Algorithm for 2-Component Mixtures of Univariate Normals

Description

Return EM algorithm output for mixtures of univariate normal distributions for the special case of 2 components, exploiting the simple structure of the problem to speed up the code.

Usage

normalmixEM2comp(x, lambda, mu, sigsqrd, eps= 1e-8, maxit = 1000, verb=FALSE)

Arguments

Details

This code is written to be very fast, sometimes more than an order of magnitude faster than normalmixEM for the same problem. It is less numerically stable that normalmixEM in the sense that it does not safeguard against underflow as carefully.

Note that when the two components are assumed to have unequal variances, the loglikelihood is unbounded. However, in practice this is rarely a problem and quite often the algorithm converges to a "nice" local maximum.

Value

normalmixEM2comp returns a list of class mixEM with items:

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

See Also

mvnormalmixEM, normalmixEM

Examples

##Analyzing the Old Faithful geyser data with a 2-component mixture of normals.

data(faithful)
attach(faithful)
set.seed(100)
system.time(out <- normalmixEM2comp(waiting, lambda=.5, 
            mu=c(50,80), sigsqrd=100))
out$all.loglik # Note:  must be monotone increasing

# Compare elapsed time with more general version
system.time(out2 <- normalmixEM(waiting, lambda=c(.5,.5), 
            mu=c(50,80), sigma=c(10,10), arbvar=FALSE))
out2$all.loglik # Values should be identical to above

EC-MM Algorithm for Mixtures of Univariate Normals with linear constraints

Description

Return EC-MM (see below) algorithm output for mixtures of normal distributions with linear constraints on the means and variances parameters, as in Chauveau and Hunter (2013). The linear constraint for the means is of the form \mu = M \beta + C, where M and C are matrix and vector specified as parameters. The linear constraints for the variances are actually specified on the inverse variances, by \pi = A \gamma, where \pi is the vector of inverse variances, and A is a matrix specified as a parameter (see below).

Usage

normalmixMMlc(x, lambda = NULL, mu = NULL, sigma = NULL, k = 2,
              mean.constr = NULL, mean.lincstr = NULL, 
              mean.constant = NULL, var.lincstr = NULL, 
              gparam = NULL, epsilon = 1e-08, maxit = 1000, 
              maxrestarts=20, verb = FALSE) 

Arguments

Details

This is a specific "EC-MM" algorithm for normal mixtures with linear constraints on the means and variances parameters. EC-MM here means that this algorithm is similar to an ECM algorithm as in Meng and Rubin (1993), except that it uses conditional MM (Minorization-Maximization)-steps instead of simple M-steps. Conditional means that it alternates between maximizing with respect to the mu and lambda while holding sigma fixed, and maximizing with respect to sigma and lambda while holding mu fixed. This ECM generalization of EM is forced in the case of linear constraints because there is no closed-form EM algorithm.

Value

normalmixMMlc returns a list of class mixEM with items:

Author(s)

Didier Chauveau

References

• McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley & Sons, Inc.

• Meng, X.-L. and Rubin, D. B. (1993) Maximum Likelihood Estimation Via the ECM Algorithm: A General Framework, Biometrika 80(2): 267-278.

• Chauveau, D. and Hunter, D.R. (2013) ECM and MM algorithms for mixtures with constrained parameters, preprint https://hal.science/hal-00625285.

• Thomas, H., Lohaus, A., and Domsch, H. (2011) Stable Unstable Reliability Theory, British Journal of Mathematical and Statistical Psychology 65(2): 201-221.

See Also

normalmixEM, mvnormalmixEM, normalmixEM2comp, tauequivnormalmixEM

Examples

## Analyzing synthetic data as in the tau equivalent model  
## From Thomas et al (2011), see also Chauveau and Hunter (2013)
## a 3-component mixture of normals with linear constraints.
lbd <- c(0.6,0.3,0.1); m <- length(lbd)
sigma <- sig0 <- sqrt(c(1,9,9))
# means constaints mu = M beta
M <- matrix(c(1,1,1,0,-1,1), 3, 2)
beta <- c(1,5) # unknown constrained mean
mu0 <- mu <- as.vector(M %*% beta)
# linear constraint on the inverse variances pi = A.g
A <- matrix(c(1,1,1,0,1,0), m, 2, byrow=TRUE)
iv0 <- 1/(sig0^2)
g0 <- c(iv0[2],iv0[1] - iv0[2]) # gamma^0 init 

# simulation and EM fits
set.seed(50); n=100; x <- rnormmix(n,lbd,mu,sigma)
s <- normalmixEM(x,mu=mu0,sigma=sig0,maxit=2000) # plain EM
# EM with var and mean linear constraints
sc <- normalmixMMlc(x, lambda=lbd, mu=mu0, sigma=sig0,
					mean.lincstr=M, var.lincstr=A, gparam=g0)
# plot and compare both estimates
dnormmixt <- function(t, lam, mu, sig){
	m <- length(lam); f <- 0
	for (j in 1:m) f <- f + lam[j]*dnorm(t,mean=mu[j],sd=sig[j])
	f}
t <- seq(min(x)-2, max(x)+2, len=200)
hist(x, freq=FALSE, col="lightgrey", 
		ylim=c(0,0.3), ylab="density",main="")
lines(t, dnormmixt(t, lbd, mu, sigma), col="darkgrey", lwd=2) # true
lines(t, dnormmixt(t, s$lambda, s$mu, s$sigma), lty=2) 
lines(t, dnormmixt(t, sc$lambda, sc$mu, sc$sigma), col=1, lty=3)
legend("topleft", c("true","plain EM","constr EM"), 
	col=c("darkgrey",1,1), lty=c(1,2,3), lwd=c(2,1,1))

Nonparametric EM-like Algorithm for Mixtures of Independent Repeated Measurements

Description

Returns nonparametric EM algorithm output (Benaglia et al, 2009) for mixtures of multivariate (repeated measures) data where the coordinates of a row (case) in the data matrix are assumed to be independent, conditional on the mixture component (subpopulation) from which they are drawn.

Usage

npEM(x, mu0, blockid = 1:ncol(x), 
     bw = bw.nrd0(as.vector(as.matrix(x))), samebw = TRUE, 
     h = bw, eps = 1e-8, 
     maxiter = 500, stochastic = FALSE, verb = TRUE)

Arguments

Value

npEM returns a list of class npEM with the following items:

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics, 18, 505-526.

• Benaglia, T., Chauveau, D., Hunter, D. R., and Young, D. (2009), mixtools: An R package for analyzing finite mixture models. Journal of Statistical Software, 32(6):1-29.

• Benaglia, T., Chauveau, D. and Hunter, D.R. (2011), Bandwidth Selection in an EM-like algorithm for nonparametric multivariate mixtures. Nonparametric Statistics and Mixture Models: A Festschrift in Honor of Thomas P. Hettmansperger. World Scientific Publishing Co., pages 15-27.

• Bordes, L., Chauveau, D., and Vandekerkhove, P. (2007), An EM algorithm for a semiparametric mixture model, Computational Statistics and Data Analysis, 51: 5429-5443.

See Also

plot.npEM, normmixrm.sim, spEMsymloc, spEM, plotseq.npEM

Examples

## Examine and plot water-level task data set.

## First, try a 3-component solution where no two coordinates are
## assumed i.d.
data(Waterdata)
set.seed(100)
## Not run: 
a <- npEM(Waterdata[,3:10], mu0=3, bw=4) # Assume indep but not iid
plot(a) # This produces 8 plots, one for each coordinate

## End(Not run)

## Next, same thing but pairing clock angles that are directly opposite one
## another (1:00 with 7:00, 2:00 with 8:00, etc.)
## Not run: 
b <- npEM(Waterdata[,3:10], mu0=3, blockid=c(4,3,2,1,3,4,1,2), bw=4) # iid in pairs
plot(b) # Now only 4 plots, one for each block

## End(Not run)

Nonparametric EM-like Algorithm for Mixtures of Independent Repeated Measurements - Maximum Smoothed Likelihood version

Description

Returns nonparametric Smoothed Likelihood algorithm output (Levine et al, 2011) for mixtures of multivariate (repeated measures) data where the coordinates of a row (case) in the data matrix are assumed to be independent, conditional on the mixture component (subpopulation) from which they are drawn.

Usage

npMSL(x, mu0, blockid = 1:ncol(x), 
      bw = bw.nrd0(as.vector(as.matrix(x))), samebw = TRUE, 
      bwmethod = "S", h = bw, eps = 1e-8, 
      maxiter=500, bwiter = maxiter, nbfold = NULL,
      ngrid=200, post=NULL, verb = TRUE)

Arguments

Value

npMSL returns a list of class npEM with the following items:

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics, 18, 505-526.

• Benaglia, T., Chauveau, D. and Hunter, D.R. (2011), Bandwidth Selection in an EM-like algorithm for nonparametric multivariate mixtures. Nonparametric Statistics and Mixture Models: A Festschrift in Honor of Thomas P. Hettmansperger. World Scientific Publishing Co., pages 15-27.

• Chauveau D., Hunter D. R. and Levine M. (2014), Semi-Parametric Estimation for Conditional Independence Multivariate Finite Mixture Models. Preprint (under revision).

• Levine, M., Hunter, D. and Chauveau, D. (2011), Maximum Smoothed Likelihood for Multivariate Mixtures, Biometrika 98(2): 403-416.

See Also

npEM, plot.npEM, normmixrm.sim, spEMsymloc, spEM, plotseq.npEM

Examples

## Examine and plot water-level task data set.
## Block structure pairing clock angles that are directly opposite one
## another (1:00 with 7:00, 2:00 with 8:00, etc.)
set.seed(111) # Ensure that results are exactly reproducible
data(Waterdata)
blockid <- c(4,3,2,1,3,4,1,2) # see Benaglia et al (2009a)

## Not run: 
a <- npEM(Waterdata[,3:10], mu0=3, blockid=blockid, bw=4)  # npEM solution
b <- npMSL(Waterdata[,3:10], mu0=3, blockid=blockid, bw=4) # smoothed version

# Comparisons on the 4 default plots, one for each block
par(mfrow=c(2,2))
for (l in 1:4){
plot(a, blocks=l, breaks=5*(0:37)-92.5,
	xlim=c(-90,90), xaxt="n",ylim=c(0,.035), xlab="")
plot(b, blocks=l, hist=FALSE, newplot=FALSE, addlegend=FALSE, lty=2,
	dens.col=1)
axis(1, at=30*(1:7)-120, cex.axis=1)
legend("topleft",c("npMSL"),lty=2, lwd=2)}

## End(Not run)

Constraint Function

Description

Constraint function for some mixture EM algorithms.

Usage

parse.constraints(constr, k = 2, allsame = FALSE)

Arguments

Details

This function is not intended to be called by the user.

Permutation Function

Description

Enumerates the possible permutations of a specified size from the elements of a vector having the same size.

Usage

perm(n, r, v = 1:n)

Arguments

Details

This function is called by segregmixEM and the associated internal functions. This is also a simplified version of the function permutations found in the package gtools.

Value

perm returns a matrix where each row contains one of the permutations of length r.

See Also

segregmixEM

Examples

perm(3, 3, 2:4)

Various Plots Pertaining to Mixture Models

Description

Takes an object of class mixEM and returns various graphical output for select mixture models.

Usage

## S3 method for class 'mixEM'
plot(x, whichplots = 1, 
     loglik = 1 %in% whichplots,
     density = 2 %in% whichplots,
     xlab1="Iteration", ylab1="Log-Likelihood",
     main1="Observed Data Log-Likelihood", col1=1, lwd1=2,
     xlab2=NULL, ylab2=NULL, main2=NULL, col2=NULL, 
     lwd2=2, alpha = 0.05, marginal = FALSE, ...)

Arguments

Value

plot.mixEM returns a plot of the log-likelihood versus the EM iterations by default for all objects of class mixEM. In addition, other plots may be produced for the following k-component mixture model functions:

See Also

post.beta

Examples

 
##Analyzing the Old Faithful geyser data with a 2-component mixture of normals.

data(faithful)
attach(faithful)
set.seed(100)
out <- normalmixEM(waiting, arbvar = FALSE, verb = TRUE,
                   epsilon = 1e-04)
plot(out, density = TRUE, w = 1.1)

##Fitting randomly generated data with a 2-component location mixture of bivariate normals.

x.1 <- rmvnorm(40, c(0, 0))
x.2 <- rmvnorm(60, c(3, 4))
X.1 <- rbind(x.1, x.2)

out.1 <- mvnormalmixEM(X.1, arbvar = FALSE, verb = TRUE,
                       epsilon = 1e-03)
plot(out.1, density = TRUE, alpha = c(0.01, 0.05, 0.10), 
     marginal = TRUE)

Various Plots Pertaining to Mixture Model Output Using MCMC Methods

Description

Takes an object of class mixMCMC and returns various graphical output for select mixture models.

Usage

 
## S3 method for class 'mixMCMC'
plot(x, trace.plots = TRUE, 
     summary.plots = FALSE, burnin = 2000, ...) 

Arguments

Value

plot.mixMCMC returns trace plots of the various parameters estimated by the MCMC methods for all objects of class mixMCMC. In addition, other plots may be produced for the following k-component mixture model functions:

See Also

regcr

Examples

 
## M-H algorithm for NOdata with acceptance rate about 40%.

data(NOdata)
attach(NOdata)
set.seed(100)
beta <- matrix(c(1.3, -0.1, 0.6, 0.1), 2, 2)
sigma <- c(.02, .05)
MH.out <- regmixMH(Equivalence, NO, beta = beta, s = sigma, 
                   sampsize = 2500, omega = .0013)
plot(MH.out, summary.plots = TRUE, burnin = 2450, 
     alpha = 0.01)

Plots of Marginal Density Estimates from the mvnpEM Algorithm Output

Description

Takes an object of class mvnpEM, as the one returned by the mvnpEM algorithm, and returns a set of plots of the density estimates for each coordinate within each multivariate block. All the components are displayed on each plot so it is possible to see the mixture structure for each coordinate and block. The final bandwidth values are also displayed, in a format depending on the bandwidth strategy .

Usage

## S3 method for class 'mvnpEM'
plot(x, truenorm = FALSE, lambda = NULL, mu = NULL, v = NULL, 
            lgdcex = 1, ...)

Arguments

Value

plot.mvnpEM currently just plots the figure.

See Also

mvnpEM, npEM, density.npEM

Examples

 
# example as in Chauveau and Hoang (2015) with 6 coordinates
## Not run: 
m=2; r=6; blockid <-c(1,1,2,2,3,3) # 3 bivariate blocks 
# generate some data x ...
a <- mvnpEM(x, mu0=2, blockid, samebw=F) # adaptive bandwidth
plot(a) # this S3 method produces 6 plots of univariate marginals
summary(a)
## End(Not run)

Plot Nonparametric or Semiparametric EM Output

Description

Takes an object of class npEM and returns a set of plots of the density estimates for each block and each component. There is one plot per block, with all the components displayed on each block so it is possible to see the mixture structure for each block.

Usage

## S3 method for class 'npEM'
plot(x, blocks = NULL, hist=TRUE, addlegend = TRUE,
      scale=TRUE, title=NULL, breaks="Sturges", ylim=NULL, dens.col,
      newplot = TRUE, pos.legend = "topright", cex.legend = 1, ...)         
## S3 method for class 'spEM'
plot(x, blocks = NULL, hist=TRUE, addlegend = TRUE,
      scale=TRUE, title=NULL, breaks="Sturges", ylim=NULL, dens.col,
      newplot = TRUE, pos.legend = "topright", cex.legend = 1, ...)

Arguments

Value

plot.npEM returns a list with two elements:

See Also

npEM, density.npEM, spEMsymloc, plotseq.npEM

Examples

 
## Examine and plot water-level task data set.

## First, try a 3-component solution where no two coordinates are
## assumed i.d.
data(Waterdata)
set.seed(100)
## Not run: 
a <- npEM(Waterdata[,3:10], 3, bw=4)
par(mfrow=c(2,4))
plot(a) # This produces 8 plots, one for each coordinate

## End(Not run)

## Not run: 
## Next, same thing but pairing clock angles that are directly opposite one
## another (1:00 with 7:00, 2:00 with 8:00, etc.)
b <- npEM(Waterdata[,3:10], 3, blockid=c(4,3,2,1,3,4,1,2), bw=4)
par(mfrow=c(2,2))
plot(b) # Now only 4 plots, one for each block

## End(Not run)

Plot mixture pdf for the semiparametric mixture model output by spEMsymlocN01

Description

Plot mixture density for the semiparametric mixture model output by spEMsymlocN01, with one component known and set to normal(0,1), and a symmetric nonparametric density with location parameter.

Usage

## S3 method for class 'spEMN01'
plot(x, bw = x$bandwidth, knownpdf = dnorm, add.plot = FALSE, ...)

Arguments

Value

A plot of the density of the mixture

Author(s)

Didier Chauveau

References

• Chauveau, D., Saby, N., Orton, T. G., Lemercier B., Walter, C. and Arrouys, D. Large-scale simultaneous hypothesis testing in soil monitoring: A semi-parametric mixture approach, preprint (2013).

See Also

spEMsymlocN01

Plot False Discovery Rate (FDR) estimates from output by EM-like strategies

Description

Plot FDR(p_i) estimates against index of sorted p-values from, e.g., normalmixEM or the semiparametric mixture model posterior probabilities output by spEMsymlocN01, or any EM-algorithm like normalmixEM which returns posterior probabilities. The function can simultaneously plot FDR estimates from two strategies for comparison. Plot of the true FDR can be added if complete data are available (typically in simulation studies).

Usage

plotFDR(post1, post2 = NULL, lg1 = "FDR 1", lg2 = NULL, title = NULL, 
        compH0 = 1, alpha = 0.1, complete.data = NULL, pctfdr = 0.3)

Arguments

Value

A plot of one or two FDR estimates, with the true FDR if available

Author(s)

Didier Chauveau

References

• Chauveau, D., Saby, N., Orton, T. G., Lemercier B., Walter, C. and Arrouys, D. Large-scale simultaneous hypothesis testing in monitoring carbon content from French soil database – A semi-parametric mixture approach, Geoderma 219-220 (2014), 117-124.

See Also

spEMsymlocN01

Plot sequences from the EM algorithm for censored mixture of exponentials

Description

Function for plotting sequences of estimates along iterations, from an object returned by the expRMM_EM, an EM algorithm for mixture of exponential distributions with randomly right censored data (see reference below).

Usage

  plotexpRMM(a, title=NULL, rowstyle=TRUE, subtitle=NULL, ...)

Arguments

Value

The plot returned

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: expRMM_EM, summary.mixEM, plot.mixEM.

Other models and algorithms for censored lifetime data (name convention is model_algorithm): weibullRMM_SEM, spRMM_SEM.

Examples

n=300 # sample size
m=2   # number of mixture components
lambda <- c(1/3,1-1/3); rate <- c(1,1/10) # mixture parameters
set.seed(1234)
x <- rexpmix(n, lambda, rate) # iid ~ exponential mixture
cs=runif(n,0,max(x)) # Censoring (uniform) and incomplete data
t <- apply(cbind(x,cs),1,min) # observed or censored data
d <- 1*(x <= cs)              # censoring indicator

###### EM for RMM, exponential lifetimes
l0 <- rep(1/m,m); r0 <- c(1, 0.5) # "arbitrary" initial values
a <- expRMM_EM(t, d, lambda=l0, rate=r0, k = m)
summary(a)             # EM estimates etc
plotexpRMM(a, lwd=2) # plot of EM sequences

Plot False Discovery Rate (FDR) estimates from output by EM-like strategies using plotly

Description

This is an updated version of plotFDR. For more technical details, please refer to plotFDR.

Usage

plotly_FDR(post1, post2=NULL, lg1="FDR 1", lg2=NULL, 
          compH0=1, alpha=0.1, complete.data =NULL, pctfdr=0.3,
          col = NULL, width = 3 ,
          title = NULL , title.size = 15 , title.x = 0.5 , title.y = 0.95,
          xlab = "Index" , xlab.size = 15 , xtick.size = 15,
          ylab = "Probability" , ylab.size = 15 , ytick.size = 15,
          legend.text = "" , legend.text.size = 15 , legend.size = 15)

Arguments

Value

A plot of one or two FDR estimates, with the true FDR if available

Author(s)

Didier Chauveau

References

• Chauveau, D., Saby, N., Orton, T. G., Lemercier B., Walter, C. and Arrouys, D. Large-scale simultaneous hypothesis testing in monitoring carbon content from French soil database – A semi-parametric mixture approach, Geoderma 219-220 (2014), 117-124.

See Also

spEMsymlocN01, plotFDR

Examples

## Probit transform of p-values
## from a Beta-Uniform mixture model
## comparion of parametric and semiparametric EM fit
## Note: in actual situations n=thousands
set.seed(50)
n=300 # nb of multiple tests
m=2 # 2 mixture components
a=c(1,0.1); b=c(1,1); lambda=c(0.6,0.4) # parameters
z=sample(1:m, n, rep=TRUE, prob = lambda)
p <- rbeta(n, shape1 = a[z], shape2 = b[z]) # p-values
o <- order(p)
cpd <- cbind(z,p)[o,] # sorted complete data, z=1 if H0, 2 if H1
p <- cpd[,2] # sorted p-values
y <- qnorm(p) # probit transform of the pvalues
# gaussian EM fit with component 1 constrained to N(0,1)
s1 <- normalmixEM(y, mu=c(0,-4),
                  mean.constr = c(0,NA), sd.constr = c(1,NA))
s2 <- spEMsymlocN01(y, mu0 = c(0,-3)) # spEM with N(0,1) fit
plotly_FDR(s1$post, s2$post, lg1 = "normalmixEM", lg2 = "spEMsymlocN01",
           complete.data = cpd) # with true FDR computed from z

Plot the Component CDF using plotly

Description

Plot the components' CDF via the posterior probabilities using plotly.

Usage

plotly_compCDF(data, weights, x=seq(min(data, na.rm=TRUE), max(data, na.rm=TRUE), 
               len=250), comp=1:NCOL(weights), makeplot=TRUE,
               cex = 3, width = 3,
               legend.text = "Composition", legend.text.size = 15, legend.size = 15,
               title = "Empirical CDF", title.x = 0.5, title.y = 0.95, title.size = 15,
               xlab = "Data", xlab.size = 15, xtick.size = 15,
               ylab = "Probability", ylab.size = 15, ytick.size = 15,
               col.comp = NULL) 

Arguments

Details

When makeplot is TRUE, a line plot is produced of the CDFs evaluated at x. The plot is not a step function plot; the points (x, CDF(x)) are simply joined by line segments.

Value

A matrix with length(comp) rows and length(x) columns in which each row gives the CDF evaluated at each point of x.

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

Elmore, R. T., Hettmansperger, T. P. and Xuan, F. (2004) The Sign Statistic, One-Way Layouts and Mixture Models, Statistical Science 19(4), 579–587.

See Also

makemultdata, multmixmodel.sel, multmixEM, compCDF.

Examples

## The sulfur content of the coal seams in Texas
set.seed(100)
A <- c(1.51, 1.92, 1.08, 2.04, 2.14, 1.76, 1.17)
B <- c(1.69, 0.64, .9, 1.41, 1.01, .84, 1.28, 1.59)
C <- c(1.56, 1.22, 1.32, 1.39, 1.33, 1.54, 1.04, 2.25, 1.49)
D <- c(1.3, .75, 1.26, .69, .62, .9, 1.2, .32)
E <- c(.73, .8, .9, 1.24, .82, .72, .57, 1.18, .54, 1.3)
dis.coal <- makemultdata(A, B, C, D, E,
                         cuts = median(c(A, B, C, D, E)))
temp <- multmixEM(dis.coal)
## Now plot the components' CDF via the posterior probabilities
plotly_compCDF(dis.coal$x, temp$posterior, xlab="Sulfur")

Draw Two-Dimensional Ellipse Based on Mean and Covariance using plotly

Description

This is an updated version of ellipse. For more technical details, please refer to ellipse.

Usage

plotly_ellipse(mu, sigma, alpha=.05, npoints=250,
               draw=TRUE, cex = 3, col = "#1f77b4", lwd = 3,
               title = "", title.x = 0.5, title.y = 0.95, title.size = 15,
               xlab = "X", xlab.size = 15, xtick.size = 15,
               ylab = "Y", ylab.size = 15, ytick.size = 15)

Arguments

Value

plotly_ellipse returns an npointsx2 matrix of the points forming the border of the ellipse.

References

Johnson, R. A. and Wichern, D. W. (2002) Applied Multivariate Statistical Analysis, Fifth Edition, Prentice Hall.

See Also

regcr, ellipse

Examples

## Produce a 95% ellipse with the specified mean and covariance structure.
mu <- c(1, 3)
sigma <- matrix(c(1, .3, .3, 1.5), 2, 2)
plotly_ellipse(mu, sigma, npoints = 200)

Plot sequences from the EM algorithm for censored mixture of exponentials using plotly

Description

This is an updated function of plotexpRMM. For more technical details, please refer to plotexpRMM.

Usage

  plotly_expRMM(a , title = NULL , rowstyle = TRUE , subtitle=NULL,
  width = 2 , cex = 2 , col.comp = NULL,
  legend.text = NULL, legend.text.size = 15, legend.size = 15,
  title.x = 0.5, title.y = 0.95, title.size = 15,
  xlab.size = 15, xtick.size = 15, 
  ylab.size = 15, ytick.size = 15)

Arguments

Value

The plot returned

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: expRMM_EM, summary.mixEM, plot.mixEM, plotexpRMM.

Other models and algorithms for censored lifetime data (name convention is model_algorithm): weibullRMM_SEM, spRMM_SEM.

Examples

n=300 # sample size
m=2 # number of mixture components
lambda <- c(1/3,1-1/3); rate <- c(1,1/10) # mixture parameters
set.seed(1234)
x <- rexpmix(n, lambda, rate) # iid ~ exponential mixture
cs=runif(n,0,max(x)) # Censoring (uniform) and incomplete data
t <- apply(cbind(x,cs),1,min) # observed or censored data
d <- 1*(x <= cs) # censoring indicator
###### EM for RMM, exponential lifetimes
l0 <- rep(1/m,m); r0 <- c(1, 0.5) # "arbitrary" initial values
a <- expRMM_EM(t, d, lambda=l0, rate=r0, k = m)
summary(a) # EM estimates etc
plotly_expRMM(a , rowstyle = TRUE) # plot of EM sequences

Visualization of Integrated Squared Error for a selected density from npEM output using plotly

Description

This is an updated visualization function for ise.npEM. For more technical details, please refer to ise.npEM.

Usage

plotly_ise.npEM(npEMout, component=1, block=1, truepdf=dnorm, lower=-Inf,
                upper=Inf, plots = TRUE ,
                col = NULL , width = 3,
                title = NULL , title.size = 15 , title.x = 0.5 , title.y = 0.95,
                xlab = "t" , xlab.size = 15 , xtick.size = 15,
                ylab = "" , ylab.size = 15 , ytick.size = 15,
                legend.text = "" , legend.text.size = 15 , legend.size = 15, ...)

Arguments

Details

This function calls the wkde (weighted kernel density estimate) function.

Value

Just as for the integrate function, a list of class "integrate" with components

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics, 18, 505-526.

• Benaglia, T., Chauveau, D., Hunter, D. R., and Young, D. (2009), mixtools: An R package for analyzing finite mixture models. Journal of Statistical Software, 32(6):1-29.

See Also

npEM, wkde, integrate, ise.npEM

Examples

## Not run: 
data(Waterdata)
set.seed(100)
a <- npEM(Waterdata[,3:10], mu0=3, bw=4) # Assume indep but not iid
plotly_ise.npEM(a , plots = TRUE)

## End(Not run)

Visualization of output of mixEM function using plotly

Description

This is an updated version of plot.mixEM. For more technical details, please refer to plot.mixEM.

Usage

plotly_mixEM(x, 
             loglik = TRUE,
             density = FALSE,
             xlab1="Iteration", xlab1.size=15 , xtick1.size=15,
             ylab1="Log-Likelihood", ylab1.size=15 , ytick1.size=15,
             title1="Observed Data Log-Likelihood", title1.size=15,
             title1.x = 0.5,title1.y=0.95,
             col1="#1f77b4", lwd1=3, cex1=6,
             xlab2=NULL, xlab2.size=15 , xtick2.size=15,
             ylab2=NULL, ylab2.size=15 , ytick2.size=15,
             title2=NULL, title2.size=15,
             title2.x = 0.5,title2.y=0.95, col.hist = "#1f77b4",
             col2=NULL, lwd2=3, cex2=6,
             alpha = 0.05, marginal = FALSE)

Arguments

Value

A plot of the output of mixEM function is presented depends on output type.

See Also

post.beta

Examples

## Not run: 
## EM output for data generated from a 2-component binary logistic regression model.
beta <- matrix(c(-10, .1, 20, -.1), 2, 2)
x <- runif(500, 50, 250)
x1 <- cbind(1, x)
xbeta <- x1
w <- rbinom(500, 1, .3)
y <- w*rbinom(500, size = 1, prob = (1/(1+exp(-xbeta[, 1]))))+
  (1-w)*rbinom(500, size = 1, prob =
                 (1/(1+exp(-xbeta[, 2]))))
out.2 <- logisregmixEM(y, x, beta = beta, lambda = c(.3, .7),
                       verb = TRUE, epsilon = 1e-01)
plotly_mixEM(out.2 , col2 = c("red" , "green") , density = TRUE)

## Fitting randomly generated data with a 2-component location mixture of bivariate normals.
set.seed(100)
x.1 <- rmvnorm(40, c(0, 0))
x.2 <- rmvnorm(60, c(3, 4))
X.1 <- rbind(x.1, x.2)
mu <- list(c(0, 0), c(3, 4))
out.1 <- mvnormalmixEM(X.1, arbvar = FALSE, mu = mu,
                       epsilon = 1e-02)
plotly_mixEM(out.1 , col2 = c("brown" , "blue") ,
             alpha = c(0.01 , 0.05 , 0.1),
             density = TRUE , marginal = FALSE)

## Fitting randomly generated data with a 2-component scale mixture of bivariate normals.
x.3 <- rmvnorm(40, c(0, 0), sigma =
                 matrix(c(200, 1, 1, 150), 2, 2))
x.4 <- rmvnorm(60, c(0, 0))
X.2 <- rbind(x.3, x.4)
lambda <- c(0.40, 0.60)
sigma <- list(diag(1, 2), matrix(c(200, 1, 1, 150), 2, 2))
out.2 <- mvnormalmixEM(X.2, arbmean = FALSE,
                       sigma = sigma, lambda = lambda,
                       epsilon = 1e-02)
plotly_mixEM(out.1 , col2 = c("brown" , "blue") ,
             alpha = c(0.01 , 0.05 , 0.1),
             density = TRUE , marginal = TRUE)

## EM output for simulated data from 2-component mixture of random effects.
data(RanEffdata)
set.seed(100)
x <- lapply(1:length(RanEffdata), function(i)
  matrix(RanEffdata[[i]][, 2:3], ncol = 2))
x <- x[1:20]
y <- lapply(1:length(RanEffdata), function(i)
  matrix(RanEffdata[[i]][, 1], ncol = 1))
y <- y[1:20]
lambda <- c(0.45, 0.55)
mu <- matrix(c(0, 4, 100, 12), 2, 2)
sigma <- 2
R <- list(diag(1, 2), diag(1, 2))
em.out <- regmixEM.mixed(y, x, sigma = sigma, arb.sigma = FALSE,
                         lambda = lambda, mu = mu, R = R,
                         addintercept.random = FALSE,
                         epsilon = 1e-02, verb = TRUE)
plotly_mixEM(em.out , col2 = c("gold" , "purple") , 
             density = TRUE , lwd2 = 1 , cex2 =9)

## Analyzing the Old Faithful geyser data with a 2-component mixture of normals.
data(faithful)
attach(faithful)
set.seed(100)
out <- normalmixEM(waiting, arbvar = FALSE, verb = TRUE,
                   epsilon = 1e-04)
plotly_mixEM(out, density = TRUE , col2 = c("gold" , "purple"))

## EM output for the water-level task data set.
data(Waterdata)
set.seed(100)
water <- t(as.matrix(Waterdata[,3:10]))
em.out <- repnormmixEM(water, k = 2, verb = TRUE, epsilon = 1e-03)
plotly_mixEM(em.out, density = TRUE , col2 = c("gold" , "purple"))

## End(Not run)

Various Plots Pertaining to Mixture Model Output Using MCMC Methods using plotly

Description

This is an updated version of plot.mixMCMC. For technical details, please refer to plot.mixMCMC.

Usage

 
plotly_mixMCMC(x, trace.plot = TRUE, summary.plot = FALSE, burnin = 2000, 
               credit.region = 0.95, col.cr = NULL,
               cex.trace = 3, width.trace = 3, 
               cex.summary = 3, width.summary = 1,
               title.trace = "", title.trace.x = 0.5, 
               title.trace.y = 0.95, title.trace.size = 15,
               xlab.trace = "Index", xlab.trace.size = 15, xtick.trace.size = 15,
               ylab.trace = NULL, ylab.trace.size = 15, ytick.trace.size = 15,
               title.summary = "Credible Regions", title.summary.x = 0.5, 
               title.summary.y = 0.95, title.summary.size = 15,
               xlab.summary = "Predictor", xlab.summary.size = 15, 
               xtick.summary.size = 15,
               ylab.summary = "Response", ylab.summary.size = 15, 
               ytick.summary.size = 15
) 

Arguments

Value

plotly_mixMCMC returns trace plots of the various parameters estimated by the MCMC methods for all objects of class mixMCMC. In addition, other plots may be produced for the following k-component mixture model functions:

See Also

regcr, plot.mixMCMC

Examples

 
## Not run: 
data(NOdata)
attach(NOdata)
set.seed(100)
beta <- matrix(c(1.3, -0.1, 0.6, 0.1), 2, 2)
sigma <- c(.02, .05)
MH.out <- regmixMH(Equivalence, NO, beta = beta, s = sigma,
                   sampsize = 2500, omega = .0013)
plotly_mixMCMC(x = MH.out, summary.plot = TRUE, col.cr = c("red", "green"))

## End(Not run)

Mixturegrams

Description

Construct a mixturegram for determining an apporpriate number of components using plotly.

Usage

plotly_mixturegram(data, pmbs, method=c("pca","kpca","lda"), 
                   all.n=FALSE, id.con=NULL, score=1, iter.max=50, 
                   nstart=25, xlab = "K", xlab.size = 15, 
                   xtick.size = 15, ylab = NULL, ylab.size = 15, 
                   ytick.size = 15, cex = 12, col.dot = "red", 
                   width = 1, title = "Mixturegram", title.size = 15, 
                   title.x = 0.5, title.y = 0.95)

Arguments

Value

plotly_mixturegram returns a mixturegram where the profiles are plotted over component values of k = 1,...,K.

References

Young, D. S., Ke, C., and Zeng, X. (2018) The Mixturegram: A Visualization Tool for Assessing the Number of Components in Finite Mixture Models, Journal of Computational and Graphical Statistics, 27(3), 564–575.

See Also

boot.comp, mixturegram

Examples

## Not run: 
##Data generated from a 2-component mixture of normals.
set.seed(100)
n <- 100
w <- rmultinom(n,1,c(.3,.7))
y <- sapply(1:n,function(i) w[1,i]*rnorm(1,-6,1) +
              w[2,i]*rnorm(1,0,1))
selection <- function(i,data,rep=30){
  out <- replicate(rep,normalmixEM(data,epsilon=1e-06,
                                   k=i,maxit=5000),simplify=FALSE)
  counts <- lapply(1:rep,function(j)
    table(apply(out[[j]]$posterior,1,
                which.max)))
  counts.length <- sapply(counts, length)
  counts.min <- sapply(counts, min)
  counts.test <- (counts.length != i)|(counts.min < 5)
  if(sum(counts.test) > 0 & sum(counts.test) < rep)
    out <- out[!counts.test]
  l <- unlist(lapply(out, function(x) x$loglik))
  tmp <- out[[which.max(l)]]
}
all.out <- lapply(2:5, selection, data = y, rep = 2)
pmbs <- lapply(1:length(all.out), function(i)
  all.out[[i]]$post)
plotly_mixturegram(y, pmbs, method = "pca", all.n = TRUE,
                   id.con = NULL, score = 1,
                   title = "Mixturegram (Well-Separated Data)")

## End(Not run)

Plot Nonparametric or Semiparametric EM Output

Description

This is an updater version of plot.npEM function by using plotly. For technical details, please refer to plot.npEM.

Usage

plotly_npEM(x, blocks = NULL, hist=TRUE, addlegend=TRUE,
            scale = TRUE, title=NULL, breaks="Sturges", 
            dens.col = NULL, newplot=TRUE, ylim = NULL ,
            col.hist = "#1f77b4",
            width = 3, title.x = 0.5 , title.y = 0.95, title.size = 15,
            xlab = "X" , xlab.size = 15 , xtick.size = 15,
            ylab = "Density" , ylab.size = 15 , ytick.size = 15,
            legend.text = "Posteriors",
            legend.text.size = 15,
            legend.size = 15)         
plotly_spEM(x, blocks = NULL, hist=TRUE, addlegend=TRUE,
            scale = TRUE, title=NULL, breaks="Sturges", 
            dens.col = NULL, newplot=TRUE, ylim = NULL ,
            col.hist = "#1f77b4",
            width = 3, title.x = 0.5 , title.y = 0.95, title.size = 15,
            xlab = "X" , xlab.size = 15 , xtick.size = 15,
            ylab = "Density" , ylab.size = 15 , ytick.size = 15,
            legend.text = "Posteriors",
            legend.text.size = 15,
            legend.size = 15)

Arguments

Value

plotly_npEM returns a list with two elements:

See Also

npEM, density.npEM, spEMsymloc, plotseq.npEM, plot.npEM

Examples

 
## Not run: 
## Examine and plot water-level task data set.

## First, try a 3-component solution where no two coordinates are
## assumed i.d.
data(Waterdata)
set.seed(100)
a <- npEM(Waterdata[,3:10], 3, bw=4)
plotly_npEM(a , newplot = FALSE)

## Next, same thing but pairing clock angles that are directly opposite one
## another (1:00 with 7:00, 2:00 with 8:00, etc.)
b <- npEM(Waterdata[,3:10], 3, blockid=c(4,3,2,1,3,4,1,2), bw=4)
plotly_npEM(b , newplot = FALSE)

## End(Not run)

Visualization of Posterior Regression Coefficients in Mixtures of Random Effects Regressions using plotly

Description

Returns a 2x2 matrix of plots summarizing the posterior intercept and slope terms in a mixture of random effects regression with arbitrarily many components.

Usage

plotly_post.beta(y, x, p.beta, p.z,
                 cex = 6,lwd=1,
                 title.size = 15,
                 xlab.size = 15 , xtick.size = 15,
                 ylab.size = 15 , ytick.size = 15,
                 col.data = "#1f77b4",
                 col.comp = NULL)

Arguments

Details

This is primarily used for within plot.mixEM.

Value

Plots returned.

References

Young, D. S. and Hunter, D. R. (2015) Random Effects Regression Mixtures for Analyzing Infant Habituation, Journal of Applied Statistics, 42(7), 1421–1441.

See Also

regmixEM.mixed, plot.mixEM, post.beta

Examples

data(RanEffdata)
set.seed(100)
x <- lapply(1:length(RanEffdata), function(i)
  matrix(RanEffdata[[i]][, 2:3], ncol = 2))
x <- x[1:20]
y <- lapply(1:length(RanEffdata), function(i)
  matrix(RanEffdata[[i]][, 1], ncol = 1))
y <- y[1:20]
lambda <- c(0.45, 0.55)
mu <- matrix(c(0, 4, 100, 12), 2, 2)
sigma <- 2
R <- list(diag(1, 2), diag(1, 2))
em.out <- regmixEM.mixed(y, x, sigma = sigma, arb.sigma = FALSE,
                         lambda = lambda, mu = mu, R = R,
                         addintercept.random = FALSE,
                         epsilon = 1e-02, verb = TRUE)

x.1 = em.out$x
n = sum(sapply(x.1, nrow))
x.1.sum = sum(sapply(1:length(x.1), function(i) length(x.1[[i]][,1])))
if (x.1.sum == n) {
  x = lapply(1:length(x.1), function(i) matrix(x.1[[i]][,-1], ncol = 1))
} else {
  x = x.1
}

plotly_post.beta(x = x, y = em.out$y, p.beta = em.out$posterior.beta, 
                 p.z = em.out$posterior.z)

Plotting sequences of estimates from non- or semiparametric EM-like Algorithm using plotly

Description

This is an updated version of plotseq.npEM. For technical details, please refer to plotseq.npEM.

Usage

  plotly_seq.npEM (x, col = '#1f77b4' , width = 6,
                   xlab = "Iteration" , xlab.size = 15 , xtick.size = 15,
                   ylab.size = 15 , ytick.size = 15,
                   title.size = 15 , title.x = 0.5 , title.y = 0.95)

Arguments

Value

plotly_seq.npEM returns a figure with one plot for each component proportion, and, in the case of spEMsymloc, one plot for each component mean.

Author(s)

Didier Chauveau

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics (to appear).

• Bordes, L., Chauveau, D., and Vandekerkhove, P. (2007), An EM algorithm for a semiparametric mixture model, Computational Statistics and Data Analysis, 51: 5429-5443.

See Also

plot.npEM, rnormmix, npEM, spEMsymloc, plotly_seq.npEM

Examples

## Not run: 
## Examine and plot water-level task data set.
## First, try a 3-component solution where no two coordinates are
## assumed i.d.
data(Waterdata)
set.seed(100)
## Not run:
a <- npEM(Waterdata[,3:10], mu0=3, bw=4) # Assume indep but not iid
plotly_seq.npEM(a)

## End(Not run)

Plot mixture pdf for the semiparametric mixture model output by spEMsymlocN01 using plotly.

Description

This is an updated version of plotlspEMN01 function by using plotly. For technical details, please refer to plot.spEMN01.

Usage

plotly_spEMN01(x, bw=x$bandwidth, knownpdf=dnorm, add.plot=FALSE,
               width = 3 , col.dens = NULL, col.hist =  '#1f77b4',
               title = NULL , title.size = 15 , 
               title.x = 0.5 , title.y = 0.95,
               xlab = "t" , xlab.size = 15 , xtick.size = 15,
               ylab = "Density" , ylab.size = 15 , ytick.size = 15,
               legend.text = "Densities" , legend.text.size = 15 , 
               legend.size = 15)

Arguments

Value

A plot of the density of the mixture

Author(s)

Didier Chauveau

References

• Chauveau, D., Saby, N., Orton, T. G., Lemercier B., Walter, C. and Arrouys, D. Large-scale simultaneous hypothesis testing in soil monitoring: A semi-parametric mixture approach, preprint (2013).

See Also

spEMsymlocN01, plot.spEMN01

Examples

## Probit transform of p-values
## from a Beta-Uniform mixture model
## comparion of parametric and semiparametric EM fit
## Note: in actual situations n=thousands
set.seed(50)
n=300 # nb of multiple tests
m=2 # 2 mixture components
a=c(1,0.1); b=c(1,1); lambda=c(0.6,0.4) # parameters
z=sample(1:m, n, rep=TRUE, prob = lambda)
p <- rbeta(n, shape1 = a[z], shape2 = b[z]) # p-values
o <- order(p)
cpd <- cbind(z,p)[o,] # sorted complete data, z=1 if H0, 2 if H1
p <- cpd[,2] # sorted p-values
y <- qnorm(p) # probit transform of the pvalues
# gaussian EM fit with component 1 constrained to N(0,1)
s1 <- normalmixEM(y, mu=c(0,-4),
                  mean.constr = c(0,NA), sd.constr = c(1,NA))
s2 <- spEMsymlocN01(y, mu0 = c(0,-3)) # spEM with N(0,1) fit
plotly_spEMN01(s2 , add.plot = FALSE)

Plot output from Stochastic EM algorithm for semiparametric scaled mixture of censored data using plotly.

Description

This is an updated version of plotspRMM function. For technical details, please refer to plotspRMM.

Usage

  plotly_spRMM(sem, tmax = NULL,
               width = 3 , col = '#1f77b4', cex = 3,
               title.size = 15 , 
               title.x = 0.5 , title.y = 0.95,
               xlab.size = 15 , xtick.size=15 ,
               ylab.size = 15 , ytick.size=15)

Arguments

Value

The four plots returned.

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: spRMM_SEM , plotspRMM.

Other models and algorithms for censored lifetime data (name convention is model_algorithm): expRMM_EM, weibullRMM_SEM.

Examples

## Not run: 
n=500 # sample size
m=2 # nb components
lambda=c(0.4, 0.6) # parameters
meanlog=3; sdlog=0.5; scale=0.1
set.seed(12)
# simulate a scaled mixture of lognormals
x <- rlnormscalemix(n, lambda, meanlog, sdlog, scale)
cs=runif(n,20,max(x)+400) # Censoring (uniform) and incomplete data
t <- apply(cbind(x,cs),1,min)
d <- 1*(x <= cs)
tauxc <- 100*round( 1-mean(d),3)
cat(tauxc, "percents of data censored.\n")

c0 <- c(25, 180) # data-driven initial centers (visible modes)
sc0 <- 25/180    # and scaling
s <- spRMM_SEM(t, d, scaling = sc0, centers = c0, bw = 15, maxit = 100)

plotly_spRMM(s) # default
summary(s)   # S3 method for class "spRMM"

## End(Not run)

Plot sequences from the Stochastic EM algorithm for mixture of Weibull using plotly

Description

This is an updated version of plotweibullRMM function by using plotly function. For technical details, please refer to plotweibullRMM.

Usage

  plotly_weibullRMM(a, title=NULL, rowstyle=TRUE, subtitle=NULL,
                    width = 3 , col = NULL , 
                    title.size = 15 , title.x = 0.5 , title.y = 0.95,
                    xlab = "Iterations" , xlab.size = 15 , xtick.size = 15,
                    ylab = "Estimates" , ylab.size = 15 , ytick.size = 15,
                    legend.size = 15)

Arguments

Value

The plot returned.

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: weibullRMM_SEM, summary.mixEM, plotweibullRMM.

Other models and algorithms for censored lifetime data (name convention is model_algorithm): expRMM_EM, spRMM_SEM .

Examples

n = 500 # sample size
m = 2 # nb components
lambda=c(0.4, 0.6)
shape <- c(0.5,5); scale <- c(1,20) # model parameters
set.seed(321)
x <- rweibullmix(n, lambda, shape, scale) # iid ~ weibull mixture
cs=runif(n,0,max(x)+10) # iid censoring times
t <- apply(cbind(x,cs),1,min) # censored observations
d <- 1*(x <= cs) # censoring indicator
## set arbitrary or "reasonable" (e.g., data-driven) initial values
l0 <- rep(1/m,m); sh0 <- c(1, 2); sc0 <- c(2,10)
# Stochastic EM algorithm
a <- weibullRMM_SEM(t, d, lambda = l0, shape = sh0, scale = sc0, maxit = 200)
summary(a) # Parameters estimates etc
plotly_weibullRMM(a , legend.size = 20) # plot of St-EM sequences

Plotting sequences of estimates from non- or semiparametric EM-like Algorithm

Description

Returns plots of the sequences of scalar parameter estimates along iterations from an object of class npEM.

Usage

  ## S3 method for class 'npEM'
plotseq(x, ...)

Arguments

Details

plotseq.npEM returns a figure with one plot for each component proportion, and, in the case of spEMsymloc, one plot for each component mean.

References

• Benaglia, T., Chauveau, D., and Hunter, D. R. (2009), An EM-like algorithm for semi- and non-parametric estimation in multivariate mixtures, Journal of Computational and Graphical Statistics (to appear).

• Bordes, L., Chauveau, D., and Vandekerkhove, P. (2007), An EM algorithm for a semiparametric mixture model, Computational Statistics and Data Analysis, 51: 5429-5443.

See Also

plot.npEM, rnormmix, npEM, spEMsymloc

Examples

## Example from a normal location mixture
n <- 200
set.seed(100)
lambda <- c(1/3,2/3)
mu <- c(0, 4); sigma<-rep(1, 2)
x <- rnormmix(n, lambda, mu, sigma)
b <- spEMsymloc(x, mu0=c(-1, 2), stochastic=FALSE)
plotseq(b)
bst <- spEMsymloc(x, mu0=c(-1, 2), stochastic=TRUE)
plotseq(bst)

Plot output from Stochastic EM algorithm for semiparametric scaled mixture of censored data

Description

Function for plotting various results from an object returned by spRMM_SEM, a Stochastic EM algorithm for semiparametric scaled mixture of randomly right censored lifetime data. Four plots of sequences of estimates along iterations, survival and density estimates (see reference below).

Usage

  plotspRMM(sem, tmax = NULL)

Arguments

Value

The four plots returned

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: spRMM_SEM.

Other models and algorithms for censored lifetime data (name convention is model_algorithm): expRMM_EM, weibullRMM_SEM.

Examples

# See example(spRMM_SEM)

Plot sequences from the Stochastic EM algorithm for mixture of Weibull

Description

Function for plotting sequences of estimates along iterations, from an object returned by weibullRMM_SEM, a Stochastic EM algorithm for mixture of Weibull distributions with randomly right censored data (see reference below).

Usage

  plotweibullRMM(a, title = NULL, rowstyle = TRUE, subtitle = NULL, ...)

Arguments

Value

The plot returned

Author(s)

Didier Chauveau

References

• Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7

See Also

Related functions: weibullRMM_SEM, summary.mixEM.

Other models and algorithms for censored lifetime data (name convention is model_algorithm): expRMM_EM, spRMM_SEM .

Examples

n = 500 # sample size
m = 2 # nb components
lambda=c(0.4, 0.6)
shape <- c(0.5,5); scale <- c(1,20) # model parameters
set.seed(321)
x <- rweibullmix(n, lambda, shape, scale) # iid ~ weibull mixture
cs=runif(n,0,max(x)+10) # iid censoring times
t <- apply(cbind(x,cs),1,min) # censored observations
d <- 1*(x <= cs)              # censoring indicator

## set arbitrary or "reasonable" (e.g., data-driven) initial values
l0 <- rep(1/m,m); sh0 <- c(1, 2); sc0 <- c(2,10)
# Stochastic EM algorithm 
a <- weibullRMM_SEM(t, d, lambda = l0, shape = sh0, scale = sc0, maxit = 200)

summary(a) # Parameters estimates etc
plotweibullRMM(a) # default plot of St-EM sequences

EM Algorithm for Mixtures of Poisson Regressions

Description

Returns EM algorithm output for mixtures of Poisson regressions with arbitrarily many components.

Usage

poisregmixEM(y, x, lambda = NULL, beta = NULL, k = 2,
             addintercept = TRUE, epsilon = 1e-08, 
             maxit = 10000, verb = FALSE)

Arguments

Value

poisregmixEM returns a list of class mixEM with items:

References

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

Wang, P., Puterman, M. L., Cockburn, I. and Le, N. (1996) Mixed Poisson Regression Models with Covariate Dependent Rates, Biometrics, 52(2), 381–400.

See Also

logisregmixEM

Examples

## EM output for data generated from a 2-component model.

set.seed(100)
beta <- matrix(c(1, .5, .7, -.8), 2, 2)
x <- runif(50, 0, 10)
xbeta <- cbind(1, x)%*%beta
w <- rbinom(50, 1, .5)
y <- w*rpois(50, exp(xbeta[, 1]))+(1-w)*rpois(50, exp(xbeta[, 2]))
out <- poisregmixEM(y, x, verb = TRUE,  epsilon = 1e-03)
out

Summary of Posterior Regression Coefficients in Mixtures of Random Effects Regressions

Description

Returns a 2x2 matrix of plots summarizing the posterior intercept and slope terms in a mixture of random effects regression with arbitrarily many components.

Usage

post.beta(y, x, p.beta, p.z)

Arguments

Details

This is primarily used for within plot.mixEM.

Value

post.beta returns a 2x2 matrix of plots giving:

References

Young, D. S. and Hunter, D. R. (2015) Random Effects Regression Mixtures for Analyzing Infant Habituation, Journal of Applied Statistics, 42(7), 1421–1441.

See Also

regmixEM.mixed, plot.mixEM

Examples

  ## Not run: 
## EM output for simulated data from 2-component mixture of random effects.

data(RanEffdata)
set.seed(100)
x <- lapply(1:length(RanEffdata), function(i) 
            matrix(RanEffdata[[i]][, 2:3], ncol = 2))
x <- x[1:20]
y <- lapply(1:length(RanEffdata), function(i) 
            matrix(RanEffdata[[i]][, 1], ncol = 1))
y <- y[1:20]
lambda <- c(0.45, 0.55)
mu <- matrix(c(0, 4, 100, 12), 2, 2)
sigma <- 2
R <- list(diag(1, 2), diag(1, 2))
em.out <- regmixEM.mixed(y, x, sigma = sigma, arb.sigma = FALSE,
                         lambda = lambda, mu = mu, R = R,
                         addintercept.random = FALSE,
                         epsilon = 1e-02, verb = TRUE)

## Obtaining the 2x2 matrix of plots.

x.ran <- lapply(1:length(x), function(i) x[[i]][, 2])
p.beta <- em.out$posterior.beta
p.z <- em.out$posterior.z
post.beta(y, x.ran, p.beta = p.beta, p.z = p.z)

## End(Not run)

Printing of Results from the mvnpEM Algorithm Output

Description

print method for class mvnpEM.

Usage

## S3 method for class 'mvnpEM'
print(x, ...)

Arguments

Details

print.mvnpEM prints the elements of an mvnpEM object without printing the data or the posterior probabilities. (These may still be accessed as x$data and x$posteriors.)

Value

print.mvnpEM returns (invisibly) the full value of x itself, including the data and posteriors elements.

See Also

mvnpEM, plot.mvnpEM summary.mvnpEM

Examples

# Example as in Chauveau and Hoang (2015) with 6 coordinates
## Not run: 
m=2; r=6; blockid <-c(1,1,2,2,3,3) # 3 bivariate blocks 
# generate some data x ...
a <- mvnpEM(x, mu0=2, blockid, samebw=F) # adaptive bandwidth
print(a)
## End(Not run)

Printing non- and semi-parametric multivariate mixture model fits

Description

print method for class npEM.

Usage

## S3 method for class 'npEM'
print(x, ...)

Arguments

Details

print.npEM prints the elements of an npEM object without printing the data or the posterior probabilities. (These may still be accessed as x$data and x$posteriors.)

Value

print.npEM returns (invisibly) the full value of x itself, including the data and posteriors elements.

See Also

npEM, plot.npEM summary.npEM

Examples

data(Waterdata)
set.seed(100)
## Not run: npEM(Waterdata[,3:10], 3, bw=4, verb=FALSE) # Assume indep but not iid

Add a Confidence Region or Bayesian Credible Region for Regression Lines to a Scatterplot

Description

Produce a confidence or credible region for regression lines based on a sample of bootstrap beta values or posterior beta values. The beta parameters are the intercept and slope from a simple linear regression.

Usage

regcr(beta, x, em.beta = NULL, em.sigma = NULL, alpha = .05, 
      nonparametric = FALSE, plot = FALSE, xyaxes = TRUE, ...)

Arguments

Value

regcr returns a list containing the following items:

See Also

regmixEM, regmixMH

Examples

## Nonparametric credible regions fit to NOdata. 

data(NOdata)
attach(NOdata)
set.seed(100)
beta <- matrix(c(1.3, -0.1, 0.6, 0.1), 2, 2)
sigma <- c(.02, .05)
MH.out <- regmixMH(Equivalence, NO, beta = beta, s = sigma, 
                   sampsize = 2500, omega = .0013)
attach(data.frame(MH.out$theta))
beta.c1 <- cbind(beta0.1[2400:2499], beta1.1[2400:2499])
beta.c2 <- cbind(beta0.2[2400:2499], beta1.2[2400:2499])
plot(NO, Equivalence)
regcr(beta.c1, x = NO, nonparametric = TRUE, plot = TRUE, 
      col = 2)
regcr(beta.c2, x = NO, nonparametric = TRUE, plot = TRUE, 
      col = 3)

EM Algorithm for Mixtures of Regressions

Description

Returns EM algorithm output for mixtures of multiple regressions with arbitrarily many components.

Usage

regmixEM(y, x, lambda = NULL, beta = NULL, sigma = NULL, k = 2,
         addintercept = TRUE, arbmean = TRUE, arbvar = TRUE, 
         epsilon = 1e-08, maxit = 10000, verb = FALSE)

Arguments

Value

regmixEM returns a list of class mixEM with items:

References

de Veaux, R. D. (1989), Mixtures of Linear Regressions, Computational Statistics and Data Analysis 8, 227-245.

Hurn, M., Justel, A. and Robert, C. P. (2003) Estimating Mixtures of Regressions, Journal of Computational and Graphical Statistics 12(1), 55–79.

McLachlan, G. J. and Peel, D. (2000) Finite Mixture Models, John Wiley and Sons, Inc.

See Also

regcr, regmixMH

Examples

## EM output for NOdata.
 
data(NOdata)
attach(NOdata)
set.seed(100)
em.out <- regmixEM(Equivalence, NO, verb = TRUE, epsilon = 1e-04)
em.out[3:6]

EM Algorithm for Mixtures of Regressions with Local Lambda Estimates

Description

Returns output for one step of an EM algorithm output for mixtures of multiple regressions where the mixing proportions are estimated locally.

Usage

regmixEM.lambda(y, x, lambda = NULL, beta = NULL, sigma = NULL, 
                k = 2, addintercept = TRUE, arbmean = TRUE,
                arbvar = TRUE, epsilon = 1e-8, maxit = 10000,
                verb = FALSE)

Arguments

Details

Primarily used within regmixEM.loc.

Value

regmixEM.lambda returns a list of class mixEM with items:

See Also

regmixEM.loc

Examples

## Compare a 2-component and 3-component fit to NOdata.

data(NOdata)
attach(NOdata)
set.seed(100)
out1 <- regmixEM.lambda(Equivalence, NO)
out2 <- regmixEM.lambda(Equivalence, NO, k = 3)
c(out1$loglik, out2$loglik)

Iterative Algorithm Using EM Algorithm for Mixtures of Regressions with Local Lambda Estimates

Description

Iterative algorithm returning EM algorithm output for mixtures of multiple regressions where the mixing proportions are estimated locally.

Usage

regmixEM.loc(y, x, lambda = NULL, beta = NULL, sigma = NULL, 
             k = 2, addintercept = TRUE, kern.l = c("Gaussian",
             "Beta", "Triangle", "Cosinus", "Optcosinus"), 
             epsilon = 1e-08, maxit = 10000, kernl.g = 0, 
             kernl.h = 1, verb = FALSE) 

Arguments

Value

regmixEM.loc returns a list of class mixEM with items:

See Also

regmixEM.lambda

Examples

## Compare a 2-component and 3-component fit to NOdata.

data(NOdata)
attach(NOdata)
set.seed(100)
out1 <- regmixEM.loc(Equivalence, NO, kernl.h = 2, 
                     epsilon = 1e-02, verb = TRUE)
out2 <- regmixEM.loc(Equivalence, NO, kernl.h = 2, k = 3,
                     epsilon = 1e-02, verb = TRUE)
c(out1$loglik, out2$loglik)

EM Algorithm for Mixtures of Regressions with Random Effects

Description

Returns EM algorithm output for mixtures of multiple regressions with random effects and an option to incorporate fixed effects and/or AR(1) errors.

Usage

regmixEM.mixed(y, x, w = NULL, sigma = NULL, arb.sigma = TRUE,
               alpha = NULL, lambda = NULL, mu = NULL, 
               rho = NULL, R = NULL, arb.R = TRUE, k = 2, 
               ar.1 = FALSE, addintercept.fixed = FALSE, 
               addintercept.random = TRUE, epsilon = 1e-08, 
               maxit = 10000, verb = FALSE)

Arguments

Value

regmixEM returns a list of class mixEM with items:

References

Xu, W. and Hedeker, D. (2001) A Random-Effects Mixture Model for Classifying Treatment Response in Longitudinal Clinical Trials, Journal of Biopharmaceutical Statistics, 11(4), 253–273.

Young, D. S. and Hunter, D. R. (2015) Random Effects Regression Mixtures for Analyzing Infant Habituation, Journal of Applied Statistics, 42(7), 1421–1441.

See Also

regmixEM, post.beta

Examples

## EM output for simulated data from 2-component mixture of random effects.

data(RanEffdata)
set.seed(100)
x <- lapply(1:length(RanEffdata), function(i) 
            matrix(RanEffdata[[i]][, 2:3], ncol = 2))
x <- x[1:20]
y <- lapply(1:length(RanEffdata), function(i) 
            matrix(RanEffdata[[i]][, 1], ncol = 1))
y <- y[1:20]
lambda <- c(0.45, 0.55)
mu <- matrix(c(0, 4, 100, 12), 2, 2)
sigma <- 2
R <- list(diag(1, 2), diag(1, 2))
em.out <- regmixEM.mixed(y, x, sigma = sigma, arb.sigma = FALSE,
                         lambda = lambda, mu = mu, R = R,
                         addintercept.random = FALSE,
                         epsilon = 1e-02, verb = TRUE)
em.out[4:10]

Metropolis-Hastings Algorithm for Mixtures of Regressions

Description

Return Metropolis-Hastings (M-H) algorithm output for mixtures of multiple regressions with arbitrarily many components.

Usage

regmixMH(y, x, lambda = NULL, beta = NULL, s = NULL, k = 2,
         addintercept = TRUE, mu = NULL, sig = NULL, lam.hyp = NULL,
         sampsize = 1000, omega = 0.01, thin = 1)

Arguments

Value

regmixMH returns a list of class mixMCMC with items:

References

Hurn, M., Justel, A. and Robert, C. P. (2003) Estimating Mixtures of Regressions, Journal of Computational and Graphical Statistics 12(1), 55–79.

See Also

regcr

Examples

## M-H algorithm for NOdata with acceptance rate about 40%.

data(NOdata)
attach(NOdata)
set.seed(100)
beta <- matrix(c(1.3, -0.1, 0.6, 0.1), 2, 2)
sigma <- c(.02, .05)
MH.out <- regmixMH(Equivalence, NO, beta = beta, s = sigma, 
                   sampsize = 2500, omega = .0013)
MH.out$theta[2400:2499,]

Model Selection in Mixtures of Regressions

Description

Assess the number of components in a mixture of regressions model using the Akaike's information criterion (AIC), Schwartz's Bayesian information criterion (BIC), Bozdogan's consistent AIC (CAIC), and Integrated Completed Likelihood (ICL).

Usage

regmixmodel.sel(x, y, w = NULL, k = 2, type = c("fixed", 
                "random", "mixed"), ...)

Arguments

Value

regmixmodel.sel returns a matrix of the AIC, BIC, CAIC, and ICL values along with the winner (i.e., the highest value given by the model selection criterion) for various types of regression mixtures.

References

Biernacki, C., Celeux, G. and Govaert, G. (2000) Assessing a Mixture Model for Clustering with the Integrated Completed Likelihood, IEEE Transactions on Pattern Analysis and Machine Intelligence 22(7), 719–725.

Bozdogan, H. (1987) Model Selection and Akaike's Information Criterion (AIC): The General Theory and its Analytical Extensions, Psychometrika 52, 345–370.

See Also

regmixEM, regmixEM.mixed

Examples

## Assessing the number of components for NOdata.

data(NOdata)
attach(NOdata)
set.seed(100)
regmixmodel.sel(x = NO, y = Equivalence, k = 3, type = "fixed")

EM Algorithm for Mixtures of Normals with Repeated Measurements

Description

Returns EM algorithm output for mixtures of normals with repeated measurements and arbitrarily many components.

Usage

repnormmixEM(x, lambda = NULL, mu = NULL, sigma = NULL, k = 2, 
             arbmean = TRUE, arbvar = TRUE, epsilon = 1e-08, 
             maxit = 10000, verb = FALSE)

Arguments

Value

repnormmixEM returns a list of class mixEM with items:

References

Hettmansperger, T. P. and Thomas, H. (2000) Almost Nonparametric Inference for Repeated Measures in Mixture Models, Journal of the Royals Statistical Society, Series B 62(4) 811–825.

See Also

normalmixEM

Examples

## EM output for the water-level task data set.

data(Waterdata)
set.seed(100)
water <- t(as.matrix(Waterdata[,3:10]))
em.out <- repnormmixEM(water, k = 2, verb = TRUE, epsilon = 1e-03)
em.out

Model Selection in Mixtures of Normals with Repeated Measures

Description

Assess the number of components in a mixture model with normal components and repeated measures using the Akaike's information criterion (AIC), Schwartz's Bayesian information criterion (BIC), Bozdogan's consistent AIC (CAIC), and Integrated Completed Likelihood (ICL).

Usage

repnormmixmodel.sel(x, k = 2, ...)

Arguments

Value

repnormmixmodel.sel returns a matrix of the AIC, BIC, CAIC, and ICL values along with the winner (i.e., the highest value given by the model selection criterion) for a mixture of normals with repeated measures.

References

Biernacki, C., Celeux, G., and Govaert, G. (2000). Assessing a Mixture Model for Clustering with the Integrated Completed Likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(7):719-725.

Bozdogan, H. (1987). Model Selection and Akaike's Information Criterion (AIC): The General Theory and its Analytical Extensions. Psychometrika, 52:345-370.

See Also

repnormmixEM