A mixed effects model for variances of binomialprobability as a responsee

Dear Nimrod,

Just my 2 cents (even less than 2):

On Thu, 01-01-1970, at 01:00, nimrod.rubinstein <
nimrod.rubinstein at gmail.com> wrote:

Hi all,

I have the following data:

50 subjects that belong to 2 groups (25 in each group). For each subject

I

have replicated measurements (in the range of 15-30 replicates per

subject)

of the number of trails as well as the number successes out of these
trails. Each single trial in each single replicate for each subject is
assumed to be a Bernoulli experiment, therefore the number of successes

in

each replicate is assumed to be a binomially distributed random variable.
To complicate things even more, for each subject from each group I have

not

1 but 3 series of these replicated measurements, from 3 different
environments (the 3 environments are the same for all subjects).

A few example lines of how my data looks like:
subject group environment     number.trails  number.successes

    S1    G1          E1 10100,13209,17621  8500,11622,14684

    S1    G1          E2 13322,12225,10024  10351,10222,7951

    S1    G1          E3  9812,13500,13214  7245,10622,10684

    .

    .

    .

    S50    G2          E1 60547,75345,68412 12014,13214,15244

    S50    G2          E2 50897,62541,61247 10141,13251,30645

    S50    G2          E3 59214,84264,74156 17245,16425,11123


Where the comma separated number of successes correspond to the comma
separated number of trials in each replicate.

I want to test whether the variances of the success probabilities of
subjects from group G1 are significantly smaller than the variances of

the

success probabilities of subjects from group G2.

But, as you say in point 1. below (and in another reply), if this is really
binomial, then you already know that if probabilities in G1 are different
from probabilities in G2 your variances will differ.

So I am not sure I even understand this hypothesis/question.

In other words, for each subject in each replicate a success probability
can be estimated from the observed data. Then the variance of those

success

probabilities over all replicates from a specific environment can be
estimated. I would like to know whether these variances are smaller in
subjects from G1 vs. subjects from G2.

The issues are:


   1. The variance of a success probability depends on the success
   probability itself (natural to binomially distributed data). I need to
   account for this effect.
   2. The number of trails (although at least in the hundreds for all
   subjects) varies several orders of magnitude among subjects. That by

itself

   affects the accuracy of estimating the success probability, but in
   addition, I cannot rule out an inherent dependency between the number

of

   observed trials and the success probability (there is no control over

the

   number of trials in the experiment - it is observed). In other words,
   having a lower number of trails by itself may affect the success
   probability.

That last issue seems very problematic to me. If I understand correctly, we
could think of a perverse situation where, say, you can never observe low
success probabilities if those are related to number of trials around 0.


Can you do a simple plot of probability of success (or its logit) vs number
of trials, to rule out perverse scenarios? Is probability of success always
larger than 0? And number of trials always sufficiently different from 0?


And you say "varies by subject", so the number of trials is always the same
for a subject? And the potential relationship prob <-> number of trials is
due to among-subject variation? So you could approach the problem via a
random effect for subject? If you compare G1 and G2 for number of trials:
are there differences?

I though a simple way to test my hypothesis is to compute the success
probability in each replicate as number.successes/number.trails, and then
the sample mean and sample variance of these success probabilities over

all

replicates, as well as the mean number of trials over all replicates

(added

to my data as sample.mean, sample.var and mean.trials fields). Then I can
use the lmer function in the lme4 package this way:

null.model = lmer(sample.var ~ sample.mean + mean.trials +

(1|environment),

my.data, REML = F)
alt.model = lmer(sample.var ~ sample.mean + mean.trials + group +
(1|environment), my.data, REML = F)

So fixed effects are: sample.mean, mean.trials and group, and random
effects is environment since the three environments in my data are

simply a

random sample from many other possible environments.

I do not see how this would answer your questions, given the above
concerns/issues and other possible approaches (just compare prob. between
G1 and G2, accounting for environment).


Best,


R.

To test whether group has a significant effect I'll do:

anova(null.model, alt.model).

My concerns are:


   1. Should I add subject as a random effect? My intuition is not to

since

   subjects are unique to groups, which is the effect I'm after.

However, I do

   want to account for the fact that sample variances of the same

subject from

   the three different environments are correlated.
   2. The response - sample variances, cannot be assumed to follow a

normal

   distribution (they only take positive values to begin with and as far

as I

   know (n-1)*sample.variance/population.variance ~ chi.square with df =

n-1)

   and therefore the assumptions of a linear mixed effects model do not
   strictly hold. Is there a generalized linear model that would be
   appropriate for this?

Thanks a lot,
nimrod

      [[alternative HTML version deleted]]

--
Ramon Diaz-Uriarte
Department of Biochemistry, Lab B-25
Facultad de Medicina
Universidad Aut?noma de Madrid
Arzobispo Morcillo, 4
28029 Madrid
Spain

Phone: +34-91-497-2412

Email: rdiaz02 at gmail.com
       ramon.diaz at iib.uam.es

http://ligarto.org/rdiaz