5
$\begingroup$

Here is my question.

20 (worm) mothers produce males + females and hermafrodites (many every hour).

I want to test if the proportion of males + females born on day #1 is the same as in days #2 and #3.

I believe the data of the three days is dependent, because the mothers are the same.

Doing a 3-sample $\chi^2$ test (prop.test function in R) the variability of each mother is not taken into account.

With a Wilcoxon-Mann-Whitney I would not be actually testing proportions, right?

Which test would be suitable?

Any suggestions would be highly appreciated.

$\endgroup$
1
$\begingroup$

Doing a 3-sample $\chi^2$ test (prop.test function in R) the variability of each mother is not taken into account.

Right. The $\chi^2$ test requires independent sampling for each cell.

With a Wilcoxon-Mann-Whitney I would not be actually testing proportions, right?

You would be testing proportions, but after they are ranked. Also, this test is a two-sample test, so you would have to go through some additional work.

With these data, you are starting to leave the world of specific statistical tests.

If the proportion of male+female is not too low or high, you can model this using a generalized linear mixed effect model as a starting point.

We can set up some fake data to see. First, we'll get some tools.

# Load handy packages.

library(ggplot2)
library(plyr)
library(lme4)
library(lsmeans)
library(car)
library(multcomp)

Then, set up the structure of the data.

# Generate some fake data.  Let's say that
# the proportion of male+female is higher on
# day #1 than on days #2 and #3.  Let's also
# say that there is variation among mothers.

set.seed(875487)

K <- 20     ### the number of mothers

# Specify the mean number of male+female per day.

mf.day1 <- rnorm(K, mean=33, sd=3)
mf.day2 <- rnorm(K, mean=25, sd=3)
mf.day3 <- rnorm(K, mean=24, sd=3)

# Specify the mean number of hermaphrodites per day.

h.day1 <- rnorm(K, mean=20, sd=3)
h.day2 <- rnorm(K, mean=20, sd=3)
h.day3 <- rnorm(K, mean=20, sd=3)

# Generate the data for three days.

Data <- data.frame(
    Mother = rep(1:K, 3),
    Day = rep(1:3, each=K),
    MF = c(rpois(K, mf.day1), rpois(K, mf.day2), rpois(K, mf.day3)),
    H  = c(rpois(K,  h.day1), rpois(K,  h.day2), rpois(K,  h.day3))
)

Now calculate proportions and fix up the data set so the analysis of variance will work correctly.

Data$p <- Data$MF / ( Data$MF + Data$H )

# Make factors.  Be careful when making factors of integers.

Data$Mother <- factor(Data$Mother)
Data$Day <- factor(Data$Day)

Just to see if it makes sense:

# Look at the data.

ggplot(Data, aes(x=Day, y=p, group=Mother)) +
    geom_line()

ddply(Data, .(Day), summarise, Mean=mean(p))

#  Day      Mean
#1   1 0.6223710
#2   2 0.5607464
#3   3 0.5853103

Now, analyze the data. One way to do it is to have a binomial response using the logit link function. Here is an implementation that includes Day as a fixed effect, Mother as a random effect, and Mother by Day interaction as a random effect.

# Analyze the data with a random effect model.

fit <- glmer(cbind(MF, H) ~ Day + (1|Mother) + (1|Mother:Day),
             data=Data, family=binomial)

Anova(fit)

#Analysis of Deviance Table (Type II Wald chisquare tests)
#
#Response: cbind(MF, H)
#     Chisq Df Pr(>Chisq)  
#Day 6.9907  2    0.03034 *
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a Day effect. One way to go is to look at pairwise comparisons.

# Perform pairwise comparisons.

tuk <- glht(fit, linfct = mcp(Day = "Tukey"))
tuk

#    General Linear Hypotheses
#
#Multiple Comparisons of Means: Tukey Contrasts
#
#
#Linear Hypotheses:
#           Estimate
#2 - 1 == 0 -0.24605
#3 - 1 == 0 -0.17032
#3 - 2 == 0  0.07573

cld(tuk)

#   1    2    3
# "b"  "a" "ab"

However, in your case you specified that you were interested in comparing day #1 versus the day #2 and day #3 --- presumably, the average of the two is similar to what you meant.

# Test the contrast of interest, day #1 - (day #2 + day #3)/2.
# The contrast coefficients look a little weird, but it's due to the
# coding used.

con <- glht(fit, linfct = matrix(c(0, -0.5, -0.5), 1))
summary(con)

#    Simultaneous Tests for General Linear Hypotheses
#
#Fit: glmer(formula = cbind(MF, H) ~ Day + (1 | Mother) + (1 | Mother:Day),
#    data = Data, family = binomial)
#
#Linear Hypotheses:
#       Estimate Std. Error z value Pr(>|z|)
#1 == 0   0.2082     0.0823    2.53   0.0114 *
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#(Adjusted p values reported -- single-step method)

One source of confusion with the estimates and least squares means is the scale on which the analysis is performed.

# Look at the estimated marginal means.

lsmeans(fit, "Day")

# Day    lsmean         SE df  asymp.LCL asymp.UCL
# 1   0.4895782 0.07094520 NA 0.35051135 0.6286451
# 2   0.2435273 0.07391394 NA 0.09864109 0.3884135
# 3   0.3192616 0.07409900 NA 0.17401260 0.4645105
#
#Confidence level used: 0.95

These are all on the logistic scale. We can transform them back.

lsm <- c(0.4895782, 0.2435273, 0.3192616)

exp(lsm) / (1 + exp(lsm))

#[1] 0.6200071 0.5605827 0.5791443

Pretty good agreement with the raw means. That does not have to happen.

You could also use the nlme package and model the $R$-side correlation structure, perhaps with an AR(1) process as an example.

$\endgroup$
  • 1
    $\begingroup$ On the last part with the lsmeans result. You can also back-transform them using summary(lsmeans(fit, "Day"), type = "response") $\endgroup$ – rvl Jun 6 '15 at 13:48
  • $\begingroup$ Thank you! I am digesting the answer. I will accept it as soon as that is done. $\endgroup$ – Tássio Jun 10 '15 at 8:15

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.