Nothing
inst/doc/Guidelines_GORIC_output.R
In restriktor: Restricted Statistical Estimation and Inference for Linear Models
## ----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 = 1200) # 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_", 1:ncol(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), 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))
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, 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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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
Try the restriktor package in your browser
Any scripts or data that you put into this service are public.
restriktor documentation built on Oct. 4, 2023, 9:13 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----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 = 1200) # 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_", 1:ncol(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), 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))
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, 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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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_", 1:ncol(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
Try the restriktor package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.