Estoy trabajando en un problema en el que tengo varias parejas de la actualidad que viven los varones i de que cada uno tienen un presunto antepasado paterno ni generaciones atrás (basado en genealógico de la evidencia) y donde tengo información sobre si existe un desajuste en su Y cromosómicas genotipo (exclusivamente paterna heredada, xi = 1 para el desajuste, 0 si no hay un partido). Si no hay ninguna coincidencia, la verdad ellos tienen en común un antepasado paterno, pero si hay un debe haber sido una torcedura en la cadena como resultado de uno o más asuntos extramaritales (yo sólo puede detectar, aunque si ninguno o al menos uno de esos paternidad extra-pareja, evento que sucedió, es decir, la variable dependiente es censurado). Lo que me interesa es la obtención de una estimación de máxima verosimilitud (más del 95% límites de confianza), y no sólo de la media de la paternidad extra-pareja (PPE) tasa (la probabilidad de que por la generación de un niño podría ser derivada de un extra-pareja de la cópula), pero también para tratar de inferir cómo la paternidad extra-pareja, la tasa puede haber cambiado como una función del tiempo (como el número de generaciones que separan el común antepasado paterno debe tener algo de información sobre esto - cuando hay un desajuste yo no sé, aunque cuando la EPPs hubiera ocurrido, como podría haber sido en cualquier lugar entre la generación de ese presunto progenitor y el presente, pero cuando hay un partido, estamos seguros de que no había EPPs en ninguna de las generaciones anteriores). Por lo tanto, tanto para mi dependiente de la binomial variable independiente y covariable generación/tiempo están censurados. Basado en un poco de problema similar publicado aquí ya he averiguado cómo yo podría hacer una máxima probabilidad de la estimación de la población y el tiempo promedio de paternidad extra-pareja tasa phat más del 95% perfil de probabilidad de los intervalos de confianza en R como sigue :
# Function to make overall ML estimate of EPP rate p plus 95% profile likelihood confidence intervals,
# taking into account that for pairs with mismatches multiple EPP events could have occured
#
# input is
# x=vector of booleans or 0 and 1s specifying if there was a mismatch or not (1 or TRUE = mismatch)
# n=vector with nr of generations that separate common ancestor
# output is mle2 maximum likelihood fit with best-fit param EPP rate p
estimateP = function(x, n) {
if (!is.logical(x[[1]])) x = (x==1)
neglogL = function(p, x, n) -sum((log(1 - (1-p)^n))[x]) -sum((n*log(1-p))[!x]) # negative log likelihood, see see http://stats.stackexchange.com/questions/152111/censored-binomial-model-log-likelihood
require(bbmle)
fit = mle2(neglogL, start=list(p=0.01), data=list(x=x, n=n))
return(fit)
}
Ejemplo, con algunos datos experimentales (de Larmuseau et al. ProcB de 2010):
n = c(7, 7, 7, 7, 7, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 15, 15, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 18, 18, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 22, 23, 23, 24, 24, 25, 27, 31) # number of generations/meioses that separate presumed common paternal ancestor
x = c(rep(0,6), 1, rep(0,7), 1, 1, 1, 0, 1, rep(0,20), 1, rep(0,13), 1, 1, rep(0,5)) # whether pair of individuals had non-matching Y chromosomal genotypes
El de máxima probabilidad de la estimación de la población y el tiempo promedio de paternidad extra-pareja, la tasa de más del 95% límites de confianza :
fit = estimateP(x,n)
c(coef(fit),confint(fit))*100 # estimated p and profile likelihood confidence intervals
# p 2.5 % 97.5 %
# 0.9415172 0.4306652 1.7458847
es decir, el 0,9% [0.43-1.75% 95% C. L. s] de todos los niños fueron derivados de un padre que era diferente a la de la supuesta uno.
Yo quería ir un paso más allá, y también tratar de estimar una posible tendencia temporal en la paternidad extra-pareja tasa p como una función de generación de ni (para simplificar, suponiendo una relación lineal entre las probabilidades de registro de observación de una paternidad extra-pareja de eventos y generación), teniendo en cuenta que si se produce una discrepancia el PPE eventos podría haber tenido lugar en cualquier lugar entre la generación de un ancestro común ni y el presente (generación 0), y que si no hay un desajuste que no PPE evento podría tener lugar en cualquiera de las generaciones anteriores para cada par de individuos.
Si antes de que asumimos que la probabilidad de que un niño sea derivado de un extra-pareja, la cópula $p$ a ser constante, y si $X$ es una variable aleatoria igual a $1$ cuando un cromosoma Y desajuste se observó (correspondiente a 1 o más PPE eventos) y $0$ lo contrario, entonces la probabilidad de observar sin falta de coincidencia (es decir, $X=0$) cuando la paternal antepasado vivió $n$ generaciones atrás ($n = 1, 2, 3, \ldots$ )$(1-p)^n$, mientras que la probabilidad de observar un PPE evento fue
$$\Pr(X=1\,|\, n) = 1 - (1-p)^n.$$
In a dataset of independent observations $\mathbf{x} = (x_1, x_2, \ldots)$ of $X$ with paternal ancestors living $\mathbf{n} = (n_1, n_2, \ldots)$ generations ago, the likelihood therefore was
$$L(p; \mathbf{x}, \mathbf{n}) = \prod_{x_i=1} \left(1 - (1-p)^{n_i}\right)\prod_{x_j=0} (1-p)^{n_j},$$
resulting in a log likelihood of
$$\Lambda(p) = \sum_{x_i=1} \log\left(1 - (1-p)^{n_i}\right) + \sum_{x_j=0} {n_j} \log (1-p).$$
Taking into account that in my more complex model incorporating temporal dynamics I want $p$ to be a function of $$ n ahora, con $$log (p/(1-p)) = a + b.n$$ , es decir, $$p(a,b,n) = \exp(a+b.n) / (1+\exp(a+b.n))$$
I then changed the definition of the likelihood function above accordingly and maximized it using function mle2 from package bbmle :
# ML estimation, assuming that EPP rate p shows a temporal trend
# where log(p/(1-p))=a+b*n
# input is
# x=vector of booleans or 0 and 1s specifying if there was a mismatch or not (1 or TRUE = mismatch)
# n=vector with nr of generations that separate common ancestor
# output is mle2 maximum likelihood fit with best-fit params a and b
estimatePtemp = function(x, n) {
if (!is.logical(x[[1]])) x = (x==1)
pfun = function(a, b, n) exp(a+b*n)/(1+exp(a+b*n)) # we now write p as a function of a, b and n
logL = function(a, b, x, n) sum((log(1 - (1-pfun(a, b, n))^n))[x]) +
sum((n*log(1-pfun(a, b, n)))[!x]) # a and b are params to be estimated, modified from http://stats.stackexchange.com/questions/152111/censored-binomial-model-log-likelihood
neglogL = function(a, b, x, n) -logL(a, b, x, n) # negative log-likelihood
require(bbmle)
fit = mle2(neglogL, start=list(a=-3, b=-0.1), data=list(x=x, n=n))
return(fit)
}
# fitted coefficients
estfit = estimatePtemp(x, n)
cbind(coef(estfit),confint(estfit)) # parameter estimates and profile likelihood confidence intervals
# 2.5 % 97.5 %
# a -3.09054167 -5.3191406 -1.12078519
# b -0.09870851 -0.2396262 0.02848305
summary(estfit)
# Coefficients:
# Estimate Std. Error z value Pr(z)
# a -3.090542 1.057382 -2.9228 0.003469 **
# b -0.098709 0.067361 -1.4654 0.142819
This gives me a reasonable looking historical estimate of the evolution of the EPP rate $p$ a lo largo del tiempo :
pfun = function(a, b, n) exp(a+b*n)/(1+exp(a+b*n))
xvals=1:max(n)
par(mfrow=c(1,1))
plot(xvals,sapply(xvals,function (n) pfun(coef(estfit)[1], coef(estfit)[2], n)),
type="l", xlab="Generations ago (n)", ylab="EPP rate (p)")
Sin embargo, todavía estoy un poco atascado en la forma de calcular el 95% de intervalos de confianza en la predicción general de este modelo. Podría alguien sabe cómo hacer que por casualidad? Tal vez el uso de la población intervalos de predicción (por remuestreo de los parámetros según el ajuste siguiendo una distribución normal multivariante) (o el método delta también?)? Y alguien podría comentar si mi lógica anterior es correcto? También me pregunto si este tipo de censurados modelo binomial es conocido bajo algún nombre estándar en la literatura, y si alguien sabe de alguna obra publicada en hacer este tipo de ML cálculos en este tipo de modelo? (Tengo la sensación de que el problema debe ser bastante estándar y corresponden a algo que ha sido hecho ya, pero parece que no puede encontrar cualquier cosa...)
[PS Papeles con más antecedentes sobre este tema/problema están disponibles aquí y aquí]
