---
title: "OMT and maximin solutions for K=3 independent means"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Our setting is:  K=3 independent normal means. 

We provide the code to find the OMT and maximin solution when the objective is $\Pi_{3}$ (the average power when all hypotheses are false) and the constraint is strong FWER control. Note that the code can be easily modified to implement other objective functions and other constraints (e.g, strong FDR control). 

This annotated example reproduces entries in the last row of Table 2 and in the first row of Table 5 in the main manuscript. 





# The R packages

```{r, eval=TRUE, message = FALSE, warning = FALSE}

library(Rcpp)
library(devtools)
library(DEoptim)
```

#  The functions

The functions we use for finding OMT and maximin solutions are detailed next. 

#### The density of the p-value  for a paramater value
The arguments are:  a vector of test statistics, $x$, and the parameter value $\theta$, denoted in the code by $cv$. 

The output is the vector of values of the density evaluated at $x$.  

```{r, eval=TRUE}
# density of p value under alternative
ff = function (x, cv=-2){
  return (dnorm(qnorm(x)-cv)/dnorm(qnorm(x)))
}

```

#### Auxiliary functions for the decision functions that maximize the Lagrangian

The arguments: (ffx,ffy,ffz) are the alternative densities for the objective; (mu1,mu2,mu3) are the Lagrange multiplies;  and (ggx, ggy, ggz) are the alternative densities for the constraints. Note that ggi=ffi for simple alternative OMT problems, but they differ when seeking the maximin OMT solution. 

The output in each of A,B, and C corresponds to $R_1^{\mu}, R_2^{\mu}, R_3^{\mu}$ in the main manuscript (in which 
$R_k^{\mu} = a_k - \sum_{L=0}^{K-1}\mu_L \times b_{L,k}$), for the specific setting with $K=3$ hypotheses, objective $\Pi_3$ (which defines $a_k$), and strong control of the FWER (which defines the  $b_{L,k}$'s).  


```{r, eval=TRUE}
A = function (ffx,ffy,ffz,mu1,mu2,mu3,ggx,ggy,ggz){
  ffx*ffy*ffz - (mu3 + mu2 *(ggx+ggy) + mu1*ggx*ggy)
}

B = function (ffx,ffy,ffz,mu1,mu2,ggx,ggy,ggz){
  ffx*ffy*ffz - (mu2 *ggz + mu1*ggx*ggz)
}

C = function (ffx,ffy,ffz,mu1,ggx,ggy,ggz){
  ffx*ffy*ffz - (mu1*ggy*ggz)
}
```

#### The  integral constraints for given Lagrange multipliers

This function is used in order to find the optimal Lagrange multipliers. Prior to calling it, the global variables f and g have to be defined: f is the non-null density of a p-value in the object, and g is the non-null density of a p-value in the constraint. 

The arguments (mu1,mu2,mu3) are the Lagrange multiplies. The three integral constraints using the  decision functions that maximize the Lagrangian with these multiplers are computed. 

The output is a vector of three scores, one for each integral constraint. Each score has value zero  if the integral  is equal to its allowed upper bound, and is otherwise positive. The score penalizes more severely (i.e., is more positive) if the integral  is larger than its allowed upper bound  than if it is smaller than its allower upper bound. The allowed upper bound is determined by the requirement that the FWER be controlled at level $\alpha$ (the feasibility constraints in the OMT problem). 


```{r, eval=TRUE}
#The actual implementation of the computation of the integral constraint is in cpp, called from the R function following it. 

cppFunction('NumericVector int123(NumericVector f, NumericVector g,  
            NumericVector len, double mu1, double mu2, double mu3) {
            int N = g.size();
            NumericVector ints(3);
            for(int iz = 0;iz<N;iz++){
            double fzi = f[iz];
            double gzi = g[iz];
            double vz = len[iz];
            for (int iy=iz;iy<N;iy++){
            double fyi = f[iy];
            double gyi = g[iy];
            double vy = vz*len[iy];
            for (int ix=iy;ix<N;ix++){
            double w=((iz==iy)? 0.5:1.)*((iy==ix)? 0.5:1);
            double v = vy*len[ix];
            double fxi=f[ix];
            double gxi=g[ix];
            double A = fxi*fyi*fzi - (mu3 + mu2 *(gxi+gyi) + mu1*gxi*gyi);
            double B = fxi*fyi*fzi- (mu2 *gzi + mu1*gxi*gzi);
            double C = fxi*fyi*fzi - (mu1*gyi*gzi);
            bool a = (A > 0) || (A+B > 0) || (A+B+C>0);
            bool b = a && ((B > 0) || (B+C>0));
            bool c = b && (C>0);
            ints[2] += a*w*v;
            ints[1] += (a*(gxi+gyi) + b*gzi)*w*v;
            ints[0] += (a*gxi*gyi+b*gxi*gzi+c*gyi*gzi)*w*v;
            }
            }
            }
            return ints;
            }')


integrate.ints.formaxmin = function(mus,alpha=0.05){
  mu1 = mus[1]
  mu2 = mus[2]
  mu3 = mus[3]

  ints = int123 (f,g,grid.size, mu1,mu2,mu3)
  int3null = ints[3]
  int2null = ints[2]
  int1null = ints[1]
  
  
  delta3 = alpha/6-int3null
  f.delta3 = ifelse(delta3<0, 1000*(exp(-delta3)-1), mu3*delta3)
  delta2 = alpha/2-int2null
  f.delta2 = ifelse(delta2<0, 1000*(exp(-delta2)-1), mu2*delta2)
  delta1 = alpha/2-int1null
  f.delta1 = ifelse(delta1<0, 1000*(exp(-delta1)-1), mu1*delta1)
  score = max(f.delta3, f.delta2,f.delta1)
  
  return ( score)
}
```
#### Diagnostic check for a policy that it achieves strong FWER control
The arguments are: policy (a,b,c) indicating whether the hypothesis corresponding to the (smallest, second,largest) p-value is rejected; the parameter value $\theta$ for the non-null density in the objective, denoted by $cv$; the vector of negative parameter values for which the integral constraints may be evaluated, denoted by $thetavec$. 

Ouptut report includes the maximum FWER when the integral constraint has two negative values; the FWER versus thetavec when the integral constraint has exactly one negative value; and the maximum FWER when the integral constraint has exactly one negative value. A policy with strong FWER control should have maximum FWER at most $\alpha$ for all configurations considered. 


```{r, eval=TRUE}
strongcontrolCheck = function(a, b, c, cv=-2, thetavec = seq(-3,0,0.1)){
  theta = thetavec
  ff.x.theta = matrix(NA,nrow = length(ffx), ncol = length(theta))
  ff.y.theta = matrix(NA,nrow = length(ffx), ncol = length(theta))
  ff.z.theta = matrix(NA,nrow = length(ffx), ncol = length(theta))
  for (i in 1:length(theta)){
    ff.x.theta[,i] = ff(xx,theta[i])
    ff.y.theta[,i] = ff(yy,theta[i])
    ff.z.theta[,i] = ff(zz,theta[i])
  }
  
  theta1 = rep (theta,length(theta))
  theta2  = rep (theta,each = length(theta))
  #print(cbind(theta1,theta2) )
  
  FWER = TRUE
  FWER.unequal = rep(0, length(theta1))
  for (i in 1:length(theta1)){
    #    print(i)  
    ffx.1 = ff.x.theta[,which(theta==theta1[i])]
    ffx.2 = ff.x.theta[,which(theta==theta2[i])]
    ffy.1 = ff.y.theta[,which(theta==theta1[i])]
    ffy.2 = ff.y.theta[,which(theta==theta2[i])]
    ffz.1 = ff.z.theta[,which(theta==theta1[i])]
    ffz.2 = ff.z.theta[,which(theta==theta2[i])]
    
    FWER.unequal[i] = sum( (a*(ffx.1*ffy.2+ffx.2*ffy.1)+b*(ffx.1*ffz.2+ffx.2*ffz.1)+c*(ffy.1*ffz.2+ffy.2*ffz.1))*wwgg)
  }
  
  #cbind(theta1,theta2,round(FWER.unequal,3))
  cat("The maximum value of FWER for exactly one null for theta combinations (",length(thetavec)^2," pairs of alternative values), \n", theta, "\n is:", max(FWER.unequal),"\n")
  
  
  
  
  #what if exactly one is nonnull? 
  FWERvec = rep(NA, length(theta))
  for (i in 1:length(theta)){
    ff.x.theta = ff(xx,theta[i])
    ff.y.theta = ff(yy,theta[i])
    ff.z.theta = ff(zz,theta[i])
    FWERvec[i] = (sum( (a*(1-b)*(ff.x.theta+ff.y.theta)+b*(ff.x.theta+ff.z.theta+ff.y.theta))*wwgg))
  }
  
  plot(theta, FWERvec*2)
  abline(0.05,0)
  
  #cbind(theta, FWERvec*2)
  cat("The maximum value of FWER for exactly one null for theta vector\n ", theta, "\n is",max(FWERvec*2))
}



```


# The grid used for numerical integration

We find the OMT policy by approximating the relevant integrals on the following grid of values. The grid is more dense near zero since the alternative density is concentrated in this region. 

```{r, eval = TRUE}

grid.flex = c(exp (log(0.0195)*101/(1:101)),((8:399)+0.5)/400)
grid.size = c(grid.flex[1],0.5*(grid.flex[-1] - grid.flex[-length(grid.flex)]))
grid.size = grid.size + c(0.5*(grid.flex[-1] - grid.flex[-length(grid.flex)]),(1-max(grid.flex)))

ng = length(grid.flex)
xx = rep (grid.flex,each = ng*ng)
yy  = rep(rep (grid.flex,each = ng),ng)
zz = rep(grid.flex,ng*ng)
gg = (grid.size %o% grid.size)%o%grid.size
nn = ng*ng*ng
jj = xx<yy | xx<zz | yy<zz
xx=xx[!jj]
yy=yy[!jj]
zz=zz[!jj]
gg=gg[!jj]
nn = length(xx)
ww = rep (1,nn)
ww[xx==yy | xx==zz | yy==zz] = 0.5
ww[xx==yy & yy==zz] = 0.25
wwgg = ww*gg

```
# The Example

We use the parameters for the last row in Table 2: 
```{r, eval = TRUE}
theta = -2 #the parameter in the last row of Table 2
alpha = 0.05

```

We seek the Lagrange multipliers that satisfy the feasibility constraints, thus solving the OMT problem. 

For average power $\Pi_3$: 

```{r, eval = TRUE}





f = ff(grid.flex, cv=theta) 
g = ff(grid.flex, cv=theta) #the density for the constraint is the same as for the objective
set.seed(1236)
#the following step takes few minus. For expository purpose, the range we seek is as follows:  mu1 between 3 and 3.5, mu2 between 0 and 0.5, and mu3 between 0 and 0.5 (we already know from a prior run with a wider range that the solution for the Lagrange multipliers should be in this range -  had this been unknown, a wider range that takes much longer to run would have been taken).  
res=DEoptim(integrate.ints.formaxmin, c(3,0,0), c(3.5,0.5,0.5), alpha=alpha, DEoptim.control(trace=FALSE)) 


mu1=res$optim$bestmem[1]
mu2=res$optim$bestmem[2]
mu3=res$optim$bestmem[3]
print(c(mu1,mu2,mu3)) #The Lagrange multipliers
```

The optimal policy is therefore given below, and  the average power $\Pi_3$ when the data is indeed generated from the objective is computed. 

```{r, eval = TRUE}
ffx = ff(xx, cv=theta)
ffy = ff(yy, cv=theta)
ffz = ff(zz, cv=theta)



A.now = A(ffx,ffy,ffz,mu1,mu2,mu3, ffx,ffy,ffz)
B.now = B(ffx,ffy,ffz,mu1,mu2,ffx,ffy,ffz)
C.now = C(ffx,ffy,ffz,mu1, ffx,ffy,ffz)                        
a = A.now>0 | (A.now+B.now)>0 | (A.now+B.now+C.now)>0
b = a & (B.now>0 | (B.now+C.now)>0)
c = b & C.now>0            
powavg= sum((a+b+c)*ffx*ffy*ffz*wwgg)*6/3
print(powavg) #the entry in the last row, third column, of Table 2. 

```

For the same policy, we  also compute the minimal power ($\Pi_{any}$):

```{r, eval = TRUE}

powmin = sum(a*ffx*ffy*ffz*wwgg)*6
print(powmin)#the entry in the last row, sixth column, of Table 2. 

```

Hommel's procedure has the following power for the same data generation

```{r, eval = FALSE}
c = xx<=alpha 
b = c | (yy<alpha/2)  
a = b |  zz<=alpha/3
powavg.hommel = sum((a+b+c)*ffx*ffy*ffz*wwgg)*6/3
powmin.hommel = sum(a*ffx*ffy*ffz*wwgg)*6
  
print(c(powavg.hommel, powmin.hommel)) #the entries in the last row, second and fifth columns, of Table 2.   

```


Next, we find the maximin OMT policy, which provides strong FWER control for the same objective above.  
For expository purpose, we already use prior runs to know that the worst case constraint (ie, the constraint that provides the policy with least value of the objective) is at the parameter value $\theta =  -2.13$. Thus, we ommit the search for the $\theta$ value in the constraint that will lead to smallest average power (since it will take too long to run). We verify that this is indeed the optimal (worst case) $\theta$  constraint, by checking that the FWER of the achieved policy is indeed at most 0.05 for any value of the parameter space, $\theta_1,\theta_2, \theta_3$, for which at least one of the $\theta_i$'s is zero.  

```{r, eval = TRUE}


#the following step takes few minus. 

theta.constraint = -2.13
g = ff(grid.flex, cv= theta.constraint) 
#seek the Lagrange multiplies when the constraint densities have expectation under the alternative of "theta.constraint"
res=DEoptim(integrate.ints.formaxmin, c(3, 0.2,0), c(3.1,0.6,0), alpha=alpha, DEoptim.control(trace=FALSE)) 
mu1=res$optim$bestmem[1]
mu2=res$optim$bestmem[2]
mu3=res$optim$bestmem[3]
print(c(mu1,mu2,mu3)) #The Lagrange multipliers


#The maximin policy is therefore given below, and  the average power, when the data is indeed generated from the objective, can be computed. 
ggx = ff(xx, cv=theta.constraint)
ggy = ff(yy, cv=theta.constraint)
ggz = ff(zz, cv=theta.constraint)



A.now = A(ffx,ffy,ffz,mu1,mu2,mu3, ggx,ggy,ggz)
B.now = B(ffx,ffy,ffz,mu1,mu2,ggx,ggy,ggz)
C.now = C(ffx,ffy,ffz,mu1, ggx,ggy,ggz)                        
a = A.now>0 | (A.now+B.now)>0 | (A.now+B.now+C.now)>0
b = a & (B.now>0 | (B.now+C.now)>0)
c = b & C.now>0            
powavg= sum((a+b+c)*ffx*ffy*ffz*wwgg)*6/3
print(3*powavg) #the entry in the first row, second column, of Table 2. 



#check that this policy provides indeed strong FWER control. This may take awhile. 
strongcontrolCheck(a,b,c, cv=theta, thetavec = seq(-3,0,0.2))

```