Grouping of base graphics commands

[nokome]

The following code chunk is from https://hub.stenci.la/elife/article-43154/snapshots/15

```{r}
#' @width 18
#' @height 20

require(RColorBrewer)

# The number of malaria patients
N = 10^6
P_coma = 0.34
P_anaemia = 0.24
P_G6PDdef = 0.15

PIs = seq(P_anaemia, 0.5, length.out = 100)
ORcoma = ORanaemia = array(dim=length(PIs))
TrueOR_anaemia = array(dim=length(PIs))

# odds of G6PDd in controls
O1 = P_G6PDdef/(1-P_G6PDdef)

for(i in 1:length(PIs)){
  
  # The Probability of being anaemic if you are G6PD deficient
  P_anaemia_def = PIs[i] 
  
  # We solve the equation to work out the Probability of being anaemic if you are G6PD normal
  # This is dependent on the previous probabilities (simple algebra)
  P_anaemia_norm = (P_anaemia - P_G6PDdef*P_anaemia_def)/(1-P_G6PDdef)
  
  # We assume that G6PD status and Coma status are independent
  G6PDstatus = sample(c('Normal','Def'), size = N, replace = T, 
                      prob = c(1-P_G6PDdef, P_G6PDdef))
  Comastatus = sample(c('No Coma','Coma'), size = N, replace = T, 
                      prob = c(1-P_coma, P_coma))
  # Generate anaemia status dependent on G6PD status
  Anaemiastatus = array(dim = N)
  normals = G6PDstatus=='Normal'
  defs = !normals
  Anaemiastatus[normals] = sample(x = c('No Anaemia','Anaemia'), 
                                  size = sum(normals), 
                                  replace = T, 
                                  prob = c(1-P_anaemia_norm,P_anaemia_norm))
  Anaemiastatus[defs] = sample(x = c('No Anaemia','Anaemia'), 
                               size = sum(defs),
                               replace = T, 
                               prob = c(1-P_anaemia_def,P_anaemia_def))
  
  Study_dat = data.frame(Coma = Comastatus,
                         G6PD = G6PDstatus,
                         Anaemia = Anaemiastatus)
  
  P_G6PDd_Anaemia = PIs[i]*P_G6PDdef/P_anaemia
  TrueOR_anaemia[i] =  P_G6PDd_Anaemia/(1-P_G6PDd_Anaemia)/O1
  # now subselect only those without both
  study_patients = xor(Comastatus=='Coma' , Anaemiastatus == 'Anaemia')
  
  Study_dat = Study_dat[study_patients,]
  
  # odds of G6PDd in cerebral malaria group
  O2 = (sum(Study_dat$Coma=='Coma' & Study_dat$G6PD=='Def')/
          sum(Study_dat$Coma=='Coma' & Study_dat$G6PD=='Normal'))
  
  # odds ratio for G6PDd between cases and controls
  ORcoma[i] = O2/O1
  
  
  # odds of G6PDd in SMA group
  O2 = sum(Study_dat$Anaemia=='Anaemia' & Study_dat$G6PD=='Def')/
    sum(Study_dat$Anaemia=='Anaemia' & Study_dat$G6PD=='Normal')
  # odds ratio for G6PDd between cases and controls: SMA
  ORanaemia[i] = O2/O1
  
  
}
Results = list(ORanaemia=ORanaemia,
               ORcoma=ORcoma,
               TrueOR_anaemia=TrueOR_anaemia)
               
                 
bluecols = brewer.pal(3, 'Blues')
redcols = brewer.pal(3, 'Reds')

#reportcols = brewer.pal(3, 'Accent')
par(las = 1, bty='n', mfrow=c(1,2))


plot(Results$TrueOR_anaemia, Results$ORanaemia, 
     xlab = '',
     type='l', lwd=3, col =bluecols[3], 
     ylab='Observed odds ratio for G6PDd in SMA versus controls', 
     ylim = c(0.7,1.8), xlim=c(1,2))
title('Severe malarial anaemia (SMA)')
mtext(text = 'True odds ratio for G6PDd\nin SMA versus controls',
      side = 1,line = 3.5)
axis(2, at = c(0.7,.8), labels = c(0.7,''))
polygon(c(-1,20,20,-1),c(1.22,1.22,1.8,1.8),
        col=adjustcolor(bluecols[1],alpha.f = 0.4), border = NA)
polygon(c(1.22,1.22,1.8,1.8),c(-1,20,20,-1),
        col=adjustcolor(bluecols[1],alpha.f = 0.4), border = NA)
abline(h = 1.48, col= bluecols[2], lwd=2, lty=2)
abline(v = 1.48, col= bluecols[2], lwd=2, lty=2)

lines(Results$TrueOR_anaemia, Results$ORanaemia,
      lwd=3, col =bluecols[3])
ind = which(Results$ORanaemia>1.22 & 
              Results$ORanaemia<1.8 & 
              Results$ORcoma > 0.69 & 
              Results$ORcoma <0.98)
abline(h=1)
#xs = range(Results$TrueOR_anaemia[ind])

############
plot(Results$TrueOR_anaemia, Results$ORcoma,
     lwd=3, col =redcols[3],
     xlab = '',
     type='l', xlim=c(1,2),
     ylab='Observed odds ratio for G6PDd in CM versus controls', 
     ylim = c(0.7,1.8))
title('Cerebral malaria (CM)')
mtext(text = 'True odds ratio for G6PDd\nin SMA versus controls',
      side = 1,line = 3.5)
axis(2, at = c(0.7,.8), labels = c(0.7,''))
polygon(c(-1,20,20,-1),c(0.69,0.69,0.98,0.98),
        col=adjustcolor(redcols[1],alpha.f = 0.4), border = NA)

polygon(c(1.22,1.22,1.8,1.8),c(-1,20,20,-1),
        col=adjustcolor(bluecols[1],alpha.f = 0.4), border = NA)
abline(h = c(0.82), col= redcols[2], lwd=2, lty=2)
abline(v = 1.48, col= bluecols[2], lwd=2, lty=2)

abline(h=1)
lines(Results$TrueOR_anaemia, Results$ORcoma,
      lwd=3, col =redcols[3])
```

![](article.rmd.media/fig2.jpg)

:::
{#fig2}

When this is executed in R / RStudio, the second panel has error bacnds: https://elife.stencila.io/article-43154/v15/article.xml.media/fig2.jpg. But when executed in Rasta, these bands are missing.

This is related to a previous fix 0fe994d