2
$\begingroup$

I have been using this site for the past year to get tips on stats but I have run into a problem that I cannot overcome. It isn't too complicated, but I am not sure how to determine significance or if my ideas are correct. Hopefully some of you can help out a bit. Thanks ahead of time.

I have a situation where I have two strains of mice. In one strain of mouse, they have a high incidence of a disease while another has a lower incidence of the same disease. The disease is binary, either they have it or they do not.

If I treat a population of either mouse strain with Drug A, they both get a reduction in the incidence of the disease. The high incidence mouse strain has the disease incidence drop from 30% -> 15%. The low incidence mouse strain sees a drop from 10% -> 1% after treatment.

From this experiment, if I looked the absolute drop in the incidence it would appear that the drug is more effective in the high incidence group that has a decrease of 15%, compared to 9% in the other. However, (to me) it is clear that the drug is far more effective in the low incidence mouse strain because the fold change is a 10-fold reduction, while the high incidence mice only have a 2-fold reduction.

Hopefully the scenario made sense... My question is:

If I want to determine if the fold reduction in the lower incidence mice is significantly greater that in the high incidence, what kind of statistical test should I do?

I am thinking that I should log(2) transform the incidences and then test the effect of treatment with a linear model or ANOVA. I am not sure if this is kosher, especially with my binary response variable and no replicates (impossible to do replicates due to large number of mice ~3000).

Any of your suggestions/answers would be great and highly appreciated.

$\endgroup$
  • 1
    $\begingroup$ You have to decide what you mean by "more effective," or whether even stating your findings that way really helps others understand your results. You seem to have a situation where the absolute benefit is greater for the high-incidence group while the fold difference is greater for the low-incidence group. There's nothing wrong with that. If by "more effective" you mean the numbers that would be helped by giving everyone in the population the drug, then it's more effective in the high-incidence group. Think about that issue before you worry about the statistical test. $\endgroup$ – EdM May 30 '15 at 14:00
  • $\begingroup$ Thanks alot for the response EdM. So by more effective, I mean that the treatment is able to reduce the incidence (or risk of disease) to close to zero. Essentially in the low incidence group, almost all of the sick would be cured compared with only half of those in the high incidence group. In reality, my situation is more complicated than this, but I simplified it because it would be very time consuming to explain. $\endgroup$ – Ehi Akhirome May 30 '15 at 16:48
  • $\begingroup$ Assuming what I mentioned above, would it be possible to compare fold changes in a statistical manner? If so, what type of tests would I do? Just a point toward the right direction would be awesome. Thanks $\endgroup$ – Ehi Akhirome May 30 '15 at 16:51
  • $\begingroup$ It sounds like you are comparing the response of one treatment on two different strains of mouse. when you say "disease" is true/false then your fundamental distribution is binomial. Check out the following: cookbook-r.com/Statistical_analysis/Logistic_regression $\endgroup$ – EngrStudent - Reinstate Monica Jun 1 '15 at 16:51
1
$\begingroup$

This type of situation with binary outcomes is usually handled by logistic regression, which in your case would use the disease/no-disease states of the 3000 mice to estimate the log of the odds of having disease as a function of predictor variables. Here, you would include the mouse strain and the drug treatment as predictor variables; also, since you are interested in whether the drug is more "effective" in one mouse strain, you would include an interaction term between strain and drug. A significant interaction term would support the hypothesis that drug responses, in terms of influence on log-odds of having disease, differ between mouse strains. The results expressed in terms of log-odds can be translated back to probability scales and thus to fold differences in incidence of disease if you wish.

That said, I suggest that you think some more about your definition of "more effective"; I'll translate from mice to people to put this into a public health perspective. Say that there were 1000 people in a high-incidence group and 1000 in a low-incidence group. Say that the drug costs \$1 per dose, and you only had \$1000 dollars to spend.

Which group should you spend the money on? You would prevent disease in 150 people by giving it to the high-incidence group, but only in 90 people if you give it to the low-incidence group. From this cost-benefit perspective, the drug is more effective if used on the high-incidence group, even if the fold-change in incidence is greater in the low-incidence group.

$\endgroup$
  • 1
    $\begingroup$ It might be useful to compare what the different link functions would imply about the parameter being estimated. Log link gives risk ratios, linear link gives absolute changes in risk, and logit link gives log odds. The OP is probably most interested in the log link. $\endgroup$ – Andrew M Jun 1 '15 at 16:57
  • $\begingroup$ Thanks for the great answer! I think the log odds are what I need. I have done logistic model of my data and found that an interaction effect between strain * treatment is a significant factor. In my understanding, this is saying that the natural log of my treatment effect is significantly different between strains. If I am correct in saying that, then I believe that the logistic model is saying that the ln(absolute treatment effect) is different between mouse strains? By absolute treatment effect I mean 15% decrease and 9% decrease in the high and low incidence groups, respectively $\endgroup$ – Ehi Akhirome Jun 4 '15 at 4:38
  • $\begingroup$ Is there a way to modify the model to assess ln(relative fold change) of the treatment effect? By relative fold change, I mean the 2 fold and 10 fold decrease in the high and low incidence groups, respectively. Thanks again for all the help guys $\endgroup$ – Ehi Akhirome Jun 4 '15 at 4:44
  • $\begingroup$ There is no need to re-analyze the data in terms of fold changes (and I'm not sure how to do that directly). You can simply report the relative fold changes, citing the logistic regression as support for your contention that the effects of the drug differ between strains. If you want error estimates on fold changes, you can try translating the logistic regression coefficients and their standard errors from log-odds scales to probability scales. This requires some care, but would be a good way for you to get a solid grasp on the meanings of those coefficients. $\endgroup$ – EdM Jun 4 '15 at 14:25
  • $\begingroup$ Thanks alot EdM, the error estimates on the fold changes are exactly what I want! I will try what you said and let you know how it goes. $\endgroup$ – Ehi Akhirome Jun 4 '15 at 15:30
1
$\begingroup$

I translate your question into:

  • If I want to determine the difference in response to the same treatment of two species of mice, what is the correct approach?
  • how do I test if my idea (about the effectiveness) is correct?

Thoughts:

  • I think that this is a good places for Bayes rule. (Yudkowsky Bayes)
  • I also think that you should look at confidence intervals and sample sizes.
  • Bayes would come after the CI given sample sizes, imo.

My first question would be: "What is the lower 95% confidence interval (CI) given my sample size."

If I have 10 mice of each species then the ratios mean something substantially different than for 100 mice. Here is a great calculator for approaching the question: link. The column that I recommend is Jeffreys.

If I have 10 mice and 3 of them have the disease then

  • my lower 95% CI is 9.27%
  • my lower 95% CI is 60.58%
  • my mean is 30.00%

This means that there is good chance that the "true" value is somewhere between 9% and 61%, and poor chance that it is outside of these values given the data. The region where "true" is expected to live is about 50% - half the entire domain.

If I have 100 mice and 30 of them have the disease then

  • my lower 95% CI is 21.68%
  • my lower 95% CI is 39.45%
  • my mean is 30.00%

This means that there is good chance that the "true" value is somewhere between 21% and 40%, and poor chance that it is outside of these values given the data. The region where "true" is expected to live is 20% - one fifth the entire domain. It took 10x the samples, but reduced the domain by 2.5x.

You can see how higher sample rate substantially pulls in the width. You can also see how with only 10 samples, a value of 10% is interior to the 95% CI. That means if you only had 10 samples of each breed per treatment, the experimental approach might have a problem. It might not give enough data to clearly differentiate the two.

My second question would be: "If I have a mouse of a given species, and it does not have the disease, and I have equal counts of both mouse species in my population, what is the likelihood that it was given a treatment".

Here is some code using Bayes rule (hopefully not incorrectly):

    #set random seed
set.seed(1)

#total number of mice over 4
#    each treat, each species is this value
N <- 2000

#roll some dice
p30t      <- rbinom(N,1,1-0.1)       #treated species30pct undiseased
p30ut     <- rbinom(N,1,1-0.3)      #untreated species30pct undiseased
p10t      <- rbinom(N,1,1-0.01)      #treated species10pct
p10ut     <- rbinom(N,1,1-0.1)       #untreated species10pct

#some sums 
num_treated    <- 2*N
num_untreated  <- 2*N
num_species10  <- 2*N
num_species30  <- 2*N
num_total      <- 4*N

# prob of species10 AND undiseased given treated
p_s10_ud_gt <- (sum(p10t))/num_treated
p_s10_ud_gt

# prob(treated)
p_treated <- num_treated/(num_treated+num_untreated)
p_treated

# prob of species10 AND undiseased 
p_s10_ud <- (sum(p10t)+sum(p10ut))/num_total
p_s10_ud

#prob treated given s10 and undiseased
p_t_gs10_gud <- (p_s10_ud_gt*p_treated)/p_s10_ud
p_t_gs10_gud

####

# prob of species30 AND undiseased given treated
p_s30_ud_gt <- (sum(p30t))/num_treated
p_s30_ud_gt

# prob(treated)
p_treated <- num_treated/(num_treated+num_untreated)
p_treated

# prob of species30 AND undiseased 
p_s30_ud <- (sum(p30t)+sum(p30ut))/num_total
p_s30_ud

#prob treated given s30 and undiseased
p_t_gs30_gud <- (p_s30_ud_gt*p_treated)/p_s30_ud
p_t_gs30_gud

#### ####

# prob of species10 AND diseased given treated
p_s10_d_gt <- (sum(1-p10t))/num_treated
p_s10_d_gt

# prob(treated)
p_treated <- num_treated/(num_treated+num_untreated)
p_treated

# prob of species10 AND diseased 
p_s10_d <- (sum(1-p10t)+sum(1-p10ut))/num_total
p_s10_d

#prob treated given s10 and diseased
p_t_gs10_gd <- (p_s10_d_gt*p_treated)/p_s10_d
p_t_gs10_gd

####

# prob of species30 AND diseased given treated
p_s30_d_gt <- (sum(1-p30t))/num_treated
p_s30_d_gt

# prob(treated)
p_treated <- num_treated/(num_treated+num_untreated)
p_treated

# prob of species30 AND diseased 
p_s30_d <- (sum(1-p30t)+sum(1-p30ut))/num_total
p_s30_d

#prob treated given s30 and undiseased
p_t_gs30_gd <- (p_s30_d_gt*p_treated)/p_s30_d
p_t_gs30_gd

When I put this into a table and compare it with the analytic result I get the following:

enter image description here

where:

enter image description here

You can see where I manually entered the numeric results from running the R code and compared it to the "pure case".

Bottom lines:

What this tells me is:

  1. I got pretty consistent results using two separate modes so if I made a mistake, then I did it two separate times.
  2. If I have a critter that is both treated AND diseased, it is 10x more likely to be from Species30 than to be from Species10.
  3. If I have a critter that is diseased, the chance of it having been treated is 2.75x higher if it is from Species30 than if it is from Species10.
  4. The sample size may be high, but given the 1% rate for species 10, there is still a pretty big difference between realistic results and the limit of infinite samples. That 12% relative error when all others are between 5% and 0.2% is a warning. You might want to unset the seed from the code above and run the simulation a few hundred times, storing the estimates, and determine not only the mean value, but the tolerance around the mean given your sample sizes. This will put some error bars around the 10x and the 2.75x.

Having a clear understanding of the difference between points 2 and 3 is likely going to be important in the discussion. The difference is whether or not you know if it has been treated. For the intended outcome do you want knowledge of treatment to be a given as in you are going to treat the entire population? Are you looking at narrower scope in your solution? Every technical analysis has to support a business decision. Every business decision is about maximizing return of value. Reliable and excellent growth come from maximizing return of value to the customer. This means that many analyses could be improved in clarity and effect if the nature of the customer, the nature of value, and the nature of the decision are clearly articulated. It is better to ask the question 2 or 20 times and be sure you are asking the right question than to spend resources 2 or 20 times answering the wrong question. Most business harm, like many aircraft accidents, are a tragic sequence of repeatedly asking the wrong question.

Good management of point 4 is likely going to ensure repeatability. My estimate is a mean absolute relative error of 14.45% but a max near 92%. You are not unlikely to get an experimental rate 14.67% through 5.52% instead of 9% for the "probability of treated given species10 and diseased". If N is 4000, then the range is from 13.02 to 6.22. Doubling the sample size shaves of around 1% from either side of the tolerance band.

To be completed later. Lunch is calling.

$\endgroup$
  • $\begingroup$ Thanks for the great explanation. I am not too worried about the confidence intervals because my mice populations contain a few thousand mice in each. My question is more along the line of figuring out if the relative fold reduction of disease incidence after treatment is different for my two mice strain (as EdM noted above). $\endgroup$ – Ehi Akhirome Jun 4 '15 at 4:48
  • $\begingroup$ Then you are more interested in part 2. Is it okay if I demo my results in "R"? $\endgroup$ – EngrStudent - Reinstate Monica Jun 4 '15 at 13:49
  • 1
    $\begingroup$ Do you have time to complete now? $\endgroup$ – kjetil b halvorsen Jul 23 at 10:36
  • 1
    $\begingroup$ @kjetilbhalvorsen - I have since lost my train of thought. It was about 4 years ago. $\endgroup$ – EngrStudent - Reinstate Monica Jul 25 at 0:29

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.