## ----setup, include = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(tidy.opts = list(width.cutoff = 80), tidy = TRUE)
knitr::opts_chunk$set(comment = NA, warning = FALSE)
options(width = 2000) # width console output
library(restriktor)
library(bain)

## ----echo=FALSE, results = FALSE, message = FALSE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

n.coef <- 3 # Number of dummy variables (model without intercept)

mu <- rep(0, n.coef)
intercept <- 0

r2 <- 0.15
samplesize <- 100

b.ratios <- c(3,2,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <-uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
## Check - als steeds verschil van d in betas
## es = d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
# d = betas[1] - betas[2]
#k <- n.coef
#i = 1:k
#d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
# ES = 0.2660913
# 
# Ik doe nu niets met de correlaties van -1... Alleen de var van 1
#
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.2660913

set.seed(123) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit


## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1 vs complement
H1 <- "D1 > D2 > D3" # Denoting mu1 > mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_1c <- goric(fit, hypotheses = list(H1 = H1), comparison = "complement")
results_1c

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1 vs unconstrained (default)
H1 <- "D1 > D2 > D3" # mu1 > mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_1u <- goric(fit, hypotheses = list(H1 = H1))
results_1u

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1, H2, and unconstrained (default)
H1 <- "D1 > D2 > D3" # mu1 > mu2 > mu3
H2 <- "D1 < D2 < D3" # mu1 < mu2 < mu3

# Apply GORIC #
set.seed(123) 
results_12u <- goric(fit, hypotheses = list(H1 = H1, H2 = H2))
results_12u
round(results_12u$ratio.gw, 3)

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1, H2, and unconstrained (default)
H1 <- "D1 > D2" # # => H1: D1 > D2, D3 # mu1 > mu2, mu3
H2 <- "D1 < D2" # # => H2: D1 < D2, D3 # mu1 < mu2, mu3

# Apply GORIC #
set.seed(123) 
results_12 <- goric(fit, hypotheses = list(H1, H2), comparison = "none")
results_12
round(results_12$ratio.gw, 3)

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,2.5,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')
sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.2508602

set.seed(123) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[,1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# General hypothesis vs complement
H1g <- "D1 > D2 > D3" # mu1 >  mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_12 <- goric(fit, hypotheses = list(H1g = H1g), comparison = "complement")
results_12
round(results_12$ratio.gw, 3)

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# More specific hypothesis vs complement
H1m <- "D1 - D2 > .1, D2 > D3" # mu1 - mu2 > .1, mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_12 <- goric(fit, hypotheses = list(H1m = H1m), comparison = "complement")
results_12
round(results_12$ratio.gw, 3)

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,2,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
## Check - als steeds verschil van d in betas
## es = d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
# d = betas[1] - betas[2]
#k <- n.coef
#i = 1:k
#d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
# ES = 0.2660913

set.seed(123) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[,1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1, H2, and unconstrained (default) - subset true
H1 <- "D1 > D2 > D3"               # mu1 > mu2 > mu3
H2 <- "D1 > D2" # H2: D1 > D2, D3  # mu1 > mu2,  mu3

# Apply GORIC #
set.seed(123) 
results_12u <- goric(fit, hypotheses = list(H1, H2))
results_12u
round(results_12u$ratio.gw, 3)

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,2,3)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.0972037

set.seed(123) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[,1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1, H2, and unconstrained (default) - non-overlapping part true
H1 <- "D1 > D2 > D3"               # mu1 > mu2 > mu3
H2 <- "D1 > D2" # H2: D1 > D2, D3  # mu1 > mu2,  mu3

# Apply GORIC #
set.seed(123) 
results_12u <- goric(fit, hypotheses = list(H1, H2))
results_12u
round(results_12u$ratio.gw, 3)

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# non-overlapping part vs  its complement
Hn <- "D1 > D2 < D3" # mu1 > mu2 < mu3

# Apply GORIC #
set.seed(123) 
results_n <- goric(fit, hypotheses = list(Hn = Hn), comparison = "complement")
results_n

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,2,2)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if(sum(abs(b.ratios)) == 0){ # all zeros (thus under H0)
  betas <- b.ratios
} else{
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.1121875

set.seed(123) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1 vs complement
H1 <- "D1 > D2 > D3" # mu1 > mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_1c <- goric(fit, hypotheses = list(H1), comparison = "complement")
results_1c
#
coef(fit)

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Check on three boundary hypotheses
Hb1 <- "D1 = D2 > D3" # mu1 = mu2 > mu3
Hb2 <- "D1 > D2 = D3" # mu1 > mu2 = mu3
Hb3 <- "D1 = D2 = D3" # mu1 = mu2 = mu3
  
# Apply GORIC #
set.seed(123) 
results_b123 <- goric(fit, hypotheses = list(Hb1 = Hb1, Hb2 = Hb2, Hb3 = Hb3))
results_b123

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# H1, H2, and unconstrained (default) - non-overlapping part true
H1 <- "D1 > D2 > D3" # mu1 > mu2 > mu3
H2 <- "D1 > D2 < D3" # mu1 > mu2 < mu3

# Apply GORIC #
set.seed(123) 
results_12u <- goric(fit, hypotheses = list(H1, H2))
results_12u
round(results_12u$ratio.gw, 3)
round(results_12u$ratio.pw, 3)

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# overlap H1 and H2 versus its complement
Hb <- "D1 > D2, -.1 < D2 - D3 < .1" # mu1 > mu2, -.1 < mu2 - mu3 < .1

# Apply GORIC #
set.seed(123) 
results_b <- goric(fit, hypotheses = list(Hb = Hb), comparison = "complement")
results_b

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,3,1.5)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if(sum(abs(b.ratios)) == 0){ # all zeros (thus under H0)
  betas <- b.ratios
} else{
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.1642296

set.seed(125) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 = D2 > 0, D3 > 0" # mu1 = mu2 > 0, mu3 > 0

# Apply GORIC #
set.seed(123) 
results_0c <- goric(fit, hypotheses = list(H1), comparison = "complement")
results_0c
round(results_0c$ratio.lw, 3) # ratio of loglik weights
round(results_0c$ratio.gw, 3) # ratio of GORIC weights

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "-0.2 < D1 - D2 < .2, D1 > 0, D2 > 0, D3 > 0"
# -0.2 < mu1 - mu2 < .2, mu1 > 0, mu2 > 0, mu3 > 0

# Apply GORIC #
set.seed(123) 
results_12 <- goric(fit, hypotheses = list(H1), comparison = "complement")
results_12
round(results_12$ratio.lw, 3) # ratio of loglik weights
round(results_12$ratio.gw, 3) # ratio of GORIC weights
coef(fit)

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(1.7,1.6,1.5)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else{
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
#betas
#[1] 0.5668791 0.5335333 0.5001875
#> 0.5668791 - 0.5335333; 0.5335333 - 0.5001875
#[1] 0.0333458
#[1] 0.0333458
#d = betas[1] - betas[2]
#k <- n.coef
#i = 1:k
#d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
#[1] 0.02722676
#
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))

set.seed(124) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 > D2 > D3" # mu1 > mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_2c <- goric(fit, hypotheses = list(H1 = H1), comparison = "complement")
results_2c
round(results_2c$ratio.gw, 3)

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,2,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
## Check - als steeds verschil van d in betas
## es = d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
# d = betas[1] - betas[2]
#k <- n.coef
#i = 1:k
#d * (sqrt(sum((2*i-1-k)^2)))/(2*sqrt(k))
# ES = 0.2660913
#
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))

set.seed(125) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 - D2 > .5, D2 > D3" # mu1 - mu2 > .5, mu2 > mu3

# Apply GORIC #
set.seed(123) 
results_3c <- goric(fit, hypotheses = list(H1 = H1), comparison = "complement")
results_3c
round(results_3c$ratio.gw, 3)

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,3,1.5)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
#install.packages('fastDummies')
#library(fastDummies)
#sample <- dummy_cols(sample, select_columns = 'D')

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.1642296

set.seed(125) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(bain)

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 = D2 > 0 & D3 > 0" 

# Apply bain
set.seed(123)
bain1_1 <- bain(fit, H1)
bain1_1
set.seed(123)
bain1_2 <- bain(fit, H1, fraction = 2)     
bain1_2
set.seed(123)
bain1_3 <- bain(fit, H1, fraction = 3)     
bain1_3

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
b.ratios <- c(3,2.5,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D
sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if(sum(abs(b.ratios)) == 0){ # all zeros (thus under H0)
  betas <- b.ratios
}else{
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <-uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.2508602

set.seed(125) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 = D2 > 0 & D3 > 0" 

# Apply bain
set.seed(123)
bain1_1 <- bain(fit, H1)
bain1_1

# Apply GORIC #
set.seed(123) 
results_12 <- goric(fit, hypotheses = list(H1), comparison = "complement")
results_12

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
samplesize <- 1000

b.ratios <- c(3,2.5,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.2508602

set.seed(125) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 = D2 > 0 & D3 > 0" 

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Apply bain
set.seed(123)
bain1_1_N1000 <- bain(fit, H1)
bain1_1_N1000

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Apply GORIC #
set.seed(123) 
results_12_N1000 <- goric(fit, hypotheses = list(H1), comparison = "complement")
results_12_N1000

## ----echo = F, message = FALSE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
samplesize <- 10000

b.ratios <- c(3,2.5,1)

sample <- NULL

# Determine true / population beta coefficients in data generating process
D_ <- matrix(rep(1:n.coef), nrow = samplesize)
D <- as.factor(D_)
sample$D <- D

sample <- model.matrix(~sample$D-1)
colnames(sample) <- paste0("D_", seq_len(ncol(sample)))
sample <- as.data.frame(sample)

#sigma <- matrix(-0.99, nrow = n.coef, ncol = n.coef)
sigma <- matrix(-1, nrow = n.coef, ncol = n.coef)
diag(sigma) <- 1
        
# Define error variance
var.e <- 1 - r2 # (because all vars are standardized, var(resid)=(1-R2)*sigma_y)
        
if (sum(abs(b.ratios)) == 0) { # all zeros (thus under H0)
  betas <- b.ratios
} else {
  # Solve for x here
  fun <- function (x) {
    (t(b.ratios*x) %*% sigma %*% b.ratios*x) / (t(b.ratios*x) %*% sigma %*% b.ratios*x + var.e) - r2
  }
  
  x <- uniroot(fun, lower=0, upper=100)$root
  
  # Construct betas
  betas <- b.ratios*x
}
# k <- n.coef
# ES = (1/1) * sqrt((1/k) * sum((betas - mean(betas))^2))
# ES = 0.2508602

set.seed(125) # set seed: to obtain same results when you re-run it
epsilon <- rnorm(samplesize, sd=sqrt(var.e))
sample$y <- as.matrix(sample[, 1:n.coef]) %*% matrix(betas, nrow = n.coef) + epsilon

# Obtain fit
fit <- lm(y ~ 0 + D, data=sample)
#fit

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Equality versus its complement
H1 <- "D1 = D2 > 0 & D3 > 0" 

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Apply bain
set.seed(123)
bain1_1_N10000 <- bain(fit, H1)
bain1_1_N10000

## ----message = FALSE, warning = FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Apply GORIC #
set.seed(123) 
results_12_N10000 <- goric(fit, hypotheses = list(H1), comparison = "complement")
results_12_N10000