# Use "all" version since some haplos are missing from the "drop-1" version
hapcount <- read.table( 'haplo-freq-all.txt', header=T, row=NULL, stringsAsFactors=FALSE)
hapfreq <- hapcount$n / sum( hapcount$n)
names( hapfreq) <- hapcount$H
haps <- read.table( 'xc-hs-river-haplo.txt', header=T, row=NULL, stringsAsFactors=FALSE)
haps <- within( haps, {
h1 <- 'H' %&% h1
h2 <- 'H' %&% h2
Pr_if_mat <- ifelse( h1==h2, hapfreq[ h1], 0)
Pr_if_pat <- hapfreq[ h1] * hapfreq[ h2]
})
# Log-likelihood as a function of ppn of HSP that are Maternal
# lglk2 has different ppns for same-river and cross-river; params on logit scale
lglk2 <- function( logit_Prmat)
sum( log( inv.logit( logit_Prmat)[ samediff] * (Pr_if_mat - Pr_if_pat) + Pr_if_pat))
environment( lglk2) <- list2env( haps) # so it finds Pr_if_mat/pat
orig_samediff <- haps$samediff # for debugging
# lglk1 enforces same ppn regardless; *not* on logit scale, for dull technical reasons
lglk1 <- function( Prmat1) { lp <- logit( Prmat1); lglk2( c( S=lp, D=lp)) }
# If no "river effect":
fit1 <- optimize( NEG( lglk1), lower=0, upper=1)
Prmatsame <- fit1$minimum
# Allowing "river effect":
lp <- logit( Prmatsame)
fit2 <- optim( c( S=lp, D=lp), NEG( lglk2), method='BFGS')
Delta <- -2*(fit2$value - fit1$objective)
# GLRT test--- but sodding boundary sodding value...
# NB NEG() in optimizer calls, hence multiplier is "-2" here
# signif <- pchisq( Delta, df=1, lower=FALSE)
# Null distro via perm test
nsame <- sum( haps$samediff=='S')
npair <- nrow( haps)
Dnull <- function( iter) {
if( iter %% 10 == 0) {
cat( '\r', iter); flush.console()
}
which.same <- rsample( nsame, 1:npair, replace=FALSE)
diffo <- rep( 'D', npair)
diffo[ which.same] <- 'S'
assign( 'samediff', diffo, environment( lglk2))
sfit1 <- optimize( NEG( lglk1), lower=0, upper=1)
sfit2 <- optim( c( S=lp, D=lp), NEG( lglk2), method='BFGS')
Dsim <- -2*(sfit2$value - sfit1$objective)
return( Dsim)
}
null_samps <- do.on( 1:10000, Dnull(.))
cat( sprintf( 'Observed lglk ratio Delta: %5.3f', Delta), '\n')
cat( 'MLE(ppn mat) under H1:\n')
print( fit2$par)
cat( sprintf( 'Ppn of null-distro Deltas above observed value: %5.3f (n=%i simulations)',
mean( null_samps >= Delta), length( null_samps)), '\n')