Fit Mixture Distribution to Data and Evaluate Results

mx_metafit(breaks, counts, parms, sd_min = 0, save.data = TRUE, ...)

Arguments

breaks	class boundaries of the data
counts	frequency of observations
parms	list of initial parameters for the mixture, or a suitable `mxObj`ect with start parameters
sd_min	lower boundary value for standard deviation and rate parameter
save.data	if data should be included in the returned object
...	additional arguments passed to `mle2` and `optim`

Examples

## === preparation of data ===

zd <- 6:35
cutoff <- 6.0

counts <- c(155, 0, 8, 12, 17, 35, 37, 66, 42, 39, 13, 4, 4, 8, 19, 36,
  80, 205, 188, 219, 170, 104, 32, 13, 12, 0, 3, 0, 0, 0)

breaks <- c(zd, (max(zd+1))) - cutoff

observations <- unbin(zd, counts) - cutoff

## === quick example ===

(comp <- mx_guess_components(observations, bw=2/3, mincut=0.9))
#>          mean       sd         L
#> 1  0.03579815 0.700029 0.1043107
#> 2  6.81950690 2.022370 0.1772299
#> 3 18.50339058 2.143380 0.7184594
obj <- mxObj(comp, left="e")

ret <- mx_metafit(breaks, counts, obj)
summary(ret)
#> Maximum likelihood estimation
#> 
#> Call:
#> mle2(minuslogl = llunimix, start = parms)
#> 
#> Coefficients:
#>         Estimate Std. Error  z value     Pr(z)    
#> L1     0.1017448  0.0077600  13.1114 < 2.2e-16 ***
#> L2     0.1821405  0.0099684  18.2718 < 2.2e-16 ***
#> rate1  6.9607008  2.6184206   2.6584  0.007852 ** 
#> mean2  7.3083918  0.1311111  55.7420 < 2.2e-16 ***
#> sd2    2.0741827  0.1031195  20.1144 < 2.2e-16 ***
#> mean3 19.0217520  0.0614574 309.5114 < 2.2e-16 ***
#> sd3    1.9838012  0.0453060  43.7867 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 8162.154 

## === details ===

comp <- mx_guess_components(observations, bw=2/3, mincut=0.9)
comp # parameters of components in form of a data frame
#>          mean       sd         L
#> 1  0.03579815 0.700029 0.1043107
#> 2  6.81950690 2.022370 0.1772299
#> 3 18.50339058 2.143380 0.7184594

obj <- mxObj(comp, left="n") ## all components normal
obj <- mxObj(comp, left="e") ## left component exponential (= default)

## base function, needs parms in correct format
ret <-fit_unimix(breaks, counts, parms=pstart(obj), type="enn")

## for more difficult cases, slower
ret <-fit_unimix(breaks, counts, parms=pstart(obj), type="enn",
        method="BFGS",
        control=list(maxit=1000, ndeps=rep(1e-4, length(pstart(obj)))))

## higher level function, determines type from start parameters
ret <- mx_metafit(breaks, counts, obj)

## --- evaluation of results --
coef(ret)         # top-level mxObj-object, contains all weights
#> $e1
#> $e1$type
#> [1] "e"
#> 
#> $e1$L
#> [1] 0.1017448
#> 
#> $e1$rate
#> [1] 6.960701
#> 
#> 
#> $n2
#> $n2$type
#> [1] "n"
#> 
#> $n2$L
#> [1] 0.1821405
#> 
#> $n2$mean
#> [1] 7.308392
#> 
#> $n2$sd
#> [1] 2.074183
#> 
#> 
#> $n3
#> $n3$type
#> [1] "n"
#> 
#> $n3$mean
#> [1] 19.02175
#> 
#> $n3$sd
#> [1] 1.983801
#> 
#> $n3$L
#> [1] 0.7161147
#> 
#> 
coef(ret@fit)     # included mle2-object, last weight missing (sums up to 1.0)
#>         L1         L2      rate1      mean2        sd2      mean3        sd3 
#>  0.1017448  0.1821405  6.9607008  7.3083918  2.0741827 19.0217520  1.9838012 

summary(ret)      # from bbe::mle
#> Maximum likelihood estimation
#> 
#> Call:
#> mle2(minuslogl = llunimix, start = parms)
#> 
#> Coefficients:
#>         Estimate Std. Error  z value     Pr(z)    
#> L1     0.1017448  0.0077600  13.1114 < 2.2e-16 ***
#> L2     0.1821405  0.0099684  18.2718 < 2.2e-16 ***
#> rate1  6.9607008  2.6184206   2.6584  0.007852 ** 
#> mean2  7.3083918  0.1311111  55.7420 < 2.2e-16 ***
#> sd2    2.0741827  0.1031195  20.1144 < 2.2e-16 ***
#> mean3 19.0217520  0.0614574 309.5114 < 2.2e-16 ***
#> sd3    1.9838012  0.0453060  43.7867 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> -2 log L: 8162.154 

AIC(ret)
#> [1] 8176.154
logLik(ret)
#> 'log Lik.' -4081.077 (df=7)
vcov(ret)          # covariance matrix
#>                  L1            L2         rate1         mean2           sd2
#> L1     6.021791e-05 -1.233938e-05 -1.524291e-06  3.611191e-06 -8.170295e-06
#> L2    -1.233938e-05  9.936889e-05  1.489036e-06  3.465900e-05  5.865620e-05
#> rate1 -1.524291e-06  1.489036e-06  6.856126e+00 -2.893887e-05  1.621884e-06
#> mean2  3.611191e-06  3.465900e-05 -2.893887e-05  1.719011e-02  1.322799e-03
#> sd2   -8.170295e-06  5.865620e-05  1.621884e-06  1.322799e-03  1.063363e-02
#> mean3 -2.482953e-07  1.062269e-05 -2.880115e-07  3.137907e-04  4.063190e-04
#> sd3    3.221778e-07 -1.440805e-05  3.854576e-07 -4.190744e-04 -5.305539e-04
#>               mean3           sd3
#> L1    -2.482953e-07  3.221778e-07
#> L2     1.062269e-05 -1.440805e-05
#> rate1 -2.880115e-07  3.854576e-07
#> mean2  3.137907e-04 -4.190744e-04
#> sd2    4.063190e-04 -5.305539e-04
#> mean3  3.777007e-03 -1.205075e-04
#> sd3   -1.205075e-04  2.052633e-03

cov2cor(vcov(ret)) # correlation matrix
#>                  L1            L2         rate1         mean2           sd2
#> L1     1.000000e+00 -1.595164e-01 -7.501802e-05  3.549345e-03 -1.021020e-02
#> L2    -1.595164e-01  1.000000e+00  5.704801e-05  2.651866e-02  5.706214e-02
#> rate1 -7.501802e-05  5.704801e-05  1.000000e+00 -8.429519e-05  6.006751e-06
#> mean2  3.549345e-03  2.651866e-02 -8.429519e-05  1.000000e+00  9.783941e-02
#> sd2   -1.021020e-02  5.706214e-02  6.006751e-06  9.783941e-02  1.000000e+00
#> mean3 -5.206331e-04  1.733945e-02 -1.789767e-06  3.894277e-02  6.411395e-02
#> sd3    9.163832e-04 -3.190247e-02  3.249239e-06 -7.054987e-02 -1.135620e-01
#>               mean3           sd3
#> L1    -5.206331e-04  9.163832e-04
#> L2     1.733945e-02 -3.190247e-02
#> rate1 -1.789767e-06  3.249239e-06
#> mean2  3.894277e-02 -7.054987e-02
#> sd2    6.411395e-02 -1.135620e-01
#> mean3  1.000000e+00 -4.327975e-02
#> sd3   -4.327975e-02  1.000000e+00