powercalculation

# Pre-operative echocardiography for predicting perioperative cardiorespiratory adverse events in children with IPAH undergoing cardiac catheterisation
# Tim Dawes
# May 2023

# Code to simulate sample sizes:

library(colorRamps)
library(Hmisc)
library(RColorBrewer)

cols<- c(brewer.pal(3,"Reds")[-1], "black", rev(brewer.pal(3, "Blues")[-1]))


      incidences<- matrix(c(3,4,1,2,4,0.001,31,55,10,42,40,12,1,0.0001,3,2,0.0001,4,35,11,56,24,26,54), ncol=4, byrow=F)
      baseline.incidence<- ((incidences[,1] + incidences[,2]) / 70)
      alpha<- 0.05
      sampleDist<- function(n, bi) {sample(x = 0:6, n, replace=T, prob = c(1-bi, rep(bi/6,6)))}


# Loops to check whether effect is detected in successive trials
      df<- list()
      
      par(mfrow=c(1,1), mar=c(5,6,2,3), mgp=c(1,1,1))
      maxN<- 300
      maxTrials<- 500
      Rs<- seq(0.8,1.2, length.out=5)
      sample.size<- rep(0, length(Rs))
      options(warn=-1)
      s<- c(20,50,100,150,200,300)
      
      for (k in 1:length(Rs))
      {
        sig<- rep(0, length(s))
        
            for (i in s)
            {
              cat(".")
              
              for (j in 1:maxTrials)
              {
                
                linpred<- x<- rep(0, i)
                
                  a<- incidences[1,1]
                  b<- incidences[1,2]
                  c<- incidences[1,3]
                  d<- incidences[1,4]
                  
                  OR = (a*d) / (b*c)
                  SE = sqrt((1/a) + (1/b) + (1/c) + (1/d))
                  
                  beta = log(OR) * Rs[k]
                  
                  x = scale(sampleDist(i, baseline.incidence[1]))
            
                  
                  linpred<- x*beta + rnorm(i,0,3)
              
                pr = exp(linpred) / (1+exp(linpred))
                y = rbinom(i,1,pr)
                model = glm(y ~ x, family = binomial(link="logit"))
                model.anova<- Anova(model, type=2)
                p.value<- model.anova$`Pr(>Chisq)`
                sm = summary(model)
                if (nrow(sm$coefficients) == 1) {power == FALSE} else {power = p.value<alpha}
                
                
              sig[match(i,s)]<- sig[match(i,s)] + power
              }
            }
        

        df[[k]]<- data.frame(X=s, Y= 100 * sig / maxTrials)
        
        if (k==1) {plot(df[[k]]$X, df[[k]]$Y, type='p', lwd=1, col=cols[k], pch=19, cex=1, xlab="", ylab="", bty='n', xaxt='n', yaxt='n', xlim=c(0,maxN), ylim=c(0,100))}
        if (k>1) {for (l in 2:k) {points(df[[l]]$X, df[[l]]$Y, type='p', lwd=1, col=cols[l], pch=19, cex=1, xlab="", ylab="", bty='n', xaxt='n', yaxt='n', xlim=c(0,100), ylim=c(0,100))}}
        
      
        axis(side=1, lwd=5, at=seq(0,maxN,length.out=9), line=0, cex.axis=2, font=2)
        mtext("Number of patients in study", side=1, line=3.5, cex=3, font=2)
        axis(side=2, lwd=5, at=seq(0,100,length.out=5), line=-1, cex.axis=2, font=2, las=2)
        mtext("Power (%)", side=2, line=3, cex=3, font=2)
        lo<- loess(Y~X, span=0.75, data=df[[k]])
        sample.size[k]<- with(df[[1]], which.min(abs(predict(lo,X)-90)))
        
        new.y<- with(df[[k]], predict(lo, X))
        new.y[new.y>100]<- 100
        new.x<- df[[k]]$X
        points(new.x, new.y, type='l', lwd=7, col=cols[k])
        
      
      }
      
      
      
      
      
# Sample size needed in each group for 90% power
      
      segments(40,90,1450,90,lwd=5,col="black", lty=2)
      segments(s[sample.size[3]],0,s[sample.size[3]],100,lwd=5,col="black", lty=2)
      
      text(50, 95, "90% Power", cex=2.2, font=2)
      
      step<- 6
      text(220, (4*step)-5, "Standard deviation", cex=2, font=2)
      
      for (k in 1:length(Rs))
          {
          segments(160,(k*step)-5,280,(k*step)-5, lwd=40, col=cols[k])
          t<- c("Standard deviation 120%", "Standard deviation 110%", "Standard deviation 100%", "Standard deviation 90%", "Standard deviation 80%")[k]
          text(220,(k*step)-5, labels=t, col="white", cex=1.8)
      }