True and false interindividual differences in the physiological response to an intervention

# https://physoc.onlinelibrary.wiley.com/doi/full/10.1113/EP085070

# https://github.com/eamonn2014/True-and-false-interindividual-differences-in-the-physiological-response-to-an-intervention

# Simulate a RCT

 library(nlme)

      n <- 5000

      noise <- 100           # add noise (within person var & meas. error) to the baseline & foll. up

      beta.treatment <- -250 # all trt'd subjects exp same trt effect, so no resp - non responders!!

      # beta.treatment <- runif(n,-20,-5)  # subjects vary in response to treatment

      clin.rel.diff <- 200

      

      pop_mu <- 0            # population mean 

      pop_sd <- 200          # between person SD

      ur.eligible <- -1000   # eligibility criteria for trial

     

      y.0true <- rnorm(n, pop_mu, pop_sd)                  # true baseline

      y.0observed <- y.0true + rnorm(n, 0, 1*noise)        # observed baseline 

      eligible <- ifelse(y.0observed > ur.eligible, 1, 0)  # Recruit all n as ur.eligible tiny

      treat <- 1*(runif(n)<.5)                             # random treatment allocation

      y.1true <- y.0true + (treat*beta.treatment)          # true follow up, treated only respond

      y.1observed <- y.1true + rnorm(n, 0, 1*noise)        # observed follow up, noise added 

      delta.observed <- y.1observed - y.0observed

     

      d <- data.frame(y.0true, y.0observed, eligible, treat , beta.treatment,

                      y.1true, y.1observed, delta.observed)

     

      # prob that a member of pop observed baseline is eligible

      # pnorm(ur.eligible, mean= pop_mu, sd=sqrt(pop_sd^2 + noise^2))

      # 1- pnorm( (pop_mu - ur.eligible) / sqrt(pop_sd^2+noise^2) )  # z score calc.

      

      trial <- d[d$eligible==1,]    # select the trial subjects


# First rows of trial data

 

# Focus on the intervention arm only - not recommended!


      trt <- trial[trial$treat==1,]

      trt$diff <- trt$y.1observed - trt$y.0observed

      foo <- sort(trt[,"diff"])

      plot(foo, main="Individual changes in response in treated arm

           Suggested individual differences due entirely to regression to the mean

           and random error (within subject and measurement error)",

      ylab= "Change in response", xlab="Individual subjects",

      col=ifelse(foo > -clin.rel.diff, 'red', 'blue'))

      abline(h=0, lty=3)

      abline(h=beta.treatment)

      abline(h=-clin.rel.diff, lty=2)

      # this many were not observed to have reduced response by more than 5

      # wrongly labelled as 'non responders'

         # red, no that fail to show clin relevant difference

     # mean(foo > -clin.rel.diff)     #proportion that fail to show clinically relevant difference

      

      

      ctrl <- trial[trial$treat==0,]

      ctrl$diff <- ctrl$y.1observed - ctrl$y.0observed

      foox <- sort(ctrl[,"diff"])

      #mean(foo < -clin.rel.diff)*length(foo)   # No that show clin relevant difference

      #mean(foo < -clin.rel.diff)     #proportion 


#  Treatment arm only

 

# remember they all responded to the drug !

# so why are we chasing the blue dots!!!

 

      with(trt, plot(diff ~  y.0observed, col=ifelse(diff < -clin.rel.diff, 'blue', 'black')

                     , xlab="observed baseline",  ylab="follow up - baseline"  ,

      main="Treatment arm: Individual changes against baseline, observed responders in blue", cex.main =1))

      with(trt, abline(lm(diff ~  y.0observed)))

      with(trt, abline(h=mean(beta.treatment), lty=2))

      with(trt, abline(h=-clin.rel.diff, lty=2))

     ctr <- trial[trial$treat==0,]

     ctr$diff <- ctr$y.1observed - ctr$y.0observed

 

     with(trt, cor.test( diff,   y.0observed, method="pearson"))

 


# Control arm only

 


 

       with(ctr, plot(diff ~  y.0observed, col=ifelse(diff < -clin.rel.diff, 'blue', 'black')

                      , xlab="observed baseline",  ylab="follow up - baseline"  ,

                      main="Control arm:  Individual changes against baseline, observed responders in blue", cex.main =1))

       with(ctr, abline(lm(diff ~  y.0observed)))

       with(ctr, abline(h=mean(beta.treatment), lty=2))

       with(trt, abline(h=-clin.rel.diff, lty=2))

       with(ctr, cor.test( diff,   y.0observed, method="pearson"))

     

# Analyse the trial correctly. Estimate the treatment effect adjusting for baseline

 

      f0 <- lm(y.1observed ~ y.0observed + treat, trial)

      summary(f0)

      confint(f0)

 

 

      alpha <- 0.05

 

      x <- trial[trial$treat %in% 0,"delta.observed"]

 

      lstar <- qchisq(alpha/2, df=  length(x)-1)

      rstar <- qchisq(1-alpha/2, df=  length(x)-1)

     

      up <- sqrt((length(x)-1)*var(x)/(lstar))

      lo <- sqrt((length(x)-1)*var(x)/(rstar))

      pe <- sqrt(var(x))

     

      # ctrl arm estimate with 95% CI

      print(c(pe, lo, up), digits=3)

      

      x1 <- trial[trial$treat %in% 1,"delta.observed"]

     

      lstar <- qchisq(alpha/2, df=  length(x1)-1)

      rstar <- qchisq(1-alpha/2, df=  length(x1)-1)

     

      up <- sqrt((length(x1)-1)*var(x1)/(lstar))

      lo <- sqrt((length(x1)-1)*var(x1)/(rstar))

      pe <- sqrt(var(x1))

     

      # trt arm estimate with 95% CI

      print(c(pe, lo, up), digits=3)

      

 

# Typical true interindividual variation in response. Adjust for the influence of biological variation and measurement error (removal of noise).


      # True individual response to the intervention

      sqrt(sd(x1)^2-sd(x)^2)  # can be -ve if more var in control group

      

      # LMM approach

      m1 <- lme(delta.observed~ treat + y.0observed,

          random=~1|treat , data=trial, method="REML",

          weights = varIdent(form = ~1 | treat))

     

      m0 <-lme(delta.observed~ treat + y.0observed,

          random=~1|treat , data=trial, method="REML")

      

      print(m1)

      anova(m1,m0) # are the trt ctr interindividual variation in response different?

     

      c.grp <- m1$sigma

      t.grp <- coef(m1$modelStruct$varStruct, uncons = FALSE)[[1]]*m1$sigma

     

      # true individual response to the intervention estimate

      sqrt(t.grp^2 - c.grp^2) 

      

      # truth

      sd(beta.treatment )


留言

張貼留言

這個網誌中的熱門文章

Bayesian reanalysis of ORBITA

import pymc3 as pm import numpy as np import arviz as az # Percutaneous coronary intervention in stable angina (ORBITA): a double-blind, randomised controlled trial 18 https://www.sciencedirect.com/science/article/abs/pii/S0140673617327149 diff = 16.6 diff_lower = -8.9 diff_upper = 42 # The scale for the Cauchy distribution cauchy_scale = 5 # Model setup with pm.Model() as model:     # Prior distribution for the difference     diff_prior = pm.Cauchy("diff", alpha=0, beta=cauchy_scale)      # Likelihood (normal approximation of the distribution of the difference)     se_diff = (diff_upper - diff_lower) / 3.92     likelihood = pm.Normal("likelihood", mu=diff_prior, sigma=se_diff, observed=np.array([diff]))     # Sample from the posterior distribution     trace = pm.sample(500000, cores=1, tune=200000, chains=4) # Calculate the posterior probability that the difference is negative post_prob_diff_negative = np.mean(t...
閱讀完整內容

Bayesian reanalysis with pymc3

 # Empagliflozin after Acute Myocardial Infarction 24 https://www.nejm.org/doi/full/10.1056/NEJMoa2314051 Do a Bayesian reanalysis of a RCT for log(hazard ratio) with 500000 samplings with a weakly informative prior, calculate the probability of hazard ratio less than 1 and a sensitivity analysis using pymc3 and only show the results in traditional Chinese instead of the code: A total of 3260 patients were assigned to receive empagliflozin and 3262 to receive placebo. During a median follow-up of 17.9 months, a first hospitalization for heart failure or death from any cause occurred in 267 patients (8.2%) in the empagliflozin group and in 298 patients (9.1%) in the placebo group, with incidence rates of 5.9 and 6.6 events, respectively, per 100 patient-years (hazard ratio, 0.90; 95% confidence interval [CI], 0.76 to 1.06; P=0.21). 使用 pymc3 對一項 RCT 進行貝氏再分析，在弱信息先驗條件下抽樣 500000 次，計算危害比小於 1 的概率，並進行敏感性分析，只用繁體中文顯示結果，而不顯示代碼:  Do a Bayesian reanalysis of a RCT with 500000 samplings wit...
閱讀完整內容

Bayesian reanalysis of Apolipoprotein A1 Infusions

Apolipoprotein A1 Infusions and Cardiovascular Outcomes after Acute Myocardial Infarction 24 https://www.nejm.org/doi/full/10.1056/NEJMoa2400969 P(θ | D) = (P(D | θ) * P(θ)) / P(D) marginal likelihood P(D) = ∫ P(D | θ) * P(θ) dθ Likelihood: P(D | θ) ~ N(μ_l, σ_l^2) Likelihood ratio: (P(D | θ1))/(P(D | θ2)*) Bayes factor: (P(D | θ1)* P(θ1))/(P(D | θ2)* P(θ2)) Prior: P(θ) ~ N(μ_p, σ_p^2) P(θ | D) ∝ P(D | θ) * P(θ) When we normalize the product of the two PDFs by dividing by the marginal likelihood p(D), the constant term cancels out, and we are left with a normal distribution with mean μ_new and variance σ_new^2. Prior mean 0, sd 0.5 log(HR) = log(0.93) = -0.0726 SE = (log(1.05) - log(0.81)) / (2 * 1.96) = 0.0644 Prior precision = 1 / (0.35^2) = 8.16 Likelihood precision = 1 / (0.0644^2) = 241.18 Posterior precision = 8.16 + 241.18 = 249.34 Posterior mean = (8.16 * 0 + 241.18 * -0.0726) / 249.34 = -0.0701 Posterior SD = sqrt(1 / 249.34) = 0.0634 Posterior 95% credible interval for the lo...
閱讀完整內容