ACME: average causal mediation effects
ADE: average direct effects
So, proportion of mediation = ACME / (ACME + ADE)
I use G and R estimated using Flexsurv package because there independently estimated from t0.
===========
Mediation test on natural isolates. short version
h qin
June 13, 2017
rm(list=ls())
#setwd("~/github/0.network.aging.prj.bmc/0a.rls.fitting")
setwd("~/github/bmc_netwk_aging_manuscript/R1/0.nat.rls.fitting")
library('flexsurv')## Loading required package: survivalsource("../lifespan.r")Parse strains from files
files = list.files(path="../qinlab_rls/", pattern="rls.tab")
tmp1 = gsub("\\d{6}.", "", files)
redundant_strains = gsub(".rls.tab", "", tmp1)
strains = sort( unique( redundant_strains ))
strains## [1] "101S" "BY4716" "BY4741" "BY4742"
## [5] "BY4743" "JSBY4741" "M1-2" "M13"
## [9] "M14" "M2-8" "M22" "M32"
## [13] "M34" "M5" "M8" "RM112N"
## [17] "S288c" "SGU57" "sir2D.4741a" "sir2D.4742"
## [21] "sir2DSIR2.4742" "SK1" "W303" "YPS128"
## [25] "YPS163"
Take files from natural isolates
my.strains=c("101S", "M1-2","M13","M14","M2-8","M22","M32","M34","M5","M8","RM112N","S288c","SGU57", "YPS128","YPS163")
files2=c();
for( i in 1:length(my.strains)){
files2 = c( files2, files[grep(my.strains[i], files)]);
}
report = data.frame(cbind(my.strains))
report$samplesize = NA; report$R=NA; report$t0=NA; report$n=NA; report$G=NA; report$longfilename=NA;
files = files2;
strains = my.strains; Now, fit all RLS data sets by strains
for( i in 1:length(report[,1])){
#for( i in 3:4){
my.files = files[grep(strains[i], files)]
report$longfilename[i] = paste(my.files, collapse = "::");
tb = read.table( paste("../qinlab_rls/",my.files[1],sep=''), sep="\t")
if( length(my.files)> 1){
for( fi in 2:length(my.files)) {
tmp.tb = read.table( paste("../qinlab_rls/",my.files[fi],sep=''), sep="\t")
tb = rbind( tb, tmp.tb)
}
}
report$samplesize[i] = length(tb[,1])
GompFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'gompertz')
WeibFlex = flexsurvreg(formula = Surv(tb[,1]) ~ 1, dist = 'weibull')
report$avgLS[i] = mean(tb[,1])
report$stdLS[i] = sd(tb[,1])
report$CV[i] = report$stdLS[i] / report$avgLS[i]
report$GompGFlex[i] = GompFlex$res[1,1]
report$GompRFlex[i] = GompFlex$res[2,1]
report$GompLogLikFlex[i] = round(GompFlex$loglik, 1)
report$GompAICFlex[i] = round(GompFlex$AIC)
report$WeibShapeFlex[i] = WeibFlex$res[1,1]
report$WeibRateFlex[i] = WeibFlex$res[2,1]
report$WeibLogLikFlex[i] = round(WeibFlex$loglik, 1)
report$WeibAICFlex[i] = round(WeibFlex$AIC)
#set initial values
Rhat = report$GompRFlex[i]; # 'i' was missing. a bug costed HQ a whole afternoon.
Ghat = report$GompGFlex[i];
nhat = 6;
t0= (nhat-1)/Ghat;
fitBinom = optim ( c(Rhat, t0, nhat), llh.binomialMortality.single.run,
lifespan=tb[,1],
#method='SANN') #SANN needs control
method="L-BFGS-B",
lower=c(1E-10, 1, 4), upper=c(1,200,20) );
report[i, c("R", "t0", "n")] = fitBinom$par[1:3]
report$G[i] = (report$n[i] - 1)/report$t0[i]
}
report2 = report; Mediation test on Gflex <–t0 <– RFlex
Hong thinks the results indicate the t0 is the mediator from Flex to GFlex, but not sure.
library(mediation)## Loading required package: MASS## Loading required package: Matrix## Loading required package: mvtnorm## Loading required package: sandwich## mediation: Causal Mediation Analysis
## Version: 4.4.5set.seed(20170801)
report2$log10GompRFlex = log10(report2$GompRFlex)
med.fit = lm(t0 ~ log10GompRFlex, data=report2)
summary(med.fit)##
## Call:
## lm(formula = t0 ~ log10GompRFlex, data = report2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2238 -7.6956 -0.6106 1.5195 22.5871
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 99.564 20.269 4.912 0.000284 ***
## log10GompRFlex 19.967 7.429 2.688 0.018617 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.928 on 13 degrees of freedom
## Multiple R-squared: 0.3572, Adjusted R-squared: 0.3078
## F-statistic: 7.225 on 1 and 13 DF, p-value: 0.01862out.fit = lm(GompGFlex ~ t0 + log10GompRFlex, data=report2)
summary(out.fit)##
## Call:
## lm(formula = GompGFlex ~ t0 + log10GompRFlex, data = report2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.008037 -0.004385 -0.001282 0.001773 0.012763
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1818701 0.0236582 7.687 5.64e-06 ***
## t0 -0.0020169 0.0001916 -10.529 2.05e-07 ***
## log10GompRFlex -0.0095046 0.0063994 -1.485 0.163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006857 on 12 degrees of freedom
## Multiple R-squared: 0.9447, Adjusted R-squared: 0.9355
## F-statistic: 102.5 on 2 and 12 DF, p-value: 2.861e-08med.out <- mediate(med.fit, out.fit, treat = "log10GompRFlex", mediator = "t0", robustSE = TRUE, sims = 100)
summary(med.out)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.04383 -0.08096 -0.01 <2e-16 ***
## ADE -0.00718 -0.02469 0.01 0.46
## Total Effect -0.05101 -0.08276 -0.02 <2e-16 ***
## Prop. Mediated 0.86767 0.46191 1.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 15
##
##
## Simulations: 100Mediation test 2 on Rflex <–t0 <– GFlex
Hong thinks this is negative result, which means t0 works only in one direction.
med.fit = lm(t0 ~ GompGFlex, data=report2)
summary(med.fit)##
## Call:
## lm(formula = t0 ~ GompGFlex, data = report2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6401 -1.9424 -0.8670 -0.0658 8.1513
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.999 3.723 25.52 1.72e-12 ***
## GompGFlex -427.325 31.369 -13.62 4.50e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.168 on 13 degrees of freedom
## Multiple R-squared: 0.9345, Adjusted R-squared: 0.9295
## F-statistic: 185.6 on 1 and 13 DF, p-value: 4.503e-09out.fit = lm(log10GompRFlex ~ t0 + GompGFlex, data=report2)
summary(out.fit)##
## Call:
## lm(formula = log10GompRFlex ~ t0 + GompGFlex, data = report2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52342 -0.24509 0.04331 0.22083 0.34492
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002954 2.387477 -0.001 0.999
## t0 -0.017838 0.024884 -0.717 0.487
## GompGFlex -16.337520 10.999937 -1.485 0.163
##
## Residual standard error: 0.2843 on 12 degrees of freedom
## Multiple R-squared: 0.457, Adjusted R-squared: 0.3666
## F-statistic: 5.051 on 2 and 12 DF, p-value: 0.02562med.out <- mediate(med.fit, out.fit, treat = "GompGFlex", mediator = "t0", robustSE = TRUE, sims = 100)
summary(med.out)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 6.454 -14.468 31.46 0.74
## ADE -15.052 -45.810 8.83 0.18
## Total Effect -8.598 -18.553 -1.06 0.04 *
## Prop. Mediated -0.526 -9.794 3.54 0.74
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 15
##
##
## Simulations: 100
No comments:
Post a Comment