Tuesday, June 13, 2017

mediator effect ion t0 between GFlex and RFlex

Update on July 13, 2018:
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.



===========

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: survival
source("../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.5
set.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.01862
out.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-08
med.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: 100

Mediation 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-09
out.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.02562
med.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