# Reproduce both analyses with base R. Run from repository root: Rscript scripts/analyze.R
options(stringsAsFactors=FALSE, digits=15)
set.seed(20260905)
dir.create('results',showWarnings=FALSE);dir.create('data/derived',showWarnings=FALSE)
write_csv <- function(x,f) write.csv(x,f,row.names=FALSE,na='')
# Cluster sandwich CR1: G/(G-1)*(N-1)/(N-K), t reference G-1.
cluster_result <- function(fit, groups, term='lmaturity') {
 X<-model.matrix(fit); u<-residuals(fit); n<-nrow(X);k<-ncol(X)
 stopifnot(fit$rank==k, length(groups)==n)
 sc<-rowsum(X*as.numeric(u),groups,reorder=FALSE);g<-nrow(sc)
 bread<-solve(crossprod(X));v<-bread%*%crossprod(sc)%*%bread*g/(g-1)*(n-1)/(n-k)
 b<-coef(fit)[term];se<-sqrt(v[term,term]);ci<-b+c(-1,1)*qt(.975,g-1)*se
 data.frame(n=n,clusters=g,beta=unname(b),se=se,lo=ci[1],hi=ci[2],p=2*pt(-abs(b/se),g-1),percent_doubling=100*(2^b-1),percent_lo=100*(2^ci[1]-1),percent_hi=100*(2^ci[2]-1))
}
# AnAge: preserve raw source row for traceability.
a<-read.delim('data/raw/anage_anage_data.txt',check.names=FALSE,quote='',na.strings=c('','NA'))
a$source_row<-seq_len(nrow(a)); stopifnot(!anyDuplicated(a$HAGRID))
a$species_name<-paste(a$Genus,a$Species)
m<-a[a$Class=='Mammalia',]
m$eligible<-with(m,is.finite(`Female maturity (days)`)&`Female maturity (days)`>0 & is.finite(`Adult weight (g)`)&`Adult weight (g)`>0 & is.finite(`Maximum longevity (yrs)`)&`Maximum longevity (yrs)`>0 & !is.na(Order)&!is.na(Family)&`Data quality`%in%c('acceptable','high'))
p<-m[m$eligible,]
p$lmaturity<-log(p$`Female maturity (days)`);p$lmass<-log(p$`Adult weight (g)`);p$llife<-log(p$`Maximum longevity (yrs)`)
p$order<-p$Order;p$family<-p$Family;p$origin<-p$`Specimen origin`;p$effort<-p$`Sample size`
p$flag_maturity_exceeds_life<-p$`Female maturity (days)`/365.25 >= p$`Maximum longevity (yrs)`
qc_a<-data.frame(metric=c('raw_export','mammals','acceptable_high_mammals','positive_maturity_mammals','positive_mass_mammals','positive_longevity_mammals','eligible','families','orders','maturity_exceeds_life'),value=c(nrow(a),nrow(m),sum(m$`Data quality`%in%c('acceptable','high')),sum(m$`Female maturity (days)`>0,na.rm=TRUE),sum(m$`Adult weight (g)`>0,na.rm=TRUE),sum(m$`Maximum longevity (yrs)`>0,na.rm=TRUE),nrow(p),length(unique(p$family)),length(unique(p$order)),sum(p$flag_maturity_exceeds_life)))
write_csv(qc_a,'results/anage_qc.csv');write_csv(m,'data/derived/anage_mammals_screened.csv')
f<-lm(llife~lmaturity+lmass+factor(order),p)
p$cooks_distance<-cooks.distance(f);p$leverage<-hatvalues(f)
write_csv(p,'data/derived/anage_included.csv')
write_csv(data.frame(model.matrix(f),llife=p$llife,family=p$family,check.names=FALSE),'data/derived/anage_design.csv')
models<-list(primary=cluster_result(f,p$family),unadjusted=cluster_result(lm(llife~lmaturity,p),p$family),mass_adjusted=cluster_result(lm(llife~lmaturity+lmass,p),p$family),family_intercepts=cluster_result(lm(llife~lmaturity+lmass+factor(family),p),p$family))
for (name in c('no_humans','no_bats_primates','high_quality','captivity','no_flagged')) {
 z<-switch(name,no_humans=p[p$species_name!='Homo sapiens',],no_bats_primates=p[!p$order%in%c('Chiroptera','Primates'),],high_quality=p[p$`Data quality`=='high',],captivity=p[p$origin=='captivity',],no_flagged=p[!p$flag_maturity_exceeds_life,])
 if(nrow(z)>30 && length(unique(z$family))>10) models[[name]]<-cluster_result(lm(llife~lmaturity+lmass+factor(order),z),z$family)
}
models$effort_origin_adjusted<-cluster_result(lm(llife~lmaturity+lmass+factor(order)+factor(effort)+factor(origin),p),p$family)
ma<-do.call(rbind,lapply(names(models),function(nm)cbind(model=nm,models[[nm]])));rownames(ma)<-NULL;write_csv(ma,'results/anage_models.csv')
loo<-do.call(rbind,lapply(sort(unique(p$order)),function(o){z<-p[p$order!=o,]; cbind(excluded=o,cluster_result(lm(llife~lmaturity+lmass+factor(order),z),z$family))}));write_csv(loo,'results/anage_leave_order_out.csv')
write_csv(p[order(-p$cooks_distance),c('species_name','source_row','cooks_distance','leverage')],'results/anage_influence.csv')
# Shared deterministic family-blocked folds; no unseen order in evaluation.
nfam<-tapply(p$family,p$order,function(x)length(unique(x)));cv<-p[p$order%in%names(nfam[nfam>=2]),];cv$fold<-NA_integer_
set.seed(20260905)
for(o in sort(unique(cv$order))){fams<-sample(sort(unique(cv$family[cv$order==o])));assign<-setNames(rep(1:5,length.out=length(fams)),fams);cv$fold[cv$order==o]<-assign[cv$family[cv$order==o]]}
cv$pred_base<-cv$pred_full<-NA_real_
for(j in 1:5){tr<-cv[cv$fold!=j,];te<-cv[cv$fold==j,];if(nrow(te)==0)next
 stopifnot(all(unique(te$order)%in%tr$order),!any(unique(te$family)%in%tr$family))
 cv$pred_base[cv$fold==j]<-predict(lm(llife~lmass+factor(order),tr),te)
 cv$pred_full[cv$fold==j]<-predict(lm(llife~lmaturity+lmass+factor(order),tr),te)
}
cv$sqerr_base<-(cv$llife-cv$pred_base)^2;cv$sqerr_full<-(cv$llife-cv$pred_full)^2
write_csv(cv,'data/derived/anage_cv_predictions.csv')
cvres<-data.frame(n=nrow(cv),families=length(unique(cv$family)),orders=length(unique(cv$order)),rmse_base=sqrt(mean(cv$sqerr_base)),rmse_full=sqrt(mean(cv$sqerr_full)),relative_rmse_reduction=1-sqrt(mean(cv$sqerr_full)/mean(cv$sqerr_base)),mean_paired_sqerr_difference=mean(cv$sqerr_full-cv$sqerr_base))
write_csv(cvres,'results/anage_cv.csv')
# DrugAge: deduplicate ORIGINAL rows, preserve source row after filtering.
d<-read.csv('data/raw/drugage_drugage.csv',check.names=FALSE,na.strings=c('','NA'));dup<-duplicated(d);d$source_row<-seq_len(nrow(d));raw_n<-nrow(d);d<-d[!dup,]
d$species_clean<-tolower(trimws(d$species));d$compound_clean<-tolower(trimws(d$compound_name))
d$paired<-is.finite(d$avg_lifespan_change_percent)&is.finite(d$max_lifespan_change_percent)&!is.na(d$pubmed_id)
v<-d[d$paired,];v$delta<-v$max_lifespan_change_percent-v$avg_lifespan_change_percent
v$flag_nonphysical<-v$avg_lifespan_change_percent<= -100 | v$max_lifespan_change_percent<= -100
write_csv(d,'data/derived/drugage_screened.csv');write_csv(v,'data/derived/drugage_paired.csv')
pub<-aggregate(cbind(delta,avg_lifespan_change_percent,max_lifespan_change_percent)~pubmed_id,v,mean);pub$n_rows<-as.integer(table(v$pubmed_id)[as.character(pub$pubmed_id)])
stopifnot(nrow(pub)>=20)
write_csv(pub,'data/derived/drugage_publication_means.csv')
summarize_pub<-function(x,label,seed=20260905){set.seed(seed);n<-length(x);b<-replicate(10000,mean(sample(x,n,replace=TRUE)));tt<-t.test(x);data.frame(analysis=label,units=n,estimate=mean(x),boot_lo=quantile(b,.025,names=FALSE),boot_hi=quantile(b,.975,names=FALSE),t_lo=tt$conf.int[1],t_hi=tt$conf.int[2],se=sd(x)/sqrt(n),p=tt$p.value,median=median(x),negative=sum(x<0),zero=sum(x==0),positive=sum(x>0))}
dr<-summarize_pub(pub$delta,'publication_balanced')
# Sensitivity estimates use different weighting, not independent confirmations.
compound<-aggregate(delta~compound_clean,v,mean);rr<-v[!v$flag_nonphysical,];rr$logcontrast<-log1p(rr$max_lifespan_change_percent/100)-log1p(rr$avg_lifespan_change_percent/100);rrpub<-aggregate(logcontrast~pubmed_id,rr,mean)
sens<-rbind(dr,summarize_pub(compound$delta,'compound_balanced'),summarize_pub(rrpub$logcontrast,'log_response_ratio'))
write_csv(sens,'results/drugage_estimates.csv')
sp<-do.call(rbind,lapply(c('caenorhabditis elegans','drosophila melanogaster','mus musculus'),function(s){z<-v[v$species_clean==s,];u<-aggregate(delta~pubmed_id,z,mean);cbind(species=s,paired_rows=nrow(z),summarize_pub(u$delta,s))}));sp$p_holm<-p.adjust(sp$p,'holm');write_csv(sp,'results/drugage_species.csv')
leave<-data.frame(pubmed_id=pub$pubmed_id,estimate=(sum(pub$delta)-pub$delta)/(nrow(pub)-1));write_csv(leave,'results/drugage_leave_publication_out.csv')
qcd<-data.frame(metric=c('raw_rows','exact_duplicates_removed','deduplicated_rows','paired_rows','paired_publications','paired_compounds','paired_species','missing_average','missing_maximum','missing_pubmed','nonphysical_paired','row_weighted_delta','publication_mean_average','publication_mean_maximum','publication_median_delta'),value=c(raw_n,sum(dup),nrow(d),nrow(v),nrow(pub),length(unique(v$compound_clean)),length(unique(v$species_clean)),sum(!is.finite(d$avg_lifespan_change_percent)),sum(!is.finite(d$max_lifespan_change_percent)),sum(is.na(d$pubmed_id)),sum(v$flag_nonphysical),mean(v$delta),mean(pub$avg_lifespan_change_percent),mean(pub$max_lifespan_change_percent),median(pub$delta)))
write_csv(qcd,'results/drugage_qc.csv')
missingtab<-do.call(rbind,lapply(sort(unique(d$species_clean)),function(s){z<-d[d$species_clean==s,];data.frame(species=s,rows=nrow(z),paired=sum(z$paired),fraction_paired=mean(z$paired),mean_average_paired=mean(z$avg_lifespan_change_percent[z$paired],na.rm=TRUE),mean_average_unpaired=mean(z$avg_lifespan_change_percent[!z$paired],na.rm=TRUE))}));write_csv(missingtab,'results/drugage_missingness.csv')
primary<-data.frame(paper=c('anage','drugage'),p=c(models$primary$p,dr$p));primary$p_holm<-p.adjust(primary$p,'holm');write_csv(primary,'results/primary_tests.csv')
# Publication-ready original graphics; data points, scales, units and uncertainty labelled.
colmain<-'#006C70';colaccent<-'#B15428';colmuted<-'#607080'
plotfiles<-function(stem,fun,w=8,h=5){svg(paste0('results/',stem,'.svg'),width=w,height=h,family='sans');fun();dev.off();png(paste0('results/',stem,'.png'),width=w*180,height=h*180,res=180);fun();dev.off()}
plotfiles('anage_adjusted',function(){par(mar=c(5,5,3,1));rx<-resid(lm(lmaturity~lmass+factor(order),p));ry<-resid(lm(llife~lmass+factor(order),p));plot(rx,ry,pch=16,col=adjustcolor(colmain,.32),cex=.65,xlab='Maturity residual (natural-log days)',ylab='Longevity residual (natural-log years)',main='Later maturation and mammalian longevity');abline(a=0,b=models$primary$beta,col=colaccent,lwd=2);abline(h=0,v=0,col='#CCCCCC',lty=3);mtext(sprintf('%d species; mass and order adjusted; family-clustered inference',nrow(p)),side=3,line=.25,cex=.8)})
plotfiles('anage_robustness',function(){par(mar=c(5,10,3,1));z<-ma;labels<-gsub('_',' ',z$model);yy<-rev(seq_len(nrow(z)));plot(z$beta,yy,xlim=range(c(0,z$lo,z$hi)),ylim=c(.5,nrow(z)+.5),yaxt='n',ylab='',xlab='Log–log maturity coefficient (95% family-clustered CI)',pch=16,col=colmain,main='Sensitivity to adjustment and sample restrictions');axis(2,at=yy,labels=labels,las=1,cex.axis=.8);segments(z$lo,yy,z$hi,yy,col=colmain,lwd=2);abline(v=0,lty=2,col=colmuted)})
plotfiles('drugage_paired',function(){par(mar=c(5,5,3,1));plot(pub$avg_lifespan_change_percent,pub$max_lifespan_change_percent,pch=16,col=adjustcolor(colmain,.6),xlab='Average/median lifespan change (%)',ylab='Maximum / upper-tail change (%)',main='Lifespan endpoints differ within source groups');abline(a=0,b=1,lty=2,col=colaccent,lwd=2);mtext(sprintf('Each point averages paired records from one source group (n = %d)',nrow(pub)),side=3,line=.25,cex=.8)})
plotfiles('drugage_forest',function(){par(mar=c(5,9,3,1));z<-rbind(data.frame(label='All source groups',dr[,c('estimate','boot_lo','boot_hi')]),data.frame(label=c('C. elegans','D. melanogaster','M. musculus'),sp[,c('estimate','boot_lo','boot_hi')]));yy<-4:1;plot(z$estimate,yy,xlim=range(c(0,z$boot_lo,z$boot_hi)),ylim=c(.5,4.5),yaxt='n',ylab='',xlab='Maximum minus average/median change (percentage points)',pch=16,col=colmain,main='Source-balanced endpoint differences');axis(2,at=yy,labels=z$label,las=1);segments(z$boot_lo,yy,z$boot_hi,yy,col=colmain,lwd=2);abline(v=0,lty=2,col=colaccent);mtext('95% source bootstrap intervals; curated records, not all interventions',side=3,line=.25,cex=.8)})
capture.output(sessionInfo(),file='results/sessionInfo.txt')
print(qc_a);print(ma);print(cvres);print(qcd);print(sens);print(sp);print(primary)
