4
$\begingroup$

I am trying to fit a finite mixture model to a dependent variable which is bounded (practically) between -0.594 and 1 (theoretically, the latent variable is bounded between -Inf - 1). The data are also bimodal, with a large number of values at '1'. The objective of the analysis is prediction of the dependent variable.

My current approach has been to fit a mixture of normal distributions using the flexmix package in R, but I'd really like to account for the bounded nature of the data, as a recent study found this to be important (I also choose k=3 components based on this study). Using flexmix for truncated data appears non-trivial, as suggested here.

Is there an R package that will permit mixture models with bounded data? I've noticed that actually predicted values do not seem to fall outside the bounded range; i.e. predicted values are not in practice greater than 1. Is this just a fluke of my data, or is it a feature of the methods I've used? Is the bounding even a problem in this context?

As an alternative, I've tried transforming the data by simply taking 1-the dependent variable, thereby giving me a (zero-inflated) variable bounded by 0 and Inf which I have tried to model as a mixture of zero-inflated poisson models but I get the error:

Error in FLXfit(model = model, concomitant = concomitant, control = control,  
: 1 Log-    likelihood: NaN

Is it possible to model non-integers with the poisson family in this context? Any suggestions or thoughts would be greatly appreciated, I'm very new to mixture modelling and indeed GLMs etc.

Here's some simulated data: https://dl.dropbox.com/u/65336009/mydata.csv

Here's my code:

require(flexmix)
require(ggplot2)
mydata <- data.frame(read.csv("mydata.csv", head=T))
attach(mydata)

#Plot of y var
summary(y)
ggplot(mydata, aes(y)) + geom_histogram(binwidth = .1)

#Simplified example of my current 'best' approach####
m1 <- flexmix(y ~ x1 + x2 + x3,
              data = mydata,
              k = 3)

#Predict cluster membership
clusters <- data.frame(clusters(m1, newdata = mydata))

#Predict y
a <- data.frame(predict(m1, newdata = mydata))

#Select prediction based on predicted cluster membership
mydata$flexmix.norm <- ifelse(clusters[,1]==1, a[,1],
                                   ifelse(clusters[,1] == 2,
                                          a[,2], a[,3]))
    print(max(mydata$flexmix.norm))

#Plot predicted values
ggplot(mydata, aes(flexmix.norm)) + geom_histogram(binwidth = .1)

#Maybe it's more natural to model as 1 - y, which is bounded (0,Inf) ####
y.d <- 1 - y
ggplot(mydata, aes(y.d)) + geom_histogram(binwidth = .1)

#Error here ***
m2 <- flexmix(y.d ~ x1 + x2 + x3,
              data = mydata,
              k = 3,
              model=FLXMRziglm(family="poisson"))
rm2 = refit(m2)

#Predict cluster membership
clusters <- NULL
clusters <- data.frame(clusters(m2, newdata = mydata))

#Predict y (note back on original scale of y)
b <- 1 - data.frame(predict(m2, newdata = mydata))

#Select prediction based on predicted cluster membership
preds$flexmix.pois <- ifelse(clusters[,1]==1, b[,1],
                              ifelse(clusters[,1] == 2,
                                     b[,2], b[,3]))

ggplot(mydata, aes(flexmix.pois)) + geom_histogram(binwidth = .1)

Thanks

$\endgroup$
  • $\begingroup$ Maybe you can set k = 2 here. I also encounter this error when I try fitting model with k value it cannot converge. $\endgroup$ – Shixiang Wang Nov 14 '18 at 2:36
  • $\begingroup$ If it is normally distributed, then it isn't bounded. If its bounded, then it isn't normally distributed. One could say for a standard normal distribution that it is practically bounded somewhere starting abs(x) >4.5. Standard EM with non-pathological data, and the correct number of components shouldn't launch the mean of your fit outside a range like that, so explicitly including in the fit doesn't give an advantage that I can see. $\endgroup$ – EngrStudent Feb 9 at 21:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.