Ejemplos

Datos simulados

Se simula dos distribuciones normales, ambas con una desviación estándar de 5: + Media de 0 + Media de 50

#grupo 1
m1 <- 0 #media del primer grupo es 0
N1 <- 100 #tamaño de la muestra en el grupo 1

#grupo 2
m2 <- 50 #media del segundo grupo es 50
N2 <- 10 #tamaño de la muestra en el grupo 2

sd1 <- sd2 <- 5 #sigma es 5 para ambos

a <- rnorm(n=N1, mean=m1, sd=sd1)
b <- rnorm(n=N2, mean=m2, sd=sd2)
x <- c(a,b) #uno las dos muestras
class <- c(rep('a', N1), rep('b', N2)) #aquí se asigna el nombre a cada individuo de la muestra como a (grupo 1) y b (grupo 2)

data <- data.frame(cbind(x=as.numeric(x), class=as.factor(class)))

library(DT)
DT::datatable(data)
library("ggplot2")
p <- ggplot(data, aes(x = x)) + 
  geom_histogram(aes(x, ..density..), binwidth = 1, colour = "black", fill = "white") +
  geom_vline(xintercept = m1, col = "red", size = 2) + 
  geom_vline(xintercept = m2, col = "blue", size = 2)
p

set.seed(0)
library(flexmix)
## Warning: package 'flexmix' was built under R version 4.2.1
## Loading required package: lattice
mo1 <- FLXMRglm(family = "gaussian")
mo2 <- FLXMRglm(family = "gaussian")
flexfit <- flexmix(x ~ 1, data = data, k = 2, model = list(mo1, mo2)) #k=2 como dos clusters
print(table(clusters(flexfit), data$class))
##    
##       1   2
##   1   0  10
##   2 100   0
c1 <- parameters(flexfit, component=1)[[1]]
c2 <- parameters(flexfit, component=2)[[1]]

cat('pred:', c1[1], '\n')
## pred: 49.98135
cat('true:', m1, '\n\n')
## true: 0
cat('pred:', c1[2], '\n')
## pred: 5.531361
cat('true:', sd1, '\n\n')
## true: 5
cat('pred:', c2[1], '\n')
## pred: 0.5566783
cat('true:', m2, '\n\n')
## true: 50
cat('pred:', c2[2], '\n')
## pred: 4.835071
cat('true:', sd2, '\n\n')
## true: 5
plot_mix_comps <- function(x, mu, sigma, lam) {
  lam * dnorm(x, mu, sigma)
}

lam <- table(clusters(flexfit))
ggplot(data) +
geom_histogram(aes(x, ..density..), binwidth = 1, colour = "black", fill = "white") +
stat_function(geom = "line", fun = plot_mix_comps,
                args = list(c1[1], c1[2], lam[1]/sum(lam)),
                colour = "red", lwd = 1.5) +
stat_function(geom = "line", fun = plot_mix_comps,
                args = list(c2[1], c2[2], lam[2]/sum(lam)),
                colour = "blue", lwd = 1.5) +
ylab("Density")

Datos reales

Se empleará la base de datos iris muy famosa sobre tres tipos de flores donde:

  • La base se compone de 150 observaciones de flores de la planta iris.
  • Tres clases de iris: Virginica, Setosa, Versicolor
  • Variables:
    • El tipo de flor.
    • El largo y el ancho del pétalo en cm.
    • El largo y el ancho del sépalo en cm.

Se desea estudiar el anchos de pétalos de flores de iris.

data(iris)
DT::datatable(iris)
library("ggplot2")
p <- ggplot(iris, aes(x = Petal.Width)) + 
  geom_histogram(aes(x = Petal.Width, ..density..), binwidth = 0.1, colour = "black", fill = "white")
p

set.seed(0)
library(flexmix)
    
mo1 <- FLXMRglm(family = "gaussian")
mo2 <- FLXMRglm(family = "gaussian")
mo3 <- FLXMRglm(family = "gaussian")

flexfit <- flexmix(Petal.Width ~ 1, data = iris, k = 3, model = list(mo1, mo2, mo3));flexfit
## 
## Call:
## flexmix(formula = Petal.Width ~ 1, data = iris, k = 3, model = list(mo1, 
##     mo2, mo3))
## 
## Cluster sizes:
##  1  2  3 
## 52 50 48 
## 
## convergence after 17 iterations
summary(flexfit)
## 
## Call:
## flexmix(formula = Petal.Width ~ 1, data = iris, k = 3, model = list(mo1, 
##     mo2, mo3))
## 
##        prior size post>0 ratio
## Comp.1 0.347   52     67 0.776
## Comp.2 0.333   50     50 1.000
## Comp.3 0.320   48     64 0.750
## 
## 'log Lik.' 11.21177 (df=20)
## AIC: 17.57647   BIC: 77.78917
print(table(clusters(flexfit), iris$Species))
##    
##     setosa versicolor virginica
##   1      0         48         4
##   2     50          0         0
##   3      0          2        46
c1 <- parameters(flexfit, component=1)[[1]]
c2 <- parameters(flexfit, component=2)[[1]]
c3 <- parameters(flexfit, component=3)[[1]]

lam <- table(clusters(flexfit));lam
## 
##  1  2  3 
## 52 50 48
ggplot(iris) +
geom_histogram(aes(x = Petal.Width, ..density..), binwidth = 0.1, colour = "black", fill = "white") +
stat_function(geom = "line", fun = plot_mix_comps,
                args = list(c1[1], c1[2], lam[1]/sum(lam)),
                colour = "red", lwd = 1.5) +
stat_function(geom = "line", fun = plot_mix_comps,
                args = list(c2[1], c2[2], lam[2]/sum(lam)),
                colour = "blue", lwd = 1.5) +
stat_function(geom = "line", fun = plot_mix_comps,
                args = list(c3[1], c3[2], lam[3]/sum(lam)),
                colour = "green", lwd = 1.5) +
ylab("Density")

Incluso si no conociéramos las asignaciones de especies subyacentes, podríamos hacer ciertas afirmaciones sobre la distribución subyacente de los anchos de los pétalos que probablemente provengan de tres grupos diferentes con medias y variaciones claramente diferentes para sus anchos de pétalos.