Estou trabalhando em um problema no qual tenho vários pares de homens atualmente vivos, i
cada um com um suposto ancestral paterno ni
gerações atrás (com base em evidências genealógicas) e onde tenho informações sobre se há uma incompatibilidade no genótipo cromossômico Y (exclusivamente paternalmente herdadoxi
= 1 para incompatibilidade, 0 se houver uma correspondência). Se não há incompatibilidade, eles realmente têm um ancestral paterno comum, mas se existe, deve ter havido uma torção na cadeia como resultado de um ou mais casos extraconjugais (só posso detectar se nenhum ou pelo menos um desses eventos de paternidade extra-par aconteceu, ou seja, a variável dependente é censurada). O que me interessa é obter uma estimativa de máxima verossimilhança (mais limites de confiança de 95%) não apenas da taxa média de paternidade de par extra (EPP) (probabilidade de que por geração uma criança derivaria de uma cópula de par extra), mas também para tentar inferir como a taxa de paternidade dos pares extras pode ter mudado em função do tempo (como o número de gerações que separaram o ancestral paterno comum deve ter algumas informações sobre isso - quando há uma incompatibilidade, não Não sei, porém, quando os EPPs teriam acontecido, como poderia estar em qualquer lugar entre a geração desse ancestral presumido e o presente, mas quando houver uma correspondência, temos certeza de que não houve EPPs em nenhuma das gerações anteriores). Portanto, tanto a minha variável binomial dependente quanto a geração / tempo covariável independente são censuradas. Com base em um problema semelhante, publicadoaqui, eu já descobri como poderia fazer uma estimativa de probabilidade máxima da população e a taxa de paternidade extra-média de tempo médio mais os phat
intervalos de confiança de probabilidade de perfil de 95% em R da seguinte maneira:
# 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 /stats/152111/censored-binomial-model-log-likelihood
require(bbmle)
fit = mle2(neglogL, start=list(p=0.01), data=list(x=x, n=n))
return(fit)
}
Exemplo com alguns dados piloto (de Larmuseau et al. ProcB 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
Estimativa de máxima verossimilhança da população e taxa de paternidade extra-média de tempo médio mais limites de confiança de 95%:
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
isto é, 0,9% [0,43-1,75% 95% CLs] de todas as crianças foram derivadas de um pai diferente do suposto.
Eu queria dar um passo adiante e também tentar estimar uma possível tendência temporal na taxa de paternidade dos pares extras em p
função da geração ni
(por simplicidade, assumindo uma relação linear entre as chances de log de observar um evento de paternidade dos pares extras e geração), levando em consideração que, se ocorrer uma incompatibilidade, os eventos EPP poderiam ter ocorrido em qualquer lugar entre a geração do ancestral comum ni
e o presente (geração 0), e que, se não houvesse incompatibilidade, nenhum evento EPP poderia ter ocorrido em algum momento. qualquer das gerações anteriores para esse par específico de indivíduos.
Se antes que assumiu a probabilidade de uma criança a ser derivado a partir de um extra-par cópula para ser constante, e se foi uma variável aleatória ser igual a quando foi observado um cromossoma incompatibilidade Y (que corresponde a um ou mais eventos de EPP) e caso contrário, a probabilidade de não haver incompatibilidade (isto é, ) quando o ancestral paterno viveu gerações atrás ( ) era , enquanto a chance de observar um evento EPP foi
Em um conjunto de dados de observações independentes de com ancestrais paternos vivendo gerações atrás, a probabilidade era
resultando em uma probabilidade logarítmica de
Considerando que no meu modelo mais complexo que incorpora a dinâmica temporal, quero que seja uma função de agora, com , ou seja,
Alterei a definição da função de probabilidade acima e a maximizei usando a função mle2
do pacote 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 /stats/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
Isso me dá uma estimativa histórica razoável da evolução da taxa de PPE longo do tempo:
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)")
No entanto, ainda estou um pouco preocupado em como calcular os intervalos de confiança de 95% na previsão geral desse modelo. Alguém saberia fazer isso por acaso? Talvez usando intervalos de previsão da população (por reamostragem de parâmetros de acordo com o ajuste após uma distribuição normal multivariada) (ou o método delta também funcionaria?)? E alguém poderia comentar se minha lógica acima está correta? Também estava me perguntando se esse tipo de modelo binomial censurado é conhecido sob algum nome padrão na literatura e se alguém conhece algum trabalho publicado sobre esse tipo de cálculo de ML nesse tipo de modelo? (Tenho a sensação de que o problema deve ser razoavelmente padrão e corresponder a algo que já foi feito, mas parece que não consigo encontrar nada ...)
[Documentos do PS com mais informações sobre esse tópico / problema estão disponíveis aqui e aqui]
fonte
Respostas:
Com base no capítulo de Ben Bolker e no comentário acima, concluí que os intervalos de previsão de 95% da população são dados por
fonte