library(mlpwr)
# packages, die von mlpwr Entwickler:innen zusätzlich zu laden empfohlen werden
library(pwr)
library(tidyr)

# weitere Packages, die je nach Simulationscode nötig sind: 
library(lme4) 
library(lmerTest) 
library(tmvtnorm)
library(truncnorm)



#### ALM ####

simfun_alm <- function(N) { 
  
  var_a = 2 
  var_b = 2 #Varianzen der Prädiktoren
  cor_ab = 0.2 #Korrelation der Prädiktoren
  cov_ab = cor_ab/(sqrt(var_a*var_b)) # Korrelation in Cov umrechnen
  
  beta <- c(1, 0.2, 0.1) # Regressionsgewichte (intercept, prädiktor1, prädoktor2)
  n_pred <- length(beta)-1
  sigma <- matrix(c(var_a,cov_ab, #Kovarianzmatrix
                    cov_ab,var_b), nrow = 2, ncol = 2 )
  
  pred <- rtmvnorm(N, mean = rep(0, n_pred), sigma = sigma) #simulieren aus multivariater Normalverteilung
  
  e <- rnorm(N, mean = 0, sd = 1) # Fehlerterm
  
  x <- cbind(1, pred) #Präd + Intercept
  y <- x %*% beta  + e  # y generieren
  
  dat <-data.frame( y=y,x=x,e=e)
  
  results <- lm(y ~ x.2 + x.3, data =dat) # Modell schätzen
  sum <- summary(results)
  
  sum$coefficients["x.2","Pr(>|t|)"] < 0.05 # ist für interessierenden Parameter p < 0.05?
}

alm <- find.design(simfun = simfun_alm, 
            boundaries = list(N = c(10,500)),
            power = 0.90, surrogate = "logreg") # Funktion um Design zu finden
summary(alm) # time: 13 seconds

plot(alm)



#### log #### 

# Definition der Simulationsfunktion für logistische Regression
simfun_log <- function(N) {
  # Logistische Funktion
  logistic <- function(x) 1 / (1 + exp(-x))
  
  # Generiere Prädiktoren
  dat <- data.frame(pred1 = stats::rnorm(N), pred2 = stats::rnorm(N))
  
  # Setze Koeffizienten der Prädiktoren
  beta <- c(1.2, 0.8)
  
  # Berechne Wahrscheinlichkeiten
  prob <- logistic(as.matrix(dat) %*% beta)
  
  # Generiere abhängige Variable basierend auf Wahrscheinlichkeiten
  dat$criterion <- stats::runif(N) < prob
  
  # Fitte ein logistisches Regressionsmodell
  mod <- stats::glm(criterion ~ pred1 + pred2, data = dat, 
                    family = stats::binomial)
  
  # Überprüfe, ob der p-Wert für pred2 kleiner als 0.05 ist
  summary(mod)$coefficients["pred2", "Pr(>|z|)"] < 0.05
}

# Finde das optimale Design für die logistische Regression
log <- find.design(simfun = simfun_log, 
                   boundaries = list(N = c(10, 500)),
                   power = 0.95, surrogate = "logreg")

# Zusammenfassung der Ergebnisse
summary(log) # time: 12 seconds

plot(log)




#### LMM ####

#Machbarkeitsanalyse

# Definition der Simulationsfunktion für LMM
simfun_multilevel <- function(n.per.school, n.schools) {
  
  #### Schritt 1: Parameter festsetzen ####
  # Simulationsmodell soll sein:
  # ~ 1 + pred + (1 + pred | group) -> also ein Random-Intercept-Random-Slope-Modell mit 1 Level-1-Prädiktor
  
  # Setze Parameter für die Simulation, d.h. die Fixed Effects und Random Effect Varianzen
  params <- list(theta = c(0.5, 0, 0.5),  # SD Random Intercept, Kovarianz Intercept und Slope, SD des Slopes
                                          # sind an der Level 1 SD (sigma) relativiert 
                                          # also wäre die SD vom Slope sigma * 0.5
                 beta = c(0, 1), # fixed intercept + fixed slope 
                 sigma = 1.5) # Level-1.Residual-SD
  # simulate.formula aus lme4 benötigt Benennung
  names(params$theta) <- c("group.(Intercept)", "group.prednew.(Intercept)", "group.prednew") 
  names(params$beta) <- c("(Intercept)", "prednew")
  
  
  #### Schritt 2: Generiere Daten ####
  group <- rep(1:n.schools, each = n.per.school)
  pred <- factor(rep(c("old", "new"), n.per.school * n.schools), levels = c("old", "new"))  #Level 1 Prädiktor
  dat <- data.frame(group = group, pred = pred)
  # Simuliere abhängige Variable
  dat$y <- simulate.formula(~ pred + (1 + pred | group), newdata = dat, newparams = params)[[1]]
  
  
  #### Schritt 3: Modell schätzen und p-Wert extrahieren ####
  
  # Fitte LMM mit lmer aus lme4
  mod <- lmer(y ~ pred + (1 + pred | group), data = dat)
  
  # Berechne den p-Wert für den Prädiktor prednew
  pvalue <- summary(mod)[["coefficients"]][2, "Pr(>|t|)"] #indizierung vom p-Wert
  
  # Rückgabe TRUE, wenn p-Wert < 0.01
  pvalue < 0.01
}

# Kostenfunktion für LMM
costfun_multilevel <- function(n.per.school, n.schools) {
  100 * n.per.school + 200 * n.schools
}

# Finde das optimale Design für das LMM
MLM <- find.design(simfun = simfun_multilevel, costfun = costfun_multilevel, 
                   boundaries = list(n.per.school = c(5, 25), 
                                     n.schools = c(10, 30)), #boundaries müssen sinnvoll geschätzt werden,   
                   power = 0.95)                             #sodass Modell geschätzt werden kann

# Zusammenfassung der Ergebnisse
summary(MLM) # time: 6 minutes