Skip to content

Instantly share code, notes, and snippets.

@strengejacke
Last active July 14, 2020 09:58
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save strengejacke/f4d4b95c74e343927efd01e91cec4950 to your computer and use it in GitHub Desktop.
Save strengejacke/f4d4b95c74e343927efd01e91cec4950 to your computer and use it in GitHub Desktop.
# Second Generation P-value
p_delta <- function(lb, ub, delta_lb, delta_ub) {
# special case: infinite CI and H0 bounds in the same direction
if ((delta_lb == -Inf & lb == -Inf) | (delta_ub == Inf & ub == Inf)) {
return(1)
}
# usual case: non-point CI & non-point Ho
# pdelta = |CI intersect Ho| / min{ |CI|, 2|Ho| }
if (delta_lb != delta_ub & lb != ub) {
if (lb > delta_ub | ub < delta_lb) {
return(0)
} else if(lb > delta_lb & ub < delta_ub){
return(1)
} else {
return(
(min(ub, delta_ub) - max(lb, delta_lb)) /
min(ub - lb, 2 * (delta_ub - delta_lb))
)
}
}
# special case 1: point CI, w/ or w/out a point H0
# pdelta = 0 if CI is inside the Ho
# pdelta = 1 if CI is inside the Ho
if (lb == ub) {
if (lb <= delta_ub & lb >= delta_lb){
return(1)
} else {
return(0)
}
}
# special case 2: point H0 & non-point CI
# pdelta = 1/2 if H0 is inside the CI
# pdelta = 0 if H0 is outside the CI
if (delta_lb == delta_ub & lb != ub) {
if (delta_lb <= ub & delta_lb >= lb) {
return(1/2)
} else {
return(0)
}
}
}
library(parameters)
data(qol_cancer)
model <- lm(QoL ~ time + age * education, data = qol_cancer)
# ROPE range detected as -1.99, 1.99
equivalence_test(model, rule = "classic", p_values = TRUE)
# CI for parameter age:educationmid, p-value approx. % in ROPE
p_delta(.75, 2.14, -1.99, 1.99)
# CI for parameter time, p-value approx. % in ROPE
p_delta(-.74, 2.52, -1.99, 1.99)
# CI for parameter age, p-value approx. % in ROPE
p_delta(-1.57, -.41, -1.99, 1.99)
@strengejacke
Copy link
Author

Second Generation P-value

p_delta <- function(lb, ub, delta_lb, delta_ub) {
  
  # special case: infinite CI and H0 bounds in the same direction
  if ((delta_lb == -Inf & lb == -Inf) | (delta_ub == Inf & ub == Inf)) {
    return(1)
  }
  
  # usual case: non-point CI & non-point Ho
  # pdelta = |CI intersect Ho| / min{ |CI|, 2|Ho| }
  if (delta_lb != delta_ub & lb != ub) {
    if (lb > delta_ub | ub < delta_lb) {
      return(0)
    } else if(lb > delta_lb & ub < delta_ub){
      return(1)
    } else {
      return(
        (min(ub, delta_ub) - max(lb, delta_lb)) /
          min(ub - lb, 2 * (delta_ub - delta_lb))
      )
    }
  }
  
  # special case 1: point CI, w/ or w/out a point H0
  # pdelta = 0 if CI is inside the Ho
  # pdelta = 1 if CI is inside the Ho
  if (lb == ub) {
    if (lb <= delta_ub & lb >= delta_lb){
      return(1)
    } else {
      return(0)
    }
  }
  
  # special case 2: point H0 & non-point CI
  # pdelta = 1/2 if H0 is inside the CI
  # pdelta = 0 if H0 is outside the CI
  if (delta_lb == delta_ub & lb != ub) {
    if (delta_lb <= ub & delta_lb >= lb) {
      return(1/2)
    } else {
      return(0)
    }
  }
}


library(parameters)
data(qol_cancer)
model <- lm(QoL ~ time + age * education, data = qol_cancer)

# ROPE range detected as -1.99, 1.99
equivalence_test(model, rule = "classic", p_values = TRUE)
#> # TOST-test for Practical Equivalence
#> 
#>   ROPE: [-1.99 1.99]
#> 
#>          Parameter        H0 inside ROPE        90% CI      p
#>        (Intercept)  Rejected      0.00 % [57.05 66.31] < .001
#>               time Undecided     83.81 % [-0.74  2.52]  0.658
#>                age  Accepted    100.00 % [-1.57 -0.41] > .999
#>       educationmid  Rejected      0.00 % [ 7.12 14.61] < .001
#>      educationhigh  Rejected      0.00 % [12.07 20.87] < .001
#>   age:educationmid  Rejected     88.81 % [ 0.75  2.14]  0.235
#>  age:educationhigh  Accepted    100.00 % [ 0.18  1.93] > .999

# CI for parameter age:educationmid, p-value approx. % in ROPE 
p_delta(.75, 2.14, -1.99, 1.99)
#> [1] 0.8920863

# CI for parameter time, p-value approx. % in ROPE 
p_delta(-.74, 2.52, -1.99, 1.99)
#> [1] 0.8374233

# CI for parameter age, p-value approx. % in ROPE 
p_delta(-1.57, -.41, -1.99, 1.99)
#> [1] 1

Created on 2020-07-14 by the reprex package (v0.3.0)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment