Haciendo una muestra de 3 χ2 prueba ( prop.test en R) la variabilidad de cada madre no se tiene en cuenta.
Sí. El χ2 requiere un muestreo independiente para cada célula.
Con un Wilcoxon-Mann-Whitney no estaría probando realmente proporciones, ¿verdad?
Estarías probando las proporciones, pero después de que estén clasificadas. Además, esta prueba es una prueba de dos muestras, por lo que tendrías que realizar un trabajo adicional.
Con estos datos, se empieza a salir del mundo de las pruebas estadísticas.
Si la proporción de hombres+mujeres no es demasiado baja o alta, se puede modelar de efecto mixto lineal generalizado como punto de partida. punto de partida.
Podemos configurar algunos datos falsos para ver. En primer lugar, vamos a conseguir algunas herramientas.
# Load handy packages.
library(ggplot2)
library(plyr)
library(lme4)
library(lsmeans)
library(car)
library(multcomp)
A continuación, establezca la estructura de los datos.
# 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))
)
Ahora calcula las proporciones y arregla el conjunto de datos para que el análisis de la varianza funcione correctamente.
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)
Sólo para ver si tiene sentido:
# 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
Ahora, analiza los datos. Una forma de hacerlo es tener una respuesta binomial utilizando el enlace logit de la función logit. Esta es una implementación que incluye el Día como efecto fijo como efecto fijo, la madre como efecto aleatorio y la interacción madre-día como efecto aleatorio.
# 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
Hay un efecto Día. Una forma de hacerlo es mirar las comparaciones por pares.
# 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"
Sin embargo, en su caso especificó que estaba interesado en comparar el día 1 con el día 2 y el día 3 --- presumiblemente, la media de los dos es similar a lo que quería decir.
# 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)
Una fuente de confusión con las estimaciones y las medias de mínimos cuadrados es la escala en la que se realiza el análisis.
# 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
Todos ellos están en la escala logística. Podemos transformarlos de nuevo.
lsm <- c(0.4895782, 0.2435273, 0.3192616)
exp(lsm) / (1 + exp(lsm))
#[1] 0.6200071 0.5605827 0.5791443
Bastante buena coincidencia con las medias brutas. Eso no tiene por qué ocurrir.
También puede utilizar el nlme y modelar el R -correlación lateral estructura, quizás con un proceso AR(1) como ejemplo.
