# Additional checks documented in protocol/deviations.md; no confirmatory claims.
options(digits=15)
source('scripts/analyze.R')
z<-pub[pub$pubmed_id!=101039,];post<-summarize_pub(z$delta,'exclude_mismatched_source_101039')
for(label in c('Yes','No')){z<-v[v$species_clean=='mus musculus'&v$ITP==label,];u<-aggregate(delta~pubmed_id,z,mean);post<-rbind(post,summarize_pub(u$delta,paste0('mouse_ITP_',label)))}
write_csv(post,'results/drugage_posthoc.csv')
set.seed(20260905);famerr<-aggregate(cbind(sqerr_full,sqerr_base)~family,cv,sum);famerr$n<-as.integer(table(cv$family)[famerr$family]);bn<-nrow(famerr)
bdelta<-replicate(10000,{i<-sample.int(bn,bn,replace=TRUE);sum(famerr$sqerr_full[i]-famerr$sqerr_base[i])/sum(famerr$n[i])})
write_csv(data.frame(estimate=mean(cv$sqerr_full-cv$sqerr_base),lo=quantile(bdelta,.025),hi=quantile(bdelta,.975)),'results/anage_cv_bootstrap.csv')
# Student tail cross-check via direct integration, independently of pt().
checktail<-function(t,df)2*integrate(function(x) exp(lgamma((df+1)/2)-lgamma(df/2)-.5*log(df*pi)-(df+1)/2*log1p(x*x/df)),abs(t),Inf,rel.tol=1e-10,abs.tol=1e-20)$value
pa<-checktail(models$primary$beta/models$primary$se,models$primary$clusters-1);pd<-checktail(dr$estimate/dr$se,dr$units-1)
stopifnot(abs(pa/models$primary$p-1)<1e-8,abs(pd/dr$p-1)<1e-8)
write_csv(data.frame(paper=c('anage','drugage'),integrated_p=c(pa,pd)),'results/pvalue_integration_check.csv')
print(post);print(range(loo$beta));print(range(leave$estimate))
