Brief intro: I'm not really sure how to appropriately treat the dependent variables in a set of hierarchical linear models that I'm trying to run. In my models, Level 1 units are children and Level 2 units are schools (we don't have classroom info so there is no nesting level for this).
Background: There are five sets of models, one for each dependent variable. The DVs are all behavioral or cognitive outcome measures that are measured at Level 1. I'm using the HLM 7 software, but I could just as easily use Mplus or R. Each of the five DVs were assessed in the same way. Let's take depression as an example of one: this construct had about 6-8 items (presented to the children as scenarios), and for each item there were 2 response options that could be classified as "the depression options," and 2 response options that could be classified as "not the depression options." The 2 items in each set are treated as equivalent in terms of measurement -- that is, no ordinal ranking from least-to-most depressed. The final form of the depression outcome (as with the other 4 outcomes) is a proportion, which can be interpreted as the proportion of items/scenarios to which the child responded in a way indicative of being depressed. Personally, I hate that the assessment was done in this way, but there's nothing I can do about it now.
Problem: However, I may be able to get permission to score the assessments in a different way, provided I can get access to the item-level data. Right now I only have access to the final proportions. Obviously the default in HLM and in other software is to treat the outcomes as normally distributed. The only distribution appropriate for the [0,1] interval that I know of is the beta, and I do not know of any R packages to implement HLM for beta-distributed variables, nor do I think Mplus or HLM can handle this. I tried using the logit transformation to put the proportions back onto a scale that I could theoretically treat as normal, but with only 6-8 items per DV I still get something that is rather discrete. Not to mention that the 0 and 1 responses, even when converted first to .0001 and .9999 before logit transforming, give me these extreme values that are quite far away from the other transformed values. I've attached histograms to illustrate what I mean.
Question: What are some appropriate ways to treat an outcome variable on the [0,1] interval that takes on only a small handful of observed values when I want to run a set of hierarchical linear models? I'd REALLY like to use the HLM 7 software for this (the P.I. prefers it), but if the best way to treat these types of variables can only be implemented in R then I'm comfortable with that option as well.