---
title: "PraktickÃ¡ ÄÃ¡st"
author: "Martin Farda"
date: "2023-11-23"
output: html_document
---

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

# Libraries

```{r}
library(MASS)
library(matrixcalc)
library(maxLik)
library(ggplot2)
library(grid)
library(gridExtra)
library(Cairo)
```

# GenerovÃ¡nÃ­ dat
## GenerÃ¡tor vysvÄtlujÃ­cÃ­h promÄnnÃ½ch

NÃ¡sledujÃ­cÃ­ funkce generuje dataframe vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch (mezi sebou nezÃ¡vislÃ½ch) ze zadanÃ½ch rozdÄlenÃ­

```{r}
generate_dataframe <- function(n, distributions, intercept) {
  covariates <- data.frame(matrix(nrow = n, ncol = length(distributions)))
  colnames(covariates) <- paste0("V", 1:length(distributions))
  for (i in 1:length(distributions)) {
    dist_params <- distributions[[i]]
    dist_name <- dist_params[[1]]
    dist_args <- dist_params[-1]
    dist_args <- lapply(dist_args, function(x) as.numeric(x))
    covariates[, i] <- do.call(dist_name, c(list(n = n), dist_args)) 
  }
  if (intercept == 1) {
    covariates <- cbind(Intercept = rep(1, n), covariates)
  }
  
  return(covariates)
}
```

Funkce vÃ½Å¡e generuje dataframe na zÃ¡kladÄ zadanÃ©ho poÄtu ÅÃ¡dkÅ¯ a listu zadanÃ½ch rozdÄlenÃ­ se zvolenÃ½mi parametry. KaÅ¾dÃ½ sloupec dataframeu je pak vektorem n realizacÃ­ jednoho ze zadanÃ½ch rozdÄlenÃ­. UkÃ¡zka uÅ¾itÃ­ nÃ­Å¾e

```{r}

set.seed(12345)
n <- 5000
distributions <- list(
  c("rnorm", mean = 0, sd = 1)
  #c("runif", min = -1, max = 1)  
)

df <- generate_dataframe(n, distributions,1)
summary(df)

```

## GenerÃ¡tor korelovanÃ½ch vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch (normÃ¡lnÃ­, rovnomÄrnÃ© a alternativnÃ­)

Tato funkce generuje dataframe vysvÄtlujÃ­ch promÄnnÃ½ch z normÃ¡lnÃ­ho, rovnomÄrnÃ©ho a alternativnÃ­ho rozdÄlenÃ­ se zadanÃ½mi parametry. VysvÄtlujÃ­cÃ­ promÄnnÃ© mohou bÃ½t korelovanÃ© (je zadÃ¡na korelaÄnÃ­ struktura trojrozmÄrnÃ©ho normÃ¡lnÃ­ho rozdÄlenÃ­ s nulovou stÅednÃ­ hodnotu a jednotkovÃ½mi rozptyly ze kterÃ½ch jsou transformacemi zÃ­skÃ¡ny vysvÄtlujÃ­cÃ­ promÄnnÃ© s poÅ¾adovanÃ½mi rozdÄlenÃ­mi).

```{r}
generate_correlated_dataframe <- function(n, mu, var, l, u, p, ro12, ro13, ro23, intercept) {
  require(MASS)
  cor_mat<-matrix(c(1, ro12, ro13, ro12, 1, ro23, ro13, ro23, 1), nrow = 3, ncol = 3, byrow = TRUE) ## korelaÄnÃ­ struktura
  data <- mvrnorm(n, c(0,0,0), cor_mat) ## nagenerovÃ¡nÃ­ trojrozmÄrnÃ©ho normÃ¡lnÃ­ho s pÅedepsanou korelaÄnÃ­ strukturou
  normal <- data[,1]*sqrt(var)+mu
  uniform <- l+(u-l)*pnorm(data[,2]) ## transformace normÃ¡lnÃ­ho rozdÄlenÃ­ na rovnomÄrnÃ© U(l,u)
  bernoulli <- ifelse(pnorm(data[,3])>(1-p),1,0) ## transformace normÃ¡lnÃ­ho rozdÄlenÃ­ na alternativnÃ­ Alt(p)

  
  covariates_cor <- data.frame(V1 = normal, V2 = uniform, V3 = bernoulli)
  
  if (intercept == 1) {
    covariates_cor <- cbind(Intercept = rep(1, n), covariates_cor)
  }
  
  return(covariates_cor)
}



```

OpÄt ukÃ¡zka uÅ¾itÃ­

```{r}

set.seed(12345)
mu<-1
var<-9
l<--2
u<-6
p<-0.7
ro12<-0.2
ro13<-0.2
ro23<-0.2

is.positive.definite(matrix(c(1, ro12, ro13, ro12, 1, ro23, ro13, ro23, 1), nrow = 3, ncol = 3, byrow = TRUE)) ## je zvolenÃ¡ matice pozitivnÄ definitnÃ­?

df_cor <- generate_correlated_dataframe(n, mu, var, l, u, p, ro12, ro13, ro23, 1)
summary(df_cor)

```


## GenerÃ¡tor odezvy

Odezvu budeme generovat na zÃ¡kladÄ zadanÃ©ho dataframeu vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch (kaÅ¾dÃ½ ÅÃ¡dek pÅÃ­sluÅ¡Ã­ jednomu pozorovÃ¡nÃ­) a zadanÃ©ho vektoru koeficientÅ¯ beta.

```{r}
generate_response <- function(beta, covariates) {
  p <- plogis(as.matrix(covariates) %*% beta)
  response <- rbinom(nrow(covariates),1, prob = p)
  
  return(response)
}


```

UkÃ¡zka uÅ¾itÃ­ pro zvolenou beta

```{r}
set.seed(12345)
beta <- c(0.4, 0.2)
covariates <- df ## pouÅ¾ijeme dataframe z pÅedchozÃ­ ukÃ¡zky
response <- generate_response(beta, covariates)



```

## GenerÃ¡tor chyb v odezvÄ

```{r}
generate_response_with_mistakes<-function(gamma_0,gamma_1,covariates, response) {
  s1<-plogis(as.matrix(covariates) %*% gamma_0)
  sensitivity_factor<-rbinom(nrow(covariates),1, prob = s1)
  s2<-1-plogis(as.matrix(covariates) %*% gamma_1)
  specificity_factor<-rbinom(nrow(covariates),1, prob = (1-s2))
  
  adjusted_response<-ifelse(response==1, ifelse(sensitivity_factor==1,1,0),ifelse(specificity_factor==1,1,0))
  
  return(adjusted_response)
}

```

UkÃ¡zka uÅ¾itÃ­ pro zvolen0 gamma_0, gamma_1 a v pÅedeÅ¡lÃ½ch pÅÃ­kladech generovanou odezvu a vysvÄtlujÃ­cÃ­ promÄnnÃ©

```{r}
gamma_0 <- c(1,0)
gamma_1 <- c(-1,0)
response_with_mistakes <- generate_response_with_mistakes(gamma_0,gamma_1,covariates, response)

summary(as.factor(response))
summary(as.factor(response-response_with_mistakes))

```

# EM algoritmus
## EM algoritmus pro znÃ¡mou senzitivitu a specificitu

EM algoritmus pro logistickou regresi s chybnÄ mÄÅenou odezvou a znÃ¡mou senzitivitu a specificitu. V pÅÃ­padÄ, Å¾e jde o konstantnÃ­ znÃ¡mou senzitivitu a specificitu, staÄÃ­ dosadit za sens a spec pouze konstanty. V obecnÄjÅ¡Ã­m pÅÃ­padÄ vektor funkcÄnÃ­ch hodnot senzitivity a specificity pro vÅ¡echna pozorovÃ¡nÃ­. (Tedy vstup s konstantami za sens a spec a vstup s vektory stejnÃ½ch konstant o dÃ©lce n za sens a spec generujÃ­ stejnÃ½ vÃ½stup)

```{r}

EM_with_known<-function(V, X, sens, spec, beta0 = rep(0,length(X[1,])), tol = 10^(-7)) {
  beta_k<-beta0 ## poÄÃ¡teÄnÃ­ odhad bety
  difference<-100 ## tak aby se while cyklus rozbÄhl
  
  ## funkce pro vÃ½poÄet podmÃ­nÄnÃ© stÅednÃ­ hodnoty odezvy pÅi znalosti zkaÅ¾enÃ© odezvy
  calculate_T_hat <- function(V, X, s1, s2, beta) {
  p_k <- plogis(as.matrix(X) %*% beta)
  T_hat_k <- V * (s1 * p_k / (s1 * p_k + (1 - s2) * (1 - p_k))) + (1 - V) * ((1 - s1) * p_k / ((1 - s1) * p_k + s2 * (1 - p_k)))
  return(T_hat_k)
} 
 
  T_hat_k <- calculate_T_hat(V, X, sens, spec, beta_k)
  
  
  ## Funkce pro vÃ½poÄet zÃ¡pornÃ© pozorovanÃ© log-vÄrohodnosti pro zadanÃ© beta, pouÅ¾ito pozdÄji pro zastavovacÃ­ pravidlo
  neg_log_likelihood <- function(current_beta, V, X, s1, s2) {
    neg_log_lik <- -sum(V * log(plogis(as.matrix(X) %*% current_beta) * s1 + (1 - s2) * (1 - plogis(as.matrix(X) %*% current_beta))) + (1 - V) * log(plogis(as.matrix(X) %*% current_beta) * (1 - s1) + s2 * (1 - plogis(as.matrix(X) %*% current_beta))))
    return(neg_log_lik)
  }

  while (difference>tol) {
 
  prev_neg_log_lik<-neg_log_likelihood(beta_k, V, X, sens, spec)  
    
  beta_k<-suppressWarnings(glm(T_hat_k~.-1, family = binomial(link = "logit"), data = X)$coef)
  #print(beta_k)
  
  current_neg_log_lik<-neg_log_likelihood(beta_k, V, X, sens, spec) 
  
  
  difference<-abs(prev_neg_log_lik-current_neg_log_lik)/prev_neg_log_lik
  ## difference<-max(abs(beta_k-prev_beta))/prev_beta[which.max(abs(beta_k-prev_beta))]
  
  #print(prev_neg_log_lik)
  #print(current_neg_log_lik)
  
  T_hat_k <- calculate_T_hat(V, X, sens, spec, beta_k)
}
  
  ## odhad po poslednÃ­m provedenÃ©m kroku
  beta_k<-suppressWarnings(glm(T_hat_k~.-1, family = binomial(link = "logit"), data = X)$coef)
  
  
  ## intervaly spolehlivosti
  var_matrix<- function(X, p, T_hat) {
    
  diag_matrix <- diag(as.numeric(p * (1 - p) - T_hat * (1 - T_hat)))
  X_matrix <- as.matrix(X)
  XX <- t(X_matrix) %*% diag_matrix %*% X_matrix
  
  result<-solve(XX)
  return(result)
  }
  
  p<-plogis(as.matrix(X) %*% beta_k)
  var_mat<-diag(var_matrix(X, p, T_hat_k))
  shift_for_CI<-qnorm(0.975)*sqrt(var_mat)
  
  lower_CI<-beta_k-shift_for_CI
  upper_CI<-beta_k+shift_for_CI
  
  return(c(beta_k, lower_CI, upper_CI))
}




```


UkÃ¡zka pouÅ¾itÃ­

```{r}
starttim<-Sys.time()
beta0<-c(2,2)
s1<-plogis(as.matrix(covariates) %*% gamma_0)[1,] ## skuteÄnÃ¡ konstantÃ­ senzitivita
s2<-1-plogis(as.matrix(covariates) %*% gamma_1)[1,] ## skuteÄnÃ¡ konstantÃ­ specificita
#s1<-1
#s2<-1

(params_estimate<-EM_with_known(response_with_mistakes,covariates,s1,s2,beta0))
beta_estimate<-params_estimate[1:2]

## stejnÃ½ vÃ½stup i pÅi vloÅ¾enÃ­ vektoru konstant za sens a spec jako funkÄnÃ­ch hodnot senzitivity a specificity pro kaÅ¾dÃ© pozorovÃ¡nÃ­
(params_estimate<-EM_with_known(response_with_mistakes,covariates,rep(s1,n),rep(s2,n),beta0))
beta_estimate<-params_estimate[1:2]

model<-glm(response~.-1, family = binomial(link = "logit"), data = covariates)
model ## JednoduchÃ© glm na skuteÄnou odezvu pro srovnÃ¡nÃ­
confint(model)

endtim<-Sys.time()
endtim-starttim
```
SpoÄtenÃ© intervaly spolehlivosti pokrÃ½vajÃ­ obÄ sloÅ¾ky beta, na druhou stranu vidÃ­me, Å¾e odhady jsou mÃ©nÄ pÅesnÃ© a se znaÄnÄ Å¡irÅ¡Ã­mi intervaly spolehlivosti neÅ¾ bÄÅ¾nÃ© glm a funkce confint kdybychom znali skuteÄnou odezvu. VÅ¡e tedy v souladu s oÄekÃ¡vanÃ­m.

## EM algoritmus s validaÄnÃ­ skupinou a neznÃ¡mou konstantnÃ­ sens a spec

EM algoritmus pro logistickou regresi s chybnÄ mÄÅenou odezvou s validaÄnÃ­ skupinou a neznÃ¡mou konstantnÃ­ senzitivitou a specificitou

```{r}

EM_with_unknown_constant<-function(V, Y, validation_indicator, X, beta0 = rep(0,length(X[1,])), sens = 0.8 , spec = 0.8, tol = 10^(-7)) {
  
  ## jednoduchÃ© hlÃ¡Å¡ky pro chybnÄ zadanÃ© vstupy
  if (!all(validation_indicator %in% c(0, 1))) {
        stop("Chyba: IndikÃ¡tor validaÄnÃ­ho vÃ½bÄru mÅ¯Å¾e obsahovat pouze hodnoty 0 nebo 1.")
    }
  
  if (length(V) != length(validation_indicator)) { ## jinak ÅeÄeno, misklasifikovanou odezvu musÃ­me mÃ­t k dispozici pro vÅ¡echna pozorovÃ¡nÃ­ a stejnÄ tak ten indikÃ¡tor
        stop("Chyba: DÃ©lka moÅ¾nÃ© chybnÃ© odezvy V musÃ­ odpovÃ­dat dÃ©lce vektoru indikÃ¡torÅ¯ validaÄnÃ­ skupiny.")
  }
  
  if (length(Y) != sum(validation_indicator == 1)) {
        stop("Chyba: DÃ©lka skuteÄnÃ© odezvy Y musÃ­ odpovÃ­dat poÄtu jedniÄek vektoru indikÃ¡torÅ¯ validaÄnÃ­ skupiny.")
    }
  
  ## vektory indexÅ¯ pÅÃ­sluÅ¡nÃ© nevalidaÄnÃ­mu a validaÄnÃ­mu vÃ½bÄru
  index_V <- which(validation_indicator == 0)
  index_Y <- which(validation_indicator == 1)
  
  ## poÄÃ¡teÄnÃ­ odhady
  beta_k<-beta0
  sens_k<-sens
  spec_k<-spec
  
  difference<-100 ## tak aby se nÃ¡m while cyklus rozbÄhl
  eps<-10^(-3) ## na odsazenÃ­ od 1 nebo 0 - jinak bychom se mohli snaÅ¾it dÄlit nulou
  
  ## funkce pro vÃ½poÄet podmÃ­nÄnÃ© stÅednÃ­ hodnoty skuteÄnÃ© odezvy pÅi danÃ© zkaÅ¾enÃ© odezvÄ. MÃ¡ smysl dÄlat jen pro nevalidaÄnÃ­ vÃ½bÄr
  calculate_T_hat_general <- function(V, X, index_V, beta, s1_k, s2_k) {
  p_k <- plogis(as.matrix(X[index_V,]) %*% beta)
  V_main<-V[index_V]
  
  p_k<-ifelse(p_k>1-eps, 1-eps, p_k) ## odsazenÃ­ aby nevzniklo dÄlenÃ­ nulou, z toho NA za T_hat_k
  p_k<-ifelse(p_k<eps, eps, p_k)
  s1_k<-ifelse(s1_k>1-eps, 1-eps, s1_k)
  s2_k<-ifelse(s2_k>1-eps, 1-eps, s2_k)
  
  T_hat_k <- V_main * (s1_k * p_k / (s1_k * p_k + (1 - s2_k) * (1 - p_k))) + (1 - V_main) * ((1 - s1_k) * p_k / ((1 - s1_k) * p_k + s2_k * (1 - p_k)))
  return(T_hat_k)
}
 
  T_hat_k <- calculate_T_hat_general(V, X, index_V, beta_k,  sens_k, spec_k)
  
  joint_response <- numeric(length(validation_indicator)) ## sjednocenÃ­ odezvy, rozliÅ¡enÃ­ dle validation_indicator na Y a T_hat_k
  joint_response[index_V] <- T_hat_k
  joint_response[index_Y] <- Y
  
  ## zÃ¡pornÃ¡ pozorovanÃ¡ log vÄrohodnost
  neg_log_likelihood_general <- function(current_beta, current_sens, current_spec) {
    ## zvlÃ¡Å¡Å¥ pro validaÄnÃ­ a nevalidaÄnÃ­ skupinu, pak seÄteme, Äleny co nezÃ¡visejÃ­ na ÅeÅ¡enÃ½ch parametrech vynechme
    neg_log_lik_val <- -sum(Y * log(plogis(as.matrix(X) %*% current_beta)[index_Y]) + (1-Y) * log(1-plogis(as.matrix(X) %*% current_beta)[index_Y]) + Y * V[index_Y] * log(current_sens) + Y * (1 - V[index_Y]) * log(1-current_sens) + (1-Y) * V[index_Y] * log(current_spec) + (1 - Y) * (1 - V[index_Y]) * log(1-current_spec))
    
    ## nynÃ­ pro nevalidaÄnÃ­ skupinu
    neg_log_lik_noval <- -sum((1 - V[index_V]) * log((1-current_sens) * plogis(as.matrix(X) %*% current_beta)[index_V] + (1-current_spec) * (1 - plogis(as.matrix(X) %*% current_beta)[index_V])) +  V[index_V] * log(current_sens * plogis(as.matrix(X) %*% current_beta)[index_V] + current_spec * (1 - plogis(as.matrix(X) %*% current_beta)[index_V])))
    
    neg_log_lik <- neg_log_lik_val + neg_log_lik_noval
    
    return(neg_log_lik)
  }
    
  while (difference>tol) {
    
  prev_neg_log_lik<-neg_log_likelihood_general(beta_k, sens_k, spec_k)
    
  beta_k<-suppressWarnings(glm(joint_response~.-1, family = binomial(link = "logit"), data = X)$coef)
  
  sens_k<-sum(joint_response * V)/sum(joint_response)
  spec_k<-sum((1-joint_response) * (1-V))/sum((1-joint_response))
  

  current_neg_log_lik<-neg_log_likelihood_general(beta_k, sens_k, spec_k)
  
  difference<-abs((current_neg_log_lik-prev_neg_log_lik))/(prev_neg_log_lik)
  
  T_hat_k <- calculate_T_hat_general(V, X, index_V, beta_k,  sens_k, spec_k)
  joint_response[index_V] <- T_hat_k ## aktualizovanÃ¡ podmÃ­nÄnÃ¡ stÅednÃ­ hodnota pro indexy nevalidaÄnÃ­ho vÃ½bÄru
  
  }
  
  ## odhady z poslednÃ­ho kroku
  beta_k<-suppressWarnings(glm(joint_response~.-1, family = binomial(link = "logit"), data = X)$coef)
  sens_k<-sum(joint_response * V)/sum(joint_response)
  spec_k<-sum((1-joint_response) * (1-V))/sum((1-joint_response))
 
  
  
  
  
  ## Intervaly spolehlivosti, znaÄenÃ­ ponechÃ¡no z procedury obecnÄjÅ¡Ã­ pro sens a spec zÃ¡visejÃ­cÃ­ na vysvetlujÃ­cÃ­ch promÄnnÃ½ch, pouze vzorce upraveny
  var_matrix<- function(X, p, T_hat, index_V, V, s1, s2, joint_response) {
    
  ## XX  
  diag_elements_beta<-p * (1 - p)
  diag_elements_beta[index_V]<-(p * (1 - p))[index_V] - T_hat * (1 - T_hat)
  diag_matrix_beta <- diag(as.numeric(diag_elements_beta))
  X_matrix <- as.matrix(X)
  XX <- t(X_matrix) %*% diag_matrix_beta %*% X_matrix
  
  ## Xsens
  diag_elements_Xgam0 <- -T_hat * (1 - T_hat)*(V[index_V]-s1)/(s1*(1-s1))
  diag_matrix_Xgam0 <- as.numeric(diag_elements_Xgam0)
  X_matrix_main<-as.matrix(X[index_V,])
  Xgam0 <-  t(X_matrix_main) %*% diag_matrix_Xgam0
  
  ## Xspec
  diag_elements_Xgam1<- T_hat * (1 - T_hat)*((1-V[index_V])-s2)/(s2*(1-s2))
  diag_matrix_Xgam1 <- as.numeric(diag_elements_Xgam1)
  Xgam1 <- t(X_matrix_main) %*% diag_matrix_Xgam1
  
  ## sensspec
  diag_elements_gam0gam1<- T_hat * (1 - T_hat)*(((1-V[index_V])-s2)/(s2*(1-s2))) * (V[index_V]-s1)/(s1*(1-s1))
  diag_matrix_gam0gam1 <- sum(as.numeric(diag_elements_gam0gam1))
  gam0gam1 <- diag_matrix_gam0gam1
  
  ## senssens
  diag_elements_gam0gam0<- (joint_response*(V-s1)/(s1*(1-s1)))^2
  diag_matrix_gam0gam0 <- sum(as.numeric(diag_elements_gam0gam0))
  gam0gam0 <- diag_matrix_gam0gam0
  
  ## specspec
  diag_elements_gam1gam1<- ((1-joint_response)*(-(1-V)+s2)/(s2*(1-s2)))^2
  diag_matrix_gam1gam1 <- sum(as.numeric(diag_elements_gam1gam1))
  gam1gam1 <- diag_matrix_gam1gam1
  
  col1<-rbind(XX,t(Xgam0),t(Xgam1))
  col2<-rbind(Xgam0,gam0gam0,gam0gam1)
  col3<-rbind(Xgam1,gam0gam1,gam1gam1)
  info_matrix<-cbind(col1,col2,col3)
  
  result<-solve(info_matrix)
  return(result)
  }
  
  p<-plogis(as.matrix(X) %*% beta_k)
  
  sens_for_IM<-ifelse(sens_k>1-eps, 1-eps, sens_k) ## odsazenÃ­ odhadu od 1 abychom v pozorovanÃ© informaÄnÃ­ matici nedÄlili nulou
  spec_for_IM<-ifelse(spec_k>1-eps, 1-eps, spec_k)
  
  var_mat<-diag(var_matrix(X, p, T_hat_k,index_V, V, sens_for_IM, spec_for_IM, joint_response))
  shift_for_CI<-qnorm(0.975)*sqrt(var_mat)
  
  lower_CI_beta<-beta_k-shift_for_CI[1:length(beta_k)]
  upper_CI_beta<-beta_k+shift_for_CI[1:length(beta_k)]
  lower_CI_sens<-max(sens_k-shift_for_CI[(length(beta_k)+1)],0)
  upper_CI_sens<-min(sens_k+shift_for_CI[(length(beta_k)+1)],1)
  lower_CI_spec<-max(spec_k-shift_for_CI[(length(beta_k)+2)],0)
  upper_CI_spec<-min(spec_k+shift_for_CI[(length(beta_k)+2)],1)
  
  
  
  return(rbind(beta_k, lower_CI_beta, upper_CI_beta, sens_k, lower_CI_sens, upper_CI_sens, spec_k, lower_CI_spec, upper_CI_spec))
  
  
  
  

}




```


UkÃ¡zka pouÅ¾itÃ­

```{r}
beta0<-c(2,2)

valid_prob<- 1/4 ## pst vÃ½bÄru do validaÄnÃ­ skupiny

set.seed(1234)
val_ind<-rbinom(n, 1, valid_prob) ## pÅiÅazenÃ­ hodnoty 1 prvkÅ¯m validaÄnÃ­ skupiny (prvnÃ­ch k prvkÅ¯ odpovÃ­dajÃ­cÃ­ vÃ½Å¡e zmÃ­nÄnÃ©mu podÃ­lu rozsahu) a hodnoty 0 prvkÅ¯m nevalidaÄnÃ­ skupiny

(parameters_estimates<-EM_with_unknown_constant(response_with_mistakes, response[val_ind == 1], val_ind, covariates, beta0)) ## vÄetnÄ mezÃ­ int. spol

(beta_estimate<-parameters_estimates[1,]) ## odhad parametru beta 
(sens_estimate<-parameters_estimates[4,1])
(spec_estimate<-parameters_estimates[7,1])
```


## EM algoritmus s validaÄnÃ­ skupinou

EM algoritmus pro logistickou regresi s chybnÄ mÄÅenou odezvou s validaÄnÃ­ skupinou a neznÃ¡mou senzitivitou a specificitou zÃ¡visejÃ­cÃ­ch na vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch

```{r}

EM_general<-function(V, Y, validation_indicator, X, beta0 = rep(0,length(X[1,])), gamma_0_0 = rep(0,length(X[1,])), gamma_1_0 = rep(0,length(X[1,])), tol = 10^(-7)) {
  
  ## jednoduchÃ© hlÃ¡Å¡ky pro chybnÄ zadanÃ© vstupy
  if (!all(validation_indicator %in% c(0, 1))) {
        stop("Chyba: IndikÃ¡tor validaÄnÃ­ho vÃ½bÄru mÅ¯Å¾e obsahovat pouze hodnoty 0 nebo 1.")
    }
  
  if (length(V) != length(validation_indicator)) { ## jinak ÅeÄeno, misklasifikovanou odezvu musÃ­me mÃ­t k dispozici pro vÅ¡echna pozorovÃ¡nÃ­ a stejnÄ tak ten indikÃ¡tor
        stop("Chyba: DÃ©lka moÅ¾nÃ© chybnÃ© odezvy V musÃ­ odpovÃ­dat dÃ©lce vektoru indikÃ¡torÅ¯ validaÄnÃ­ skupiny.")
  }
  
  if (length(Y) != sum(validation_indicator == 1)) {
        stop("Chyba: DÃ©lka skuteÄnÃ© odezvy Y musÃ­ odpovÃ­dat poÄtu jedniÄek vektoru indikÃ¡torÅ¯ validaÄnÃ­ skupiny.")
    }
  
  ## vektory indexÅ¯ pÅÃ­sluÅ¡nÃ© nevalidaÄnÃ­mu a validaÄnÃ­mu vÃ½bÄru
  index_V <- which(validation_indicator == 0)
  index_Y <- which(validation_indicator == 1)
  
  ## poÄÃ¡teÄnÃ­ odhady
  beta_k<-beta0
  gamma_0_k<-gamma_0_0
  gamma_1_k<-gamma_1_0
  
  difference<-100 ## tak aby se nÃ¡m while cyklus rozbÄhl
  eps<-10^(-3) ## na odsazenÃ­ od 1 nebo 0 - jinak bychom se mohli snaÅ¾it dÄlit nulou
  
  ## funkce pro vÃ½poÄet podmÃ­nÄnÃ© stÅednÃ­ hodnoty skuteÄnÃ© odezvy pÅi danÃ© zkaÅ¾enÃ© odezvÄ. MÃ¡ smysl dÄlat jen pro nevalidaÄnÃ­ vÃ½bÄr
  calculate_T_hat_general <- function(V, X, index_V, beta, gamma0, gamma1) {
  p_k <- plogis(as.matrix(X[index_V,]) %*% beta)
  s1_k<- plogis(as.matrix(X[index_V,]) %*% gamma0)
  s2_k<-1-plogis(as.matrix(X[index_V,]) %*% gamma1)
  V_main<-V[index_V]
  
  p_k<-ifelse(p_k>1-eps, 1-eps, p_k) ## odsazenÃ­ aby nevzniklo dÄlenÃ­ nulou a z toho NA za T_hat_k
  p_k<-ifelse(p_k<eps, eps, p_k)
  s1_k<-ifelse(s1_k>1-eps, 1-eps, s1_k)
  s2_k<-ifelse(s2_k>1-eps, 1-eps, s2_k)
  
  T_hat_k <- V_main * (s1_k * p_k / (s1_k * p_k + (1 - s2_k) * (1 - p_k))) + (1 - V_main) * ((1 - s1_k) * p_k / ((1 - s1_k) * p_k + s2_k * (1 - p_k)))
  return(T_hat_k)
}
 
  T_hat_k <- calculate_T_hat_general(V, X, index_V, beta_k,  gamma_0_k, gamma_1_k)
  
  joint_response <- numeric(length(validation_indicator)) ## sjednocenÃ­ odezvy, rozliÅ¡enÃ­ dle validation_indicator na Y a T_hat_k
  joint_response[index_V] <- T_hat_k
  joint_response[index_Y] <- Y
  
  ## zÃ¡pornÃ¡ pozorovanÃ¡ log vÄrohodnost
  neg_log_likelihood_general <- function(current_beta, current_gamma0, current_gamma1) {
    ## zvlÃ¡Å¡Å¥ pro validaÄnÃ­ a nevalidaÄnÃ­ skupinu, pak seÄteme, Äleny co nezÃ¡visejÃ­ na ÅeÅ¡enÃ½ch parametrech vynechme
    neg_log_lik_val <- -sum(Y * log(plogis(as.matrix(X) %*% current_beta)[index_Y]) + (1-Y) * log(1-plogis(as.matrix(X) %*% current_beta)[index_Y]) + Y * V[index_Y] * log(plogis(as.matrix(X) %*% current_gamma0)[index_Y]) + Y * (1 - V[index_Y]) * log(1-plogis(as.matrix(X) %*% current_gamma0)[index_Y]) + (1-Y) * V[index_Y] * log(plogis(as.matrix(X) %*% current_gamma1)[index_Y]) + (1 - Y) * (1 - V[index_Y]) * log(1-plogis(as.matrix(X) %*% current_gamma1)[index_Y]))
    
    ## nynÃ­ pro nevalidaÄnÃ­ skupinu
    neg_log_lik_noval <- -sum((1 - V[index_V]) * log((1-plogis(as.matrix(X) %*% current_gamma0)[index_V]) * plogis(as.matrix(X) %*% current_beta)[index_V] + (1-plogis(as.matrix(X) %*% current_gamma1)[index_V]) * (1 - plogis(as.matrix(X) %*% current_beta)[index_V])) +  V[index_V] * log(plogis(as.matrix(X) %*% current_gamma0)[index_V] * plogis(as.matrix(X) %*% current_beta)[index_V] + plogis(as.matrix(X) %*% current_gamma1)[index_V] * (1 - plogis(as.matrix(X) %*% current_beta)[index_V])))
    
    neg_log_lik <- neg_log_lik_val + neg_log_lik_noval
    
    return(neg_log_lik)
  }
  
  
  ## zÃ¡pornÃ¡ ÃºplnÃ¡ log vÄrohodnost pro danou gamma0 - nynÃ­ netÅeba
#  neg_log_likelihood_general_gamma0 <- function(current_gamma0) {
#    neg_log_lik <- -sum((joint_response) *((V) * log(plogis(as.matrix(X) %*% current_gamma0)) + (1 - V) * log(1 - plogis(as.matrix(X) %*% current_gamma0))))
#    return(neg_log_lik)
#  }

  
  ## zÃ¡pornÃ¡ ÃºplnÃ¡ log vÄrohodnost pro danou gamma1 - nynÃ­ netÅeba
#  neg_log_likelihood_general_gamma1 <- function(current_gamma1) {
#    neg_log_lik <- -sum((1-joint_response)*((1 - V) * log(1-plogis(as.matrix(X) %*% current_gamma1)) + (V) * log( plogis(as.matrix(X) %*% current_gamma1))))
#    return(neg_log_lik)
#  }

    
  while (difference>tol) {
    
  prev_neg_log_lik<-neg_log_likelihood_general(beta_k, gamma_0_k, gamma_1_k)
    
  beta_k<-suppressWarnings(glm(joint_response~.-1, family = binomial(link = "logit"), data = X)$coef)
  
  
  #beta_k<-glm(joint_response~.-1, family = binomial(link = "logit"), data = X)$coef
  #gamma_0_k<-optim(gamma_0_k,neg_log_likelihood_general_gamma0)$par
  #gamma_1_k<-optim(gamma_1_k,neg_log_likelihood_general_gamma1)$par
  
  ## pÅes funkci glm lze poÄÃ­tat i vÃ¡Å¾enÃ¡ logistickÃ¡ regrese
  gamma_0_k<-suppressWarnings(glm(V ~.-1, family = binomial(link = "logit"), data = X, weights = joint_response)$coef)
  gamma_1_k<-suppressWarnings(glm(V ~.-1, family = binomial(link = "logit"), data = X, weights = I(1-joint_response))$coef)
  

  current_neg_log_lik<-neg_log_likelihood_general(beta_k, gamma_0_k, gamma_1_k)
  
  difference<-abs((current_neg_log_lik-prev_neg_log_lik))/(prev_neg_log_lik)
  
  T_hat_k <- calculate_T_hat_general(V, X, index_V, beta_k,  gamma_0_k, gamma_1_k)
  joint_response[index_V] <- T_hat_k ## aktualizovanÃ¡ podmÃ­nÄnÃ¡ stÅednÃ­ hodnota pro indexy nevalidaÄnÃ­ho vÃ½bÄru
  
  }
  
  beta_k<-suppressWarnings(glm(joint_response~.-1, family = binomial(link = "logit"), data = X)$coef)
  gamma_0_k<-suppressWarnings(glm(V ~.-1, family = binomial(link = "logit"), data = X, weights = joint_response)$coef)
  gamma_1_k<-suppressWarnings(glm(V ~.-1, family = binomial(link = "logit"), data = X, weights = I(1-joint_response))$coef)
 
  
  
  
  
  ## Intervaly spolehlivosti
  var_matrix<- function(X, p, T_hat, index_V, V, s1, s2, joint_response) {
    
  ## XX  
  diag_elements_beta<-p * (1 - p)
  diag_elements_beta[index_V]<-(p * (1 - p))[index_V] - T_hat * (1 - T_hat)
  diag_matrix_beta <- diag(as.numeric(diag_elements_beta))
  X_matrix <- as.matrix(X)
  XX <- t(X_matrix) %*% diag_matrix_beta %*% X_matrix
  
  ## Xgam0
  diag_elements_Xgam0<- -T_hat * (1 - T_hat)*(V[index_V]-s1[index_V])
  diag_matrix_Xgam0 <- diag(as.numeric(diag_elements_Xgam0))
  X_matrix_main<-as.matrix(X[index_V,])
  Xgam0 <- t(X_matrix_main) %*% diag_matrix_Xgam0 %*% X_matrix_main
  
  ## Xgam1
  diag_elements_Xgam1<- T_hat * (1 - T_hat)*((1-V[index_V])-s2[index_V])
  diag_matrix_Xgam1 <- diag(as.numeric(diag_elements_Xgam1))
  Xgam1 <- t(X_matrix_main) %*% diag_matrix_Xgam1 %*% X_matrix_main
  
  ## gam0gam1
  diag_elements_gam0gam1<- T_hat * (1 - T_hat)*((1-V[index_V])-s2[index_V]) * (V[index_V]-s1[index_V])
  diag_matrix_gam0gam1 <- diag(as.numeric(diag_elements_gam0gam1))
  gam0gam1 <- t(X_matrix_main) %*% diag_matrix_gam0gam1 %*% X_matrix_main
  
  ## gam0gam0
  diag_elements_gam0gam0<- joint_response*s1*(1-s1)
  diag_elements_gam0gam0[index_V]<-T_hat*((s1*(1-s1))[index_V]-(1-T_hat)*((V-s1)[index_V])^2)
  diag_matrix_gam0gam0 <- diag(as.numeric(diag_elements_gam0gam0))
  gam0gam0 <- t(X_matrix) %*% diag_matrix_gam0gam0 %*% X_matrix
  
  ## gam1gam1
  diag_elements_gam1gam1<- (1-joint_response)*s2*(1-s2)
  diag_elements_gam1gam1[index_V]<-(1-T_hat)*((s2*(1-s2))[index_V]-T_hat*((1-V-s2)[index_V])^2)
  diag_matrix_gam1gam1 <- diag(as.numeric(diag_elements_gam1gam1))
  gam1gam1 <- t(X_matrix) %*% diag_matrix_gam1gam1 %*% X_matrix
  
  col1<-rbind(XX,Xgam0,Xgam1)
  col2<-rbind(Xgam0,gam0gam0,gam0gam1)
  col3<-rbind(Xgam1,gam0gam1,gam1gam1)
  info_matrix<-cbind(col1,col2,col3)
  
  result<-solve(info_matrix)
  return(result)
  }
  
  p<-plogis(as.matrix(X) %*% beta_k)
  s1<- plogis(as.matrix(X) %*% gamma_0_k)
  s2<-1-plogis(as.matrix(X) %*% gamma_1_k)
  
  s1_for_IM<-ifelse(s1>1-eps, 1-eps, s1) ## odsazenÃ­ odhadu od 1 abychom v pozorovanÃ© informaÄnÃ­ matici nedÄlili nulou
  s2_for_IM<-ifelse(s2>1-eps, 1-eps, s2)
  
  
  var_mat<-diag(var_matrix(X, p, T_hat_k,index_V, V, s1_for_IM, s2_for_IM, joint_response))
  shift_for_CI<-qnorm(0.975)*sqrt(var_mat)
  
  lower_CI_beta<-beta_k-shift_for_CI[1:length(beta_k)]
  upper_CI_beta<-beta_k+shift_for_CI[1:length(beta_k)]
  lower_CI_gam0<-gamma_0_k-shift_for_CI[(length(beta_k)+1):(length(beta_k)+length(gamma_0_k))]
  upper_CI_gam0<-gamma_0_k+shift_for_CI[(length(beta_k)+1):(length(beta_k)+length(gamma_0_k))]
  lower_CI_gam1<-gamma_1_k-shift_for_CI[(length(beta_k)+length(gamma_0_k)+1):(length(beta_k)+length(gamma_0_k)+length(gamma_1_k))]
  upper_CI_gam1<-gamma_1_k+shift_for_CI[(length(beta_k)+length(gamma_0_k)+1):(length(beta_k)+length(gamma_0_k)+length(gamma_1_k))]
  
  
  
  return(rbind(beta_k, lower_CI_beta, upper_CI_beta, gamma_0_k, lower_CI_gam0, upper_CI_gam0, gamma_1_k, lower_CI_gam1, upper_CI_gam1))
}




```


UkÃ¡zka pouÅ¾itÃ­

```{r}
starttim<-Sys.time()
(parameters_estimates<-EM_general(response_with_mistakes, response[val_ind == 1], val_ind, covariates, beta0))

beta_estimate<-parameters_estimates[1,] ## odhad parametru beta 
gamma0_estimate<-parameters_estimates[4,]
gamma1_estimate<-parameters_estimates[7,]

endtim<-Sys.time()
endtim-starttim
```


# NumerickÃ¡ maximalizace pozorovanÃ© log-vÄrohodnosti pro srovnÃ¡nÃ­
## Maximalizace log-vÄrohodnosti pÅi znÃ¡mÃ© sens, spec

PodobnÄ jako u EM algoritmu mÅ¯Å¾eme dosadit za sens a spec jen konstantu nebo vektor funkÄnÃ­ch hodnot senzitivity a specificity pro jednotlivÃ¡ pozorovÃ¡nÃ­.

```{r}
numeric_max_no_val<-function(V, X, sens, spec, beta0 = rep(0,length(X[1,])), tol = 10^(-7)){
  
  X<-as.matrix(X)
  
  neg_log_likelihood_noval <- function(bet) {
    
    neg_log_lik <-  -sum(V * (log(plogis(X %*% bet) * sens + (1 - spec) * 
                                       (1 - plogis(X %*% bet)))) + 
                            (1 - V) * (log(plogis(X %*% bet) * (1 - sens) + 
                                           spec * (1 - plogis(X %*% bet)))))
  }

  
  init_params <- beta0
  
  beta_est<-optim(init_params, neg_log_likelihood_noval, control = list(reltol = tol))$par


  
   ## intervaly spolehlivosti
  calculate_T_hat <- function(V, X, s1, s2, beta) {
  p_k <- plogis(as.matrix(X) %*% beta)
  T_hat_k <- V * (s1 * p_k / (s1 * p_k + (1 - s2) * (1 - p_k))) + (1 - V) * ((1 - s1) * p_k / ((1 - s1) * p_k + s2 * (1 - p_k)))
  return(T_hat_k)
} 
  
  T_hat <- calculate_T_hat(V, X, sens, spec, beta_est)
  
  var_matrix<- function(X, p, T_hat) {
    
  diag_matrix <- diag(as.numeric(p * (1 - p) - T_hat * (1 - T_hat)))
  X_matrix <- as.matrix(X)
  XX <- t(X_matrix) %*% diag_matrix %*% X_matrix
  
  result<-solve(XX)
  return(result)
  }
  
  p<-plogis(as.matrix(X) %*% beta_est)
  var_mat<-diag(var_matrix(X, p, T_hat))
  shift_for_CI<-qnorm(0.975)*sqrt(var_mat)
  
  lower_CI<-beta_est-shift_for_CI
  upper_CI<-beta_est+shift_for_CI
  
  
  
  
  return(c(beta_est, lower_CI, upper_CI))
}



```

```{r}
beta0<-c(2,2)

(params_estimate<-numeric_max_no_val(response_with_mistakes,covariates,s1,s2,beta0))
(beta_estimate<-params_estimate[1:2])

## stejnÃ½ vÃ½stup i pÅi vloÅ¾enÃ­ vektoru konstant za sens a spec jako funkÄnÃ­ch hodnot senzitivity a specificity pro kaÅ¾dÃ© pozorovÃ¡nÃ­
(params_estimate<-numeric_max_no_val(response_with_mistakes,covariates,rep(s1,n),rep(s2,n),beta0))
(beta_estimate<-params_estimate[1:2])

```

## Maximalizace log-vÄrohodnosti pÅi neznÃ¡mÃ© konstantnÃ­ sens, spec

```{r}
numeric_max_constant<-function(V, Y, validation_indicator, X, beta0 = rep(0,length(X[1,])), sens0 = 0.8 , spec0 = 0.8, tol = 10^(-7)){
  
  index_V <- which(validation_indicator == 0)
  index_Y <- which(validation_indicator == 1)
  
  main<-as.matrix(X[index_V,])
  valid<-as.matrix(X[index_Y,])
  
  
  
  neg_log_likelihood <- function(params) {
    bet <- params[1:length(beta0)]
    sens <- params[(length(beta0) + 1)]
    spec <- params[(length(beta0) + 2)]
    
    V_val <- V[index_Y]
    V_main <- V[index_V]


    
    
    neg_log_lik <- ifelse(length(V_main)==length(V), - 
                        sum(V_main * (log(plogis(main %*% bet) * sens + (1 - spec) * 
                                       (1 - plogis(main %*% bet)))) + 
                            (1 - V_main) * (log(plogis(main %*% bet) * (1 - sens) + 
                                           (spec) * (1 - plogis(main %*% bet))))),
                        -sum(Y * log(plogis(valid %*% bet)) + (1 - Y) * log(1 - plogis(valid %*% bet)) + 
                       Y * (V_val * log(sens) + (1 - V_val) * log(1 - sens)) + 
                             (1 - Y) * (V_val * log(1-spec) + (1 - V_val) * log(spec))) - 
                        sum(V_main * (log(plogis(main %*% bet) * sens + (1-spec) * 
                                       (1 - plogis(main %*% bet)))) + 
                            (1 - V_main) * (log(plogis(main %*% bet) * (1 - sens) + 
                                           (spec) * (1 - plogis(main %*% bet))))))
    return(neg_log_lik)
  }

  
  init_params <- c(beta0, sens0, spec0)
  
  params<-optim(init_params, neg_log_likelihood, control = list(reltol = tol))$par
  
  beta_est<-params[1:length(beta0)]
  sens_est<-params[length(beta0)+1]
  spec_est<-params[length(beta0)+2]
  
  ## Intervaly spolehlivosti, znaÄenÃ­ ponechÃ¡no z procedury obecnÄjÅ¡Ã­ pro sens a spec zÃ¡visejÃ­cÃ­ na vysvetlujÃ­cÃ­ch promÄnnÃ½ch, pouze vzorce upraveny
  eps<-10^(-3)
  
  calculate_T_hat_general <- function(V, X, index_V, beta, s1_k, s2_k) {
  p_k <- plogis(as.matrix(X[index_V,]) %*% beta)
  V_main<-V[index_V]
  
  p_k<-ifelse(p_k>1-eps, 1-eps, p_k) ## odsazenÃ­ aby nevzniklo dÄlenÃ­ nulou, z toho NA za T_hat_k
  p_k<-ifelse(p_k<eps, eps, p_k)
  s1_k<-ifelse(s1_k>1-eps, 1-eps, s1_k)
  s2_k<-ifelse(s2_k>1-eps, 1-eps, s2_k)
  
  T_hat_k <- V_main * (s1_k * p_k / (s1_k * p_k + (1 - s2_k) * (1 - p_k))) + (1 - V_main) * ((1 - s1_k) * p_k / ((1 - s1_k) * p_k + s2_k * (1 - p_k)))
  return(T_hat_k)
  }
  
  T_hat <- calculate_T_hat_general(V, X, index_V, beta_est,  sens_est, spec_est)
  
  joint_response<-V
  joint_response[index_V]<-T_hat
  joint_response[index_Y]<-Y
  
  var_matrix<- function(X, p, T_hat, index_V, V, s1, s2, joint_response) {
    
  ## XX  
  diag_elements_beta<-p * (1 - p)
  diag_elements_beta[index_V]<-(p * (1 - p))[index_V] - T_hat * (1 - T_hat)
  diag_matrix_beta <- diag(as.numeric(diag_elements_beta))
  X_matrix <- as.matrix(X)
  XX <- t(X_matrix) %*% diag_matrix_beta %*% X_matrix
  
  ## Xsens
  diag_elements_Xgam0 <- -T_hat * (1 - T_hat)*(V[index_V]-s1)/(s1*(1-s1))
  diag_matrix_Xgam0 <- as.numeric(diag_elements_Xgam0)
  X_matrix_main<-as.matrix(X[index_V,])
  Xgam0 <-  t(X_matrix_main) %*% diag_matrix_Xgam0
  
  ## Xspec
  diag_elements_Xgam1<- T_hat * (1 - T_hat)*((1-V[index_V])-s2)/(s2*(1-s2))
  diag_matrix_Xgam1 <- as.numeric(diag_elements_Xgam1)
  Xgam1 <- t(X_matrix_main) %*% diag_matrix_Xgam1
  
  ## sensspec
  diag_elements_gam0gam1<- T_hat * (1 - T_hat)*(((1-V[index_V])-s2)/(s2*(1-s2))) * (V[index_V]-s1)/(s1*(1-s1))
  diag_matrix_gam0gam1 <- sum(as.numeric(diag_elements_gam0gam1))
  gam0gam1 <- diag_matrix_gam0gam1
  
  ## senssens
  diag_elements_gam0gam0<- (joint_response^2)*((V-s1)/(s1*(1-s1)))^2
  diag_matrix_gam0gam0 <- sum(as.numeric(diag_elements_gam0gam0))
  gam0gam0 <- diag_matrix_gam0gam0
  
  ## specspec
  diag_elements_gam1gam1<- (1-joint_response)^(2)*((-(1-V)+s2)/(s2*(1-s2)))^2
  diag_matrix_gam1gam1 <- sum(as.numeric(diag_elements_gam1gam1))
  gam1gam1 <- diag_matrix_gam1gam1
  
  col1<-rbind(XX,t(Xgam0),t(Xgam1))
  col2<-rbind(Xgam0,gam0gam0,gam0gam1)
  col3<-rbind(Xgam1,gam0gam1,gam1gam1)
  info_matrix<-cbind(col1,col2,col3)
  
  result<-solve(info_matrix)
  return(result)
  }
  
  p<-plogis(as.matrix(X) %*% beta_est)
  
  sens_for_IM<-ifelse(sens_est>1-eps, 1-eps, sens_est) ## odsazenÃ­ odhadu od 1 abychom v pozorovanÃ© informaÄnÃ­ matici nedÄlili nulou
  spec_for_IM<-ifelse(spec_est>1-eps, 1-eps, spec_est)
  
  var_mat<-diag(var_matrix(X, p, T_hat,index_V, V, sens_for_IM, spec_for_IM, joint_response))
  shift_for_CI<-qnorm(0.975)*sqrt(var_mat)
  
  lower_CI_beta<-beta_est-shift_for_CI[1:length(beta_est)]
  upper_CI_beta<-beta_est+shift_for_CI[1:length(beta_est)]
  lower_CI_sens<-max(sens_est-shift_for_CI[(length(beta_est)+1)],0)
  upper_CI_sens<-min(sens_est+shift_for_CI[(length(beta_est)+1)],1)
  lower_CI_spec<-max(spec_est-shift_for_CI[(length(beta_est)+2)],0)
  upper_CI_spec<-min(spec_est+shift_for_CI[(length(beta_est)+2)],1)
  
  
  
  return(rbind(beta_est, lower_CI_beta, upper_CI_beta, sens_est, lower_CI_sens, upper_CI_sens, spec_est, lower_CI_spec, upper_CI_spec))
}



```


```{r}
(params_est<-numeric_max_constant(response_with_mistakes, response[val_ind == 1], val_ind, covariates))
(beta_est<-params_est[1,])
(sens_est<-params_est[4,1])
(spec_est<-params_est[7,1])
```

## Maximalizace log-vÄrohodnosti pÅi neznÃ¡mÃ© sens, spec zÃ¡vislÃ© na vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch

```{r}
numeric_max<-function(V, Y, validation_indicator, X, beta0 = rep(0,length(X[1,])), gamma_0_0 = rep(0,length(X[1,])), gamma_1_0 = rep(0,length(X[1,])), tol = 10^(-7)){
  
  index_V <- which(validation_indicator == 0)
  index_Y <- which(validation_indicator == 1)
  
  main<-as.matrix(X[index_V,])
  valid<-as.matrix(X[index_Y,])
  
  
  
  neg_log_likelihood <- function(params) {
    bet <- params[1:length(beta0)]
    gam0 <- params[(length(beta0) + 1):(length(beta0) + length(gamma_0_0))]
    gam1 <- params[(length(beta0) + length(gamma_0_0) + 1):length(params)]
    
    V_val <- V[index_Y]
    V_main <- V[index_V]


    neg_log_lik_val <- -sum(Y * log(plogis(valid %*% bet)) + (1 - Y) * log(1 - plogis(valid %*% bet)) + 
                       Y * (V_val * log(plogis(valid %*% gam0)) + (1 - V_val) * log(1 - plogis(valid %*% gam0))) + 
                             (1 - Y) * (V_val * log(plogis(valid %*% gam1)) + (1 - V_val) * log(1-plogis(valid %*% gam1))))
    
    neg_log_lik_main<- -sum(V_main * (log(plogis(main %*% bet) * plogis(main %*% gam0) + ( plogis(main %*% gam1)) * 
                                       (1 - plogis(main %*% bet)))) + 
                            (1 - V_main) * (log(plogis(main %*% bet) * (1 - plogis(main %*% gam0)) + 
                                           (1-plogis(main %*% gam1)) * (1 - plogis(main %*% bet)))))
    
    neg_log_lik <- neg_log_lik_val + neg_log_lik_main
    return(neg_log_lik)
  }

  
  init_params <- c(beta0, gamma_0_0, gamma_1_0)
  
  params<-optim(init_params, neg_log_likelihood, control = list(reltol = tol))$par
  
  beta_est<-params[1:length(beta0)]
  gam0_est<-params[(length(beta0)+1):(length(beta0)+length(gamma_0_0))]
  gam1_est<-params[(length(beta0) + length(gamma_0_0) + 1):length(params)]
  
  ## Intervaly spolehlivosti
  eps<-10^(-3)
  
  calculate_T_hat_general <- function(V, X, index_V, beta, gamma0, gamma1) {
  p_k <- plogis(as.matrix(X[index_V,]) %*% beta)
  s1_k<- plogis(as.matrix(X[index_V,]) %*% gamma0)
  s2_k<-1-plogis(as.matrix(X[index_V,]) %*% gamma1)
  V_main<-V[index_V]
  
  p_k<-ifelse(p_k>1-eps, 1-eps, p_k) ## odsazenÃ­ aby nevzniklo dÄlenÃ­ nulou a z toho NA za T_hat_k
  p_k<-ifelse(p_k<eps, eps, p_k)
  s1_k<-ifelse(s1_k>1-eps, 1-eps, s1_k)
  s2_k<-ifelse(s2_k>1-eps, 1-eps, s2_k)
  
  T_hat_k <- V_main * (s1_k * p_k / (s1_k * p_k + (1 - s2_k) * (1 - p_k))) + (1 - V_main) * ((1 - s1_k) * p_k / ((1 - s1_k) * p_k + s2_k * (1 - p_k)))
  return(T_hat_k)
  }
  
  T_hat <- calculate_T_hat_general(V, X, index_V, beta_est,  gam0_est, gam1_est)
  
  joint_response<-V
  joint_response[index_V]<-T_hat
  joint_response[index_Y]<-Y
  
  var_matrix<- function(X, p, T_hat, index_V, V, s1, s2, joint_response) {
    
  ## XX  
  diag_elements_beta<-p * (1 - p)
  diag_elements_beta[index_V]<-(p * (1 - p))[index_V] - T_hat * (1 - T_hat)
  diag_matrix_beta <- diag(as.numeric(diag_elements_beta))
  X_matrix <- as.matrix(X)
  XX <- t(X_matrix) %*% diag_matrix_beta %*% X_matrix
  
  ## Xgam0
  diag_elements_Xgam0<- -T_hat * (1 - T_hat)*(V[index_V]-s1[index_V])
  diag_matrix_Xgam0 <- diag(as.numeric(diag_elements_Xgam0))
  X_matrix_main<-as.matrix(X[index_V,])
  Xgam0 <- t(X_matrix_main) %*% diag_matrix_Xgam0 %*% X_matrix_main
  
  ## Xgam1
  diag_elements_Xgam1<- T_hat * (1 - T_hat)*((1-V[index_V])-s2[index_V])
  diag_matrix_Xgam1 <- diag(as.numeric(diag_elements_Xgam1))
  Xgam1 <- t(X_matrix_main) %*% diag_matrix_Xgam1 %*% X_matrix_main
  
  ## gam0gam1
  diag_elements_gam0gam1<- T_hat * (1 - T_hat)*((1-V[index_V])-s2[index_V]) * (V[index_V]-s1[index_V])
  diag_matrix_gam0gam1 <- diag(as.numeric(diag_elements_gam0gam1))
  gam0gam1 <- t(X_matrix_main) %*% diag_matrix_gam0gam1 %*% X_matrix_main
  
  ## gam0gam0
  diag_elements_gam0gam0<- joint_response*s1*(1-s1)
  diag_elements_gam0gam0[index_V]<-T_hat*((s1*(1-s1))[index_V]-(1-T_hat)*((V-s1)[index_V])^2)
  diag_matrix_gam0gam0 <- diag(as.numeric(diag_elements_gam0gam0))
  gam0gam0 <- t(X_matrix) %*% diag_matrix_gam0gam0 %*% X_matrix
  
  ## gam1gam1
  diag_elements_gam1gam1<- (1-joint_response)*s2*(1-s2)
  diag_elements_gam1gam1[index_V]<-(1-T_hat)*((s2*(1-s2))[index_V]-T_hat*((1-V-s2)[index_V])^2)
  diag_matrix_gam1gam1 <- diag(as.numeric(diag_elements_gam1gam1))
  gam1gam1 <- t(X_matrix) %*% diag_matrix_gam1gam1 %*% X_matrix
  
  col1<-rbind(XX,Xgam0,Xgam1)
  col2<-rbind(Xgam0,gam0gam0,gam0gam1)
  col3<-rbind(Xgam1,gam0gam1,gam1gam1)
  info_matrix<-cbind(col1,col2,col3)
  
  result<-solve(info_matrix)
  return(result)
  }
  
  p<-plogis(as.matrix(X) %*% beta_est)
  s1<- plogis(as.matrix(X) %*% gam0_est)
  s2<-1-plogis(as.matrix(X) %*% gam1_est)
  
  s1_for_IM<-ifelse(s1>1-eps, 1-eps, s1) ## odsazenÃ­ odhadu od 1 abychom v pozorovanÃ© informaÄnÃ­ matici nedÄlili nulou
  s2_for_IM<-ifelse(s2>1-eps, 1-eps, s2)
  
  var_mat<-diag(var_matrix(X, p, T_hat,index_V, V, s1_for_IM, s2_for_IM, joint_response))
  shift_for_CI<-qnorm(0.975)*sqrt(var_mat)
  
  lower_CI_beta<-beta_est-shift_for_CI[1:length(beta_est)]
  upper_CI_beta<-beta_est+shift_for_CI[1:length(beta_est)]
  lower_CI_gam0<-gam0_est-shift_for_CI[(length(beta_est)+1):(length(beta_est)+length(gam0_est))]
  upper_CI_gam0<-gam0_est+shift_for_CI[(length(beta_est)+1):(length(beta_est)+length(gam0_est))]
  lower_CI_gam1<-gam1_est-shift_for_CI[(length(beta_est)+length(gam0_est)+1):(length(beta_est)+length(gam0_est)+length(gam1_est))]
  upper_CI_gam1<-gam1_est+shift_for_CI[(length(beta_est)+length(gam0_est)+1):(length(beta_est)+length(gam0_est)+length(gam1_est))]
  
  
  
  return(rbind(beta_est, lower_CI_beta, upper_CI_beta, gam0_est, lower_CI_gam0, upper_CI_gam0, gam1_est, lower_CI_gam1, upper_CI_gam1))
}



```


```{r}
starttim<-Sys.time()
(params_est<-numeric_max(response_with_mistakes, response[val_ind == 1], val_ind, covariates))
(beta_est<-params_est[1,])
(sens_est<-params_est[4,])
(spec_est<-params_est[7,])

endtim<-Sys.time()
endtim-starttim
```







# PÅÃ­klady v prÃ¡ci
## PÅÃ­klad 1

Nejprve jednoduchÃ½ pÅÃ­klad pouze s jednou realizacÃ­ situace pro mapovÃ¡nÃ­ vÄrohodnosti pÅi rÅ¯znÄ velkÃ½ch rozsazÃ­ch vÃ½bÄru a pravdÄpodobnostech vÃ½bÄrÅ¯ do validaÄnÃ­ho vÃ½bÄru.

```{r}

##GenerovÃ¡nÃ­ dat
set.seed(1234)
distributions_ex1 <- list(
  c("rnorm", mean = 0, sd = 1)
) ## vezmeme jen intercept a jednu vysvÄtlujÃ­cÃ­ promÄnnou generovanou z N(0,1)

## poÄet pozorovÃ¡nÃ­ 100, 250, 500, 1000, 5000, 20000
n_ex1 <- 20000

covariates_ex1 <- generate_dataframe(n_ex1, distributions_ex1,1)

beta_ex1 <- c(0.4, 0.8)
response_ex1 <- generate_response(beta_ex1, covariates_ex1)


gamma_0_ex1 <- c(1.5,0)
gamma_1_ex1 <- c(-1.5,0)
s1_ex1<-plogis(as.matrix(covariates_ex1) %*% gamma_0_ex1)[1,] ## senzitivita a specificita
s2_ex1<-1-plogis(as.matrix(covariates_ex1) %*% gamma_1_ex1)[1,]

response_with_mistakes_ex1 <- generate_response_with_mistakes(gamma_0_ex1,gamma_1_ex1,covariates_ex1, response_ex1)


## Aplikace metod pro situaci se znÃ¡mou konstantnÃ­ senzitivitou a specificitou

## EM algoritmus

(params_EM_ex1<-EM_with_known(response_with_mistakes_ex1,covariates_ex1,s1_ex1,s2_ex1))
(beta_est_EM_ex1<-params_EM_ex1[1:2])
round(params_EM_ex1,3)


## numerickÃ¡ maximalizace pozorovanÃ© vÄrohodnosti

(params_num_ex1<-numeric_max_no_val(response_with_mistakes_ex1, covariates_ex1, s1_ex1, s2_ex1)) ## vÄetnÄ mezÃ­ pro int. spol
(beta_est_num_ex1<-params_num_ex1[1:2])
round(params_num_ex1,3)


## vÃ½poÄet zÃ¡pornÃ© log-vÄrohodnosti na mÅÃ­Å¾ce pro rÅ¯znÃ© bety pro vykreslenÃ­


neg_log_likelihood_known_cons <- function(current_beta, V, X, s1, s2) {
    neg_log_lik <- -sum(V * log(plogis(as.matrix(X) %*% current_beta) * s1 + (1 - s2) * (1 - plogis(as.matrix(X) %*% current_beta))) + (1 - V) * log(plogis(as.matrix(X) %*% current_beta) * (1 - s1) + s2 * (1 - plogis(as.matrix(X) %*% current_beta))))
    return(neg_log_lik)
  }

(true_beta_lik_ex1<-neg_log_likelihood_known_cons(beta_ex1, response_with_mistakes_ex1, covariates_ex1, s1_ex1, s2_ex1)) ## log-vÄrohodnost ve skuteÄnÃ© betÄ
(EM_beta_lik_ex1<-neg_log_likelihood_known_cons(beta_est_EM_ex1, response_with_mistakes_ex1, covariates_ex1, s1_ex1, s2_ex1)) ## log-vÄrohodnost pro odhad beta EM algoritmem
(num_beta_lik_ex1<-neg_log_likelihood_known_cons(beta_est_num_ex1, response_with_mistakes_ex1, covariates_ex1, s1_ex1, s2_ex1)) ## log-vÄrohodnost pro odhad beta pÅÃ­mo z log-vÄrohodnosti

beta1_values <- seq(-3, 3, length.out = 100) ## tvorba mÅÃ­Å¾ky
beta2_values <- seq(-3, 3, length.out = 100)
grid <- expand.grid(beta1 = beta1_values, beta2 = beta2_values)

z <- matrix(NA, nrow = length(beta1_values), ncol = length(beta2_values))
for (i in 1:length(beta1_values)) {
  for (j in 1:length(beta2_values)) {
    z[i, j] <- neg_log_likelihood_known_cons(c(beta1_values[i], beta2_values[j]), response_with_mistakes_ex1, covariates_ex1, s1_ex1, s2_ex1)
  }
} ## vÃ½poÄet hodnot pro beta na mÅÃ­Å¾ce

filtered_data <- subset(data.frame(grid, z = c(z)), is.finite(z))
```

```{r}
model<-summary(glm(response_ex1~.-1, family = binomial(link = "logit"), data = covariates_ex1))$coefficients
model ## JednoduchÃ© glm na skuteÄnou odezvu pro srovnÃ¡nÃ­
model[,"Estimate"]-model[,"Std. Error"]*qnorm(0.975)
model[,"Estimate"]+model[,"Std. Error"]*qnorm(0.975)

```

```{r}

## vykreslenÃ­

the_plot<-ggplot(filtered_data, aes(x = beta1, y = beta2, z = z)) +
  geom_contour(aes(z = z, color = after_stat(level))) +
  scale_color_gradient(limits = range(filtered_data$z), low = "green", high = "brown") + 
  #geom_text(data = subset(filtered_data, z == min(z)), aes(label = "Max vÄrohodnost na mÅÃ­Å¾ce"), color = "black", hjust = 1.3, vjust = 0.2) +
  #geom_point(data = subset(filtered_data, z == min(z)), aes(x = beta1, y = beta2), color = "black", size = 2) +
  geom_text(data = data.frame(beta1 = beta_est_EM_ex1[1], beta2 = beta_est_EM_ex1[2], z = EM_beta_lik_ex1), aes(label = "EM algoritmus"), color = "blue", hjust = -0.1, vjust = -0.7, show.legend = TRUE) +
  geom_point(data = data.frame(beta1 = beta_est_EM_ex1[1], beta2 = beta_est_EM_ex1[2], z = EM_beta_lik_ex1), aes(x = beta1, y = beta2), color = "blue", size = 2, show.legend = TRUE) +
  geom_text(data = data.frame(beta1 = beta_est_num_ex1[1], beta2 = beta_est_num_ex1[2], z = num_beta_lik_ex1), aes(label = "NumerickÃ¡ max"), color = "red", hjust = -0.1, vjust = 0.5) +
  geom_point(data = data.frame(beta1 = beta_est_num_ex1[1], beta2 = beta_est_num_ex1[2], z = num_beta_lik_ex1), aes(x = beta1, y = beta2), color = "red", size = 2) +
 geom_text(data = data.frame(beta1 = beta_ex1[1], beta2 = beta_ex1[2], z = true_beta_lik_ex1), aes(label = "SkuteÄnÃ¡ beta"), color = "orange", hjust = 1, vjust = 1.4) +
  geom_point(data = data.frame(beta1 = beta_ex1[1], beta2 = beta_ex1[2], z = true_beta_lik_ex1), aes(x = beta1, y = beta2), color = "orange", size = 2) + xlab("Beta 1") +
  ylab("Beta 2")

the_plot


```

```{r}
## uklÃ¡dÃ¡nÃ­ grafÅ¯ pro rÅ¯znÃ¡ n pÅed slouÄenÃ­m do jednoho layoutu
plot20000<-the_plot

```

```{r}
## slouÄenÃ­ a uloÅ¾enÃ­ grafu
join_plot<-grid.arrange(plot100, plot250, plot500, plot1000, plot5000, plot20000, ncol = 2, nrow = 3)
ggsave("ex1.pdf", join_plot, width = 8, height = 11)
```

## PÅÃ­klad 2


```{r}
start_time <- Sys.time()

set.seed(1234)
distributions_sim <- list(
  c("rnorm", mean = 0, sd = 1)
) ## N(0,1) rozdÄlenÃ­ pro vysvÄtlujÃ­cÃ­ promÄnnou
n_sim<-1000 ## rozsah vÃ½bÄru n
num_samples<-1000 ## generovÃ¡nÃ­ odezvy a chyby v odezvÄ pro jednu realizaci vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch. Z toho odhady skuteÄnÃ©ho parametru, 1000 realizacÃ­ provedeme


beta_sim <- beta_ex1 ## skuteÄnÃ¡ beta
gamma_0_sim <- gamma_0_ex1 ## skuteÄnÃ© gamma
gamma_1_sim <- gamma_1_ex1
s1_sim<-s1_ex1 ## skuteÄnÃ¡ sens
s2_sim<-s2_ex1 ## skuteÄnÃ¡ spec

valid_prob_sim<- 1/4 ## pst vÃ½bÄru do validaÄnÃ­ skupiny w

## pÅedpÅipravenÃ­ dataframes pro vyplÅovÃ¡nÃ­ bÄhem simulace

## EM algoritmy

beta_est_sim_EM_noval <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_EM_noval <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_EM_noval <- data.frame(beta1 = numeric(), beta2 = numeric())

beta_est_sim_EM_con <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_EM_con <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_EM_con <- data.frame(beta1 = numeric(), beta2 = numeric())
sens_est_sim_EM_con <- data.frame(sens = numeric())
sens_lower_est_sim_EM_con <- data.frame(sens = numeric())
sens_upper_est_sim_EM_con <- data.frame(sens = numeric())
spec_est_sim_EM_con <- data.frame(spec = numeric())
spec_lower_est_sim_EM_con <- data.frame(spec = numeric())
spec_upper_est_sim_EM_con <- data.frame(spec = numeric())

beta_est_sim_EM_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_EM_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_EM_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
gamma0_est_sim_EM_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
gamma0_lower_est_sim_EM_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
gamma0_upper_est_sim_EM_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
gamma1_est_sim_EM_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())
gamma1_lower_est_sim_EM_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())
gamma1_upper_est_sim_EM_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())

## numerickÃ¡ metoda

beta_est_sim_num_noval <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_num_noval <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_num_noval <- data.frame(beta1 = numeric(), beta2 = numeric())

beta_est_sim_num_con <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_num_con <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_num_con <- data.frame(beta1 = numeric(), beta2 = numeric())
sens_est_sim_num_con <- data.frame(sens = numeric())
sens_lower_est_sim_num_con <- data.frame(sens = numeric())
sens_upper_est_sim_num_con <- data.frame(sens = numeric())
spec_est_sim_num_con <- data.frame(spec = numeric())
spec_lower_est_sim_num_con <- data.frame(spec = numeric())
spec_upper_est_sim_num_con <- data.frame(spec = numeric())

beta_est_sim_num_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_num_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_num_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
gamma0_est_sim_num_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
gamma0_lower_est_sim_num_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
gamma0_upper_est_sim_num_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
gamma1_est_sim_num_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())
gamma1_lower_est_sim_num_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())
gamma1_upper_est_sim_num_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())


## pouze validaÄnÃ­ skupina, tedy bÄÅ¾nÃ¡ logistickÃ¡ regrese ÄistÄ na validaÄnÃ­ skupinu, zbytek zahodit
#beta_est_sim_val_con <- data.frame(beta1 = numeric(), beta2 = numeric())
#beta_lower_est_sim_val_con <- data.frame(beta1 = numeric(), beta2 = numeric())
#beta_upper_est_sim_val_con <- data.frame(beta1 = numeric(), beta2 = numeric())
#sens_est_sim_val_con <- data.frame(sens = numeric())
#sens_lower_est_sim_val_con <- data.frame(sens = numeric())
#sens_upper_est_sim_val_con <- data.frame(sens = numeric())
#spec_est_sim_val_con <- data.frame(spec = numeric())
#spec_lower_est_sim_val_con <- data.frame(spec = numeric())
#spec_upper_est_sim_val_con <- data.frame(spec = numeric())

beta_est_sim_val_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_lower_est_sim_val_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
beta_upper_est_sim_val_gen <- data.frame(beta1 = numeric(), beta2 = numeric())
#gamma0_est_sim_val_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
#gamma0_lower_est_sim_val_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
#gamma0_upper_est_sim_val_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric())
#gamma1_est_sim_val_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())
#gamma1_lower_est_sim_val_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())
#gamma1_upper_est_sim_val_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric())

for (i in 1:num_samples) { ## kvÅ¯li rychlosti kÃ³d spouÅ¡tÄnÃ½ se zakomentovanÃ½mi zrovna nepouÅ¾itÃ½mi ÄÃ¡stmi
  covariates_sim <- generate_dataframe(n_sim, distributions_sim,1)
  
  response_sim <- generate_response(beta_sim, covariates_sim)
  
  response_with_mistakes_sim <- generate_response_with_mistakes(gamma_0_sim,gamma_1_sim,covariates_sim, response_sim)
  
  
  beta0_sim<-runif(2, min = -2, max = 2) ## pro rÅ¯znÃ© realizace mÅ¯Å¾eme zaÄÃ­nat z rÅ¯znÃ½ch nÃ¡hodnÃ½ch poÄÃ¡teÄnÃ­ch odhadÅ¯
  gamma_0_0_sim<-runif(2, min = -2, max = 2)
  gamma_1_0_sim<-runif(2, min = -2, max = 2)
  sens_0_sim<-runif(1, min = 0.5, max = 1)
  spec_0_sim<-runif(1, min = 0.5, max = 1)
  
  val_ind_sim<-rbinom(n_sim, 1, valid_prob_sim) ## pÅiÅazenÃ­ hodnoty 1 prvkÅ¯m validaÄnÃ­ skupiny (prvnÃ­ch k prvkÅ¯ odpovÃ­dajÃ­cÃ­ vÃ½Å¡e zmÃ­nÄnÃ©mu podÃ­lu rozsahu) a hodnoty 0 prvkÅ¯m nevalidaÄnÃ­ skupiny

  ## EM algoritmus ÄÃ¡st
  
  ## EM algoritmus s konstantnÃ­mi znÃ¡mÃ½mi sens a spec
  
  #par_est_sim_EM_noval<-EM_with_known(response_with_mistakes_sim, covariates_sim, s1_sim, s2_sim, beta0_sim)
  
  #beta_est_sim_EM_noval<-rbind(beta_est_sim_EM_noval, data.frame(beta1 = par_est_sim_EM_noval[1], beta2 = par_est_sim_EM_noval[2]))
  #beta_lower_est_sim_EM_noval<-rbind(beta_lower_est_sim_EM_noval, data.frame(beta1 = par_est_sim_EM_noval[3], beta2 = par_est_sim_EM_noval[4]))
  #beta_upper_est_sim_EM_noval<-rbind(beta_upper_est_sim_EM_noval, data.frame(beta1 = par_est_sim_EM_noval[5], beta2 = par_est_sim_EM_noval[6]))
  
   ## EM algoritmus s konstantnÃ­mi neznÃ¡mÃ½mi sens a spec
  par_est_sim_EM_con<-EM_with_unknown_constant(response_with_mistakes_sim, response_sim[val_ind_sim==1], val_ind_sim, covariates_sim, beta0_sim, sens_0_sim, spec_0_sim)
  
  
  beta_est_sim_EM_con <- rbind(beta_est_sim_EM_con, data.frame(beta1 = par_est_sim_EM_con[1, 1], beta2 = par_est_sim_EM_con[1, 2]))
  beta_lower_est_sim_EM_con <- rbind(beta_lower_est_sim_EM_con, data.frame(beta1 = par_est_sim_EM_con[2, 1], beta2 = par_est_sim_EM_con[2, 2]))
  beta_upper_est_sim_EM_con <- rbind(beta_upper_est_sim_EM_con, data.frame(beta1 = par_est_sim_EM_con[3, 1], beta2 = par_est_sim_EM_con[3, 2]))
  
  #sens_est_sim_EM_con <- rbind(sens_est_sim_EM_con, data.frame(sens = par_est_sim_EM_con[4, 1]))
  #sens_lower_est_sim_EM_con <- rbind(sens_lower_est_sim_EM_con, data.frame(sens = par_est_sim_EM_con[5, 1]))
  #sens_upper_est_sim_EM_con <- rbind(sens_upper_est_sim_EM_con, data.frame(sens = par_est_sim_EM_con[6, 1]))
  
  
  #spec_est_sim_EM_con <- rbind(spec_est_sim_EM_con, data.frame(spec = par_est_sim_EM_con[7, 1]))
  #spec_lower_est_sim_EM_con <- rbind(spec_lower_est_sim_EM_con, data.frame(spec = par_est_sim_EM_con[8, 1]))
  #spec_upper_est_sim_EM_con <- rbind(spec_upper_est_sim_EM_con, data.frame(spec = par_est_sim_EM_con[9, 1]))
  
  ## EM algoritmus s sens a spec zÃ¡visejÃ­cÃ­ na vysv. promÄnnÃ½ch
  #par_est_sim_EM_gen<-EM_general(response_with_mistakes_sim, response_sim[val_ind_sim==1], val_ind_sim, covariates_sim, beta0_sim, gamma_0_0_sim, gamma_1_0_sim)
  
  
  #beta_est_sim_EM_gen <- rbind(beta_est_sim_EM_gen, data.frame(beta1 = par_est_sim_EM_gen[1, 1], beta2 = par_est_sim_EM_gen[1, 2]))
  #beta_lower_est_sim_EM_gen <- rbind(beta_lower_est_sim_EM_gen, data.frame(beta1 = par_est_sim_EM_gen[2, 1], beta2 = par_est_sim_EM_gen[2, 2]))
  #beta_upper_est_sim_EM_gen <- rbind(beta_upper_est_sim_EM_gen, data.frame(beta1 = par_est_sim_EM_gen[3, 1], beta2 = par_est_sim_EM_gen[3, 2]))
  
  #gamma0_est_sim_EM_gen <- rbind(gamma0_est_sim_EM_gen, data.frame(gamma0_1 = par_est_sim_EM_gen[4, 1], gamma0_2 = par_est_sim_EM_gen[4, 2]))
  #gamma0_lower_est_sim_EM_gen <- rbind(gamma0_lower_est_sim_EM_gen, data.frame(gamma0_1 = par_est_sim_EM_gen[5, 1], gamma0_2 = par_est_sim_EM_gen[5, 2]))
  #gamma0_upper_est_sim_EM_gen <- rbind(gamma0_upper_est_sim_EM_gen, data.frame(gamma0_1 = par_est_sim_EM_gen[6, 1], gamma0_2 = par_est_sim_EM_gen[6, 2]))
  
  
  #gamma1_est_sim_EM_gen <- rbind(gamma1_est_sim_EM_gen, data.frame(gamma1_1 = par_est_sim_EM_gen[7, 1], gamma1_2 = par_est_sim_EM_gen[7, 2]))
  #gamma1_lower_est_sim_EM_gen <- rbind(gamma1_lower_est_sim_EM_gen, data.frame(gamma1_1 = par_est_sim_EM_gen[8, 1], gamma1_2 = par_est_sim_EM_gen[8, 2]))
  #gamma1_upper_est_sim_EM_gen <- rbind(gamma1_upper_est_sim_EM_gen, data.frame(gamma1_1 = par_est_sim_EM_gen[9, 1], gamma1_2 = par_est_sim_EM_gen[9, 2]))



## numerickÃ¡ maximalizace ÄÃ¡st

 ## numerickÃ¡ max s konstantnÃ­mi znÃ¡mÃ½mi sens a spec
  
  #par_est_sim_num_noval<-numeric_max_no_val(response_with_mistakes_sim, covariates_sim, s1_sim, s2_sim, beta0_sim)
  
  #beta_est_sim_num_noval<-rbind(beta_est_sim_num_noval, data.frame(beta1 = par_est_sim_num_noval[1], beta2 = par_est_sim_num_noval[2]))
  #beta_lower_est_sim_num_noval<-rbind(beta_lower_est_sim_num_noval, data.frame(beta1 = par_est_sim_num_noval[3], beta2 = par_est_sim_num_noval[4]))
  #beta_upper_est_sim_num_noval<-rbind(beta_upper_est_sim_num_noval, data.frame(beta1 = par_est_sim_num_noval[5], beta2 = par_est_sim_num_noval[6]))
  
   ## numerickÃ¡ max s konstantnÃ­mi neznÃ¡mÃ½mi sens a spec
  par_est_sim_num_con<-suppressWarnings(numeric_max_constant(response_with_mistakes_sim, response_sim[val_ind_sim==1], val_ind_sim, covariates_sim, beta0_sim, sens_0_sim, spec_0_sim))
  
  
  beta_est_sim_num_con <- rbind(beta_est_sim_num_con, data.frame(beta1 = par_est_sim_num_con[1, 1], beta2 = par_est_sim_num_con[1, 2]))
  beta_lower_est_sim_num_con <- rbind(beta_lower_est_sim_num_con, data.frame(beta1 = par_est_sim_num_con[2, 1], beta2 = par_est_sim_num_con[2, 2]))
  beta_upper_est_sim_num_con <- rbind(beta_upper_est_sim_num_con, data.frame(beta1 = par_est_sim_num_con[3, 1], beta2 = par_est_sim_num_con[3, 2]))
  
  #sens_est_sim_num_con <- rbind(sens_est_sim_num_con, data.frame(sens = par_est_sim_num_con[4, 1]))
  #sens_lower_est_sim_num_con <- rbind(sens_lower_est_sim_num_con, data.frame(sens = par_est_sim_num_con[5, 1]))
  #sens_upper_est_sim_num_con <- rbind(sens_upper_est_sim_num_con, data.frame(sens = par_est_sim_num_con[6, 1]))
  
  
  #spec_est_sim_num_con <- rbind(spec_est_sim_num_con, data.frame(spec = par_est_sim_num_con[7, 1]))
  #spec_lower_est_sim_num_con <- rbind(spec_lower_est_sim_num_con, data.frame(spec = par_est_sim_num_con[8, 1]))
  #spec_upper_est_sim_num_con <- rbind(spec_upper_est_sim_num_con, data.frame(spec = par_est_sim_num_con[9, 1]))
  
  ## num algoritmus s sens a spec zÃ¡visejÃ­cÃ­ na vysv. promÄnnÃ½ch
  #par_est_sim_num_gen<-suppressWarnings(numeric_max(response_with_mistakes_sim, response_sim[val_ind_sim==1], val_ind_sim, covariates_sim, beta0_sim, gamma_0_0_sim, gamma_1_0_sim))
  
  
  #beta_est_sim_num_gen <- rbind(beta_est_sim_num_gen, data.frame(beta1 = par_est_sim_num_gen[1, 1], beta2 = par_est_sim_num_gen[1, 2]))
  #beta_lower_est_sim_num_gen <- rbind(beta_lower_est_sim_num_gen, data.frame(beta1 = par_est_sim_num_gen[2, 1], beta2 = par_est_sim_num_gen[2, 2]))
  #beta_upper_est_sim_num_gen <- rbind(beta_upper_est_sim_num_gen, data.frame(beta1 = par_est_sim_num_gen[3, 1], beta2 = par_est_sim_num_gen[3, 2]))
  
  #gamma0_est_sim_num_gen <- rbind(gamma0_est_sim_num_gen, data.frame(gamma0_1 = par_est_sim_num_gen[4, 1], gamma0_2 = par_est_sim_num_gen[4, 2]))
  #gamma0_lower_est_sim_num_gen <- rbind(gamma0_lower_est_sim_num_gen, data.frame(gamma0_1 = par_est_sim_num_gen[5, 1], gamma0_2 = par_est_sim_num_gen[5, 2]))
  #gamma0_upper_est_sim_num_gen <- rbind(gamma0_upper_est_sim_num_gen, data.frame(gamma0_1 = par_est_sim_num_gen[6, 1], gamma0_2 = par_est_sim_num_gen[6, 2]))
  
  
  #gamma1_est_sim_num_gen <- rbind(gamma1_est_sim_num_gen, data.frame(gamma1_1 = par_est_sim_num_gen[7, 1], gamma1_2 = par_est_sim_num_gen[7, 2]))
  #gamma1_lower_est_sim_num_gen <- rbind(gamma1_lower_est_sim_num_gen, data.frame(gamma1_1 = par_est_sim_num_gen[8, 1], gamma1_2 = par_est_sim_num_gen[8, 2]))
  #gamma1_upper_est_sim_num_gen <- rbind(gamma1_upper_est_sim_num_gen, data.frame(gamma1_1 = par_est_sim_num_gen[9, 1], gamma1_2 = par_est_sim_num_gen[9, 2]))

## bÄÅ¾nÃ¡ logistickÃ¡ regrese na validaÄnÃ­ skupinu, princip je udÄlat si z toho jen logistickou regresei na validaÄnÃ­ skupinÄ, proto je vlastnÄ zbyteÄnÄ navÃ­c poÄÃ­tat sens a spec, proto zakomentovanÃ©
  
  indexes1<-which(response_sim==1 & val_ind_sim==1)
  indexes0<-which(response_sim==0 & val_ind_sim==1)
  
  ## numerickÃ¡ max s konstantnÃ­mi neznÃ¡mÃ½mi sens a spec
  glm_beta<-summary(glm(response_sim[val_ind_sim==1]~.-1, family = binomial(link = "logit"), data = covariates_sim[val_ind_sim==1,]))$coefficients
  

  #glm_sens_con<-summary(glm(response_with_mistakes_sim[indexes1]~1, family = binomial(link = "logit"), data = covariates_sim[indexes1,]))$coefficients
  #glm_spec_con<-summary(glm(response_with_mistakes_sim[indexes0]~1, family = binomial(link = "logit"), data = covariates_sim[indexes0,]))$coefficients
  
  
  ## pouze validaÄnÃ­ vÃ½bÄr pÅi sens a spec zÃ¡visejÃ­cÃ­ na vysv. promÄnnÃ½ch

  #glm_gamma0<-summary(glm(response_with_mistakes_sim[indexes1]~.-1, family = binomial(link = "logit"), data = covariates_sim[indexes1,]))$coefficients
  #glm_gamma1<-summary(glm(response_with_mistakes_sim[indexes0]~.-1, family = binomial(link = "logit"), data = covariates_sim[indexes0,]))$coefficients
  
  
  
  par_est_sim_beta<-glm_beta[,"Estimate"]
  par_est_sim_beta_lower<-par_est_sim_beta-glm_beta[,"Std. Error"]*qnorm(0.975)
  par_est_sim_beta_upper<-par_est_sim_beta+glm_beta[,"Std. Error"]*qnorm(0.975)
  #par_est_sim_gamma0<-glm_gamma0[,"Estimate"]
  #par_est_sim_gamma0_lower<-par_est_sim_gamma0-glm_gamma0[,"Std. Error"]*qnorm(0.975)
  #par_est_sim_gamma0_upper<-par_est_sim_gamma0+glm_gamma0[,"Std. Error"]*qnorm(0.975)
  #par_est_sim_gamma1<-glm_gamma1[,"Estimate"]
  #par_est_sim_gamma1_lower<-par_est_sim_gamma1-glm_gamma1[,"Std. Error"]*qnorm(0.975)
  #par_est_sim_gamma1_upper<-par_est_sim_gamma1+glm_gamma1[,"Std. Error"]*qnorm(0.975)
  
  beta_est_sim_val_gen<-rbind(beta_est_sim_val_gen, data.frame(beta1= par_est_sim_beta[1], beta2 = par_est_sim_beta[2]))
  beta_lower_est_sim_val_gen<-rbind(beta_lower_est_sim_val_gen, data.frame(beta1= par_est_sim_beta_lower[1], beta2 = par_est_sim_beta_lower[2]))
  beta_upper_est_sim_val_gen<-rbind(beta_upper_est_sim_val_gen, data.frame(beta1= par_est_sim_beta_upper[1], beta2 = par_est_sim_beta_upper[2]))
  #gamma0_est_sim_val_gen<-rbind(gamma0_est_sim_val_gen, data.frame(gamma0_1= par_est_sim_gamma0[1], gamma0_2 = par_est_sim_gamma0[2]))
  #gamma0_lower_est_sim_val_gen<-rbind(gamma0_lower_est_sim_val_gen, data.frame(gamma0_1= par_est_sim_gamma0_lower[1], gamma0_2 = par_est_sim_gamma0_lower[2]))
  #gamma0_upper_est_sim_val_gen<-rbind(gamma0_upper_est_sim_val_gen, data.frame(gamma0_1= par_est_sim_gamma0_upper[1], gamma0_2 = par_est_sim_gamma0_upper[2]))
  #gamma1_est_sim_val_gen<-rbind(gamma1_est_sim_val_gen, data.frame(gamma1_1= par_est_sim_gamma1[1], gamma1_2 = par_est_sim_gamma1[2]))
  #gamma1_lower_est_sim_val_gen<-rbind(gamma1_lower_est_sim_val_gen, data.frame(gamma1_1= par_est_sim_gamma1_lower[1], gamma1_2 = par_est_sim_gamma1_lower[2]))
  #gamma1_upper_est_sim_val_gen<-rbind(gamma1_upper_est_sim_val_gen, data.frame(gamma1_1= par_est_sim_gamma1_upper[1], gamma1_2 = par_est_sim_gamma1_upper[2]))
}





end_time <- Sys.time()

as.numeric(end_time - start_time, units = "secs")


```

```{r}
beta_sim
#summary(beta_est_sim_EM_noval)
#round(apply(beta_est_sim_EM_noval, 2, sd),3)
summary(beta_est_sim_EM_con)
round(apply(beta_est_sim_EM_con, 2, sd),3)
#summary(beta_est_sim_EM_gen)
#round(apply(beta_est_sim_EM_gen, 2, sd),3)
#summary(beta_est_sim_num_noval)
#round(apply(beta_est_sim_num_noval, 2, sd),3)
summary(beta_est_sim_num_con)
round(apply(beta_est_sim_num_con, 2, sd),3)
#summary(beta_est_sim_num_gen)
#round(apply(beta_est_sim_num_gen, 2, sd),3)
summary(beta_est_sim_val_gen)
round(apply(beta_est_sim_val_gen, 2, sd),3)



```

PokrytÃ­ intervalÅ¯ spolehlivosti a srovnÃ¡nÃ­ dÃ©lek intervalÅ¯ spolehlivosti

```{r}

## EM
#EM_known_constant_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_EM_noval[,1] & beta_sim[1] <= beta_upper_est_sim_EM_noval[,1])
#(EM_known_constant_coverage_CI_beta1 <- EM_known_constant_covered_beta1  / num_samples)
#beta1_sim_CI_lengths_EM_known_constant<-(beta_upper_est_sim_EM_noval[,1]-beta_lower_est_sim_EM_noval[,1])
#round(mean(beta1_sim_CI_lengths_EM_known_constant),3)
#median(beta1_sim_CI_lengths_EM_known_constant)

#EM_known_constant_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_EM_noval[,2] & beta_sim[2] <= beta_upper_est_sim_EM_noval[,2])
#(EM_known_constant_coverage_CI_beta2 <- EM_known_constant_covered_beta2  / num_samples)
#beta2_sim_CI_lengths_EM_known_constant<-(beta_upper_est_sim_EM_noval[,2]-beta_lower_est_sim_EM_noval[,2])
#round(mean(beta2_sim_CI_lengths_EM_known_constant),3)
#median(beta2_sim_CI_lengths_EM_known_constant)

not_na_indexes1<-which(!is.na(beta_lower_est_sim_EM_con[,1]) & !is.na(beta_upper_est_sim_EM_con[,1]))
length(not_na_indexes1) ## pokud je problÃ©m s invertovÃ¡nÃ­m pozorovanÃ© IM
not_na_indexes2<-which(!is.na(beta_lower_est_sim_EM_con[,2]) & !is.na(beta_upper_est_sim_EM_con[,2]))
length(not_na_indexes2)

EM_unknown_constant_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_EM_con[not_na_indexes1,1] & beta_sim[1] <= beta_upper_est_sim_EM_con[not_na_indexes1,1])
(EM_unknown_constant_coverage_CI_beta1 <- EM_unknown_constant_covered_beta1  / length(not_na_indexes1))
beta1_sim_CI_lengths_EM_unknown_constant<-(beta_upper_est_sim_EM_con[not_na_indexes1,1]-beta_lower_est_sim_EM_con[not_na_indexes1,1])
round(mean(beta1_sim_CI_lengths_EM_unknown_constant),3)
median(beta1_sim_CI_lengths_EM_unknown_constant)

EM_unknown_constant_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_EM_con[not_na_indexes2,2] & beta_sim[2] <= beta_upper_est_sim_EM_con[not_na_indexes2,2])
(EM_unknown_constant_coverage_CI_beta2 <- EM_unknown_constant_covered_beta2  / length(not_na_indexes2))
beta2_sim_CI_lengths_EM_unknown_constant<-(beta_upper_est_sim_EM_con[not_na_indexes2,2]-beta_lower_est_sim_EM_con[not_na_indexes2,2])
round(mean(beta2_sim_CI_lengths_EM_unknown_constant),3)
median(beta2_sim_CI_lengths_EM_unknown_constant)



#not_na_indexes1<-which(!is.na(beta_lower_est_sim_EM_gen[,1]) & !is.na(beta_upper_est_sim_EM_gen[,1]))
#length(not_na_indexes1) ## pokud je problÃ©m s invertovÃ¡nÃ­m pozorovanÃ© IM
#not_na_indexes2<-which(!is.na(beta_lower_est_sim_EM_gen[,2]) & !is.na(beta_upper_est_sim_EM_gen[,2]))
#length(not_na_indexes2)

#EM_unknown_gen_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_EM_gen[not_na_indexes1,1] & beta_sim[1] <= beta_upper_est_sim_EM_gen[not_na_indexes1,1])
#(EM_unknown_gen_coverage_CI_beta1 <- EM_unknown_gen_covered_beta1  / length(not_na_indexes1))
#beta1_sim_CI_lengths_EM_unknown_gen<-(beta_upper_est_sim_EM_gen[not_na_indexes1,1]-beta_lower_est_sim_EM_gen[not_na_indexes1,1])
#round(mean(beta1_sim_CI_lengths_EM_unknown_gen),3)
#median(beta1_sim_CI_lengths_EM_unknown_gen)

#EM_unknown_gen_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_EM_gen[not_na_indexes2,2] & beta_sim[2] <= beta_upper_est_sim_EM_gen[not_na_indexes2,2])
#(EM_unknown_gen_coverage_CI_beta2 <- EM_unknown_gen_covered_beta2  / length(not_na_indexes2))
#beta2_sim_CI_lengths_EM_unknown_gen<-(beta_upper_est_sim_EM_gen[not_na_indexes2,2]-beta_lower_est_sim_EM_gen[not_na_indexes2,2])
#round(mean(beta2_sim_CI_lengths_EM_unknown_gen),3)
#median(beta2_sim_CI_lengths_EM_unknown_gen)


## num
#num_known_constant_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_num_noval[,1] & beta_sim[1] <= beta_upper_est_sim_num_noval[,1])
#(num_known_constant_coverage_CI_beta1 <- num_known_constant_covered_beta1  / num_samples)
#beta1_sim_CI_lengths_num_known_constant<-(beta_upper_est_sim_num_noval[,1]-beta_lower_est_sim_num_noval[,1])
#round(mean(beta1_sim_CI_lengths_num_known_constant),3)
#median(beta1_sim_CI_lengths_num_known_constant)

#num_known_constant_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_num_noval[,2] & beta_sim[2] <= beta_upper_est_sim_num_noval[,2])
#(num_known_constant_coverage_CI_beta2 <- num_known_constant_covered_beta2  / num_samples)
#beta2_sim_CI_lengths_num_known_constant<-(beta_upper_est_sim_num_noval[,2]-beta_lower_est_sim_num_noval[,2])
#round(mean(beta2_sim_CI_lengths_num_known_constant),3)
#median(beta2_sim_CI_lengths_num_known_constant)

not_na_indexes1<-which(!is.na(beta_lower_est_sim_num_con[,1]) & !is.na(beta_upper_est_sim_num_con[,1]))
length(not_na_indexes1) ## pokud je problÃ©m s invertovÃ¡nÃ­m pozorovanÃ© IM
not_na_indexes2<-which(!is.na(beta_lower_est_sim_num_con[,2]) & !is.na(beta_upper_est_sim_num_con[,2]))
length(not_na_indexes2)

num_unknown_constant_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_num_con[not_na_indexes1,1] & beta_sim[1] <= beta_upper_est_sim_num_con[not_na_indexes1,1])
(num_unknown_constant_coverage_CI_beta1 <- num_unknown_constant_covered_beta1  / length(not_na_indexes1) )
beta1_sim_CI_lengths_num_unknown_constant<-(beta_upper_est_sim_num_con[not_na_indexes1,1]-beta_lower_est_sim_num_con[not_na_indexes1,1])
round(mean(beta1_sim_CI_lengths_num_unknown_constant),3)
median(beta1_sim_CI_lengths_num_unknown_constant)

num_unknown_constant_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_num_con[not_na_indexes2,2] & beta_sim[2] <= beta_upper_est_sim_num_con[not_na_indexes2,2])
(num_unknown_constant_coverage_CI_beta2 <- num_unknown_constant_covered_beta2  / length(not_na_indexes2))
beta2_sim_CI_lengths_num_unknown_constant<-(beta_upper_est_sim_num_con[not_na_indexes2,2]-beta_lower_est_sim_num_con[not_na_indexes2,2])
round(mean(beta2_sim_CI_lengths_num_unknown_constant),3)
median(beta2_sim_CI_lengths_num_unknown_constant)

#not_na_indexes1<-which(!is.na(beta_lower_est_sim_num_gen[,1]) & !is.na(beta_upper_est_sim_num_gen[,1]))
#length(not_na_indexes1) ## pokud je problÃ©m s invertovÃ¡nÃ­m pozorovanÃ© IM
#not_na_indexes2<-which(!is.na(beta_lower_est_sim_num_gen[,2]) & !is.na(beta_upper_est_sim_num_gen[,2]))
#length(not_na_indexes2)

#num_unknown_gen_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_num_gen[not_na_indexes1,1] & beta_sim[1] <= beta_upper_est_sim_num_gen[not_na_indexes1,1])
#(num_unknown_gen_coverage_CI_beta1 <- num_unknown_gen_covered_beta1  / length(not_na_indexes1))
#beta1_sim_CI_lengths_num_unknown_gen<-(beta_upper_est_sim_num_gen[not_na_indexes1,1]-beta_lower_est_sim_num_gen[not_na_indexes1,1])
#round(mean(beta1_sim_CI_lengths_num_unknown_gen),3)
#median(beta1_sim_CI_lengths_num_unknown_gen)

#num_unknown_gen_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_num_gen[not_na_indexes2,2] & beta_sim[2] <= beta_upper_est_sim_num_gen[not_na_indexes2,2])
#(num_unknown_gen_coverage_CI_beta2 <- num_unknown_gen_covered_beta2  / length(not_na_indexes2))
#beta2_sim_CI_lengths_num_unknown_gen<-(beta_upper_est_sim_num_gen[not_na_indexes2,2]-beta_lower_est_sim_num_gen[not_na_indexes2,2])
#round(mean(beta2_sim_CI_lengths_num_unknown_gen),3)
#median(beta2_sim_CI_lengths_num_unknown_gen)



## val
val_unknown_covered_beta1 <- sum(beta_sim[1] >= beta_lower_est_sim_val_gen[,1] & beta_sim[1] <= beta_upper_est_sim_val_gen[,1])
(val_unknown_coverage_CI_beta1 <- val_unknown_covered_beta1  / num_samples)
beta1_sim_CI_lengths_val_unknown<-(beta_upper_est_sim_val_gen[,1]-beta_lower_est_sim_val_gen[,1])
round(mean(beta1_sim_CI_lengths_val_unknown),3)
median(beta1_sim_CI_lengths_val_unknown)

val_unknown_covered_beta2 <- sum(beta_sim[2] >= beta_lower_est_sim_val_gen[,2] & beta_sim[2] <= beta_upper_est_sim_val_gen[,2])
(val_unknown_coverage_CI_beta2 <- val_unknown_covered_beta2  / num_samples)
beta2_sim_CI_lengths_val_unknown<-(beta_upper_est_sim_val_gen[,2]-beta_lower_est_sim_val_gen[,2])
round(mean(beta2_sim_CI_lengths_val_unknown),3)
median(beta2_sim_CI_lengths_val_unknown)
```


grafickÃ© znÃ¡zornÄnÃ­ jÃ¡drovÃ½ch odhadÅ¯ hustot beta1 a beta2 dle jednotlivÃ½ch metod

```{r}
par(mfrow=c(1, 2))

#dens1_em <- density(beta_est_sim_EM_noval$beta1)
dens1_em <- density(beta_est_sim_EM_con$beta1)
#dens1_em <- density(beta_est_sim_EM_gen$beta1)
plot(dens1_em, main="JÃ¡drovÃ½ odhad hustoty beta1 \n w=1/4, n=1000", xlab="beta1", ylab = "hustota", col="blue", xlim=c(-0.6, 1.4))
#dens1_num <- density(beta_est_sim_num_noval$beta1)
dens1_num <- density(beta_est_sim_num_con$beta1)
#dens1_num <- density(beta_est_sim_num_gen$beta1)
lines(dens1_num, col="green")
dens1_pure_val <- density(beta_est_sim_val_gen$beta1)
lines(dens1_pure_val, col="red")


legend("topleft", legend=c("EM", "Num", "Val"), col=c("blue", "green","red"), lty=1, cex=0.8)

#dens2_em <- density(beta_est_sim_EM_noval$beta2)
dens2_em <- density(beta_est_sim_EM_con$beta2)
#dens2_em <- density(beta_est_sim_EM_gen$beta2)
plot(dens2_em, main="JÃ¡drovÃ½ odhad hustoty beta2 \n w=1/4, n=1000", ylab = "hustota", xlab="beta2", col="blue", xlim=c(-0.2, 1.8))
#dens2_num <- density(beta_est_sim_num_noval$beta2)
dens2_num <- density(beta_est_sim_num_con$beta2)
#dens2_num <- density(beta_est_sim_num_gen$beta2)
lines(dens2_num, col="green")
dens2_pure_val <- density(beta_est_sim_val_gen$beta2)
lines(dens2_pure_val, col="red")


legend("topleft", legend=c("EM", "Num", "Val"), col=c("blue", "green", "red"), lty=1, cex=0.8)

plot_gen_1000_w4 <- recordPlot()

## replayPlot(plot_noval_500)
```

```{r}
pdf("pr2.2.8.250.pdf", width = 12, height = 8)

plot_gen_250_w8

dev.off()

pdf("pr2.2.4.250.pdf", width = 12, height = 8)

plot_gen_250_w4

dev.off()

pdf("pr2.2.8.1000.pdf", width = 12, height = 8)

plot_gen_1000_w8

dev.off()

pdf("pr2.2.4.1000.pdf", width = 12, height = 8)

plot_gen_1000_w4

dev.off()


```



## PÅÃ­klad 3

Zopakujeme postup z pÅÃ­kladu 2.3, jen tentokrÃ¡t s vÃ­ce vysvÄtlujÃ­cÃ­mi promÄnnÃ½mi a senzitivitÄ a specificitÄ na nich aktivnÄ zÃ¡visejÃ­cÃ­mi. Vezmeme kromÄ interceptu tÅi vysvÄtlujÃ­cÃ­ promÄnnÃ© postupnÄ normÃ¡lnÃ­, rovnomÄrnÃ© a alternativnÃ­ rozdÄlenÃ­. GenerovanÃ© s korelaÄnÃ­ strukturou vychÃ¡zejÃ­cÃ­ z dÅÃ­ve definovanÃ© funkce generate_correlated_dataframe.

```{r}
start_time <- Sys.time()

set.seed(1234)
mu_ex3<-0
var_ex3<-1
l_ex3<-0
u_ex3<-2
p_ex3<-0.5
ro12<-0.4
ro13<-0.3
ro23<-0.4

#is.positive.definite(matrix(c(1, ro12, ro13, ro12, 1, ro23, ro13, ro23, 1), nrow = 3, ncol = 3, byrow = TRUE))

n_ex3<-500 ## rozsah vÃ½bÄru n
num_samples_ex3<-1000 ## generovÃ¡nÃ­ odezvy a chyby v odezvÄ pro jednu realizaci vysvÄtlujÃ­cÃ­ch promÄnnÃ½ch. Z toho odhady skuteÄnÃ©ho parametru, 1000 realizacÃ­ provedeme


beta_ex3 <- c(-0.4, 0.1, 0.2, 0.6) ## skuteÄnÃ¡ beta
gamma_0_ex3 <- c(1, 0.1, 0.2, 0.5) ## skuteÄnÃ© gamma
gamma_1_ex3 <- c(-1, -0.1, -0.3, -0.4)


valid_prob_ex3<- 1/4 ## pst vÃ½bÄru do validaÄnÃ­ skupiny w

## pÅedpÅipravenÃ­ dataframes pro vyplÅovÃ¡nÃ­ bÄhem simulace

## EM algoritmus

beta_est_ex3_EM_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
beta_lower_est_ex3_EM_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
beta_upper_est_ex3_EM_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
#gamma0_est_ex3_EM_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma0_lower_est_ex3_EM_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma0_upper_est_ex3_EM_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma1_est_ex3_EM_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())
#gamma1_lower_est_ex3_EM_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())
#gamma1_upper_est_ex3_EM_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())

## numerickÃ¡ metoda

beta_est_ex3_num_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
beta_lower_est_ex3_num_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
beta_upper_est_ex3_num_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
#gamma0_est_ex3_num_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma0_lower_est_ex3_num_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma0_upper_est_ex3_num_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma1_est_ex3_num_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())
#gamma1_lower_est_ex3_num_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())
#gamma1_upper_est_ex3_num_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())


## pouze validaÄnÃ­ skupina, tedy bÄÅ¾nÃ¡ logistickÃ¡ regrese ÄistÄ na validaÄnÃ­ skupinu, zbytek zahodit

beta_est_ex3_val_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
beta_lower_est_ex3_val_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
beta_upper_est_ex3_val_gen <- data.frame(beta1 = numeric(), beta2 = numeric(), beta3 = numeric(), beta4 = numeric())
#gamma0_est_ex3_val_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma0_lower_est_ex3_val_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma0_upper_est_ex3_val_gen <- data.frame(gamma0_1 = numeric(), gamma0_2 = numeric(), gamma0_3 = numeric(), gamma0_4 = numeric())
#gamma1_est_ex3_val_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())
#gamma1_lower_est_ex3_val_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())
#gamma1_upper_est_ex3_val_gen <- data.frame(gamma1_1 = numeric(), gamma1_2 = numeric(), gamma1_3 = numeric(), gamma1_4 = numeric())

win_EM<-0
win_num<-0

for (i in 1:num_samples_ex3) { ## kvÅ¯li rychlosti kÃ³d spouÅ¡tÄnÃ½ se zakomentovanÃ½mi zrovna nepouÅ¾itÃ½mi ÄÃ¡stmi
  covariates_ex3 <- generate_correlated_dataframe(n_ex3, mu_ex3, var_ex3, l_ex3, u_ex3, p_ex3, ro12, ro13, ro23, 1)
  
  response_ex3 <- generate_response(beta_ex3, covariates_ex3)
  
  response_with_mistakes_ex3 <- generate_response_with_mistakes(gamma_0_ex3,gamma_1_ex3,covariates_ex3, response_ex3)
  
  
  beta0_ex3<-runif(4, min = -1, max = 1) ## pro rÅ¯znÃ© realizace mÅ¯Å¾eme zaÄÃ­nat z rÅ¯znÃ½ch nÃ¡hodnÃ½ch poÄÃ¡teÄnÃ­ch odhadÅ¯
  gamma_0_0_ex3<-runif(4, min = 0, max = 1)
  gamma_1_0_ex3<-runif(4, min = -1, max = 0)
  
  val_ind_ex3<-rbinom(n_ex3, 1, valid_prob_ex3) ## pÅiÅazenÃ­ hodnoty 1 prvkÅ¯m validaÄnÃ­ skupiny (prvnÃ­ch k prvkÅ¯ odpovÃ­dajÃ­cÃ­ vÃ½Å¡e zmÃ­nÄnÃ©mu podÃ­lu rozsahu) a hodnoty 0 prvkÅ¯m nevalidaÄnÃ­ skupiny

    
  
  

  ## EM algoritmus ÄÃ¡st
  ## EM algoritmus s sens a spec zÃ¡visejÃ­cÃ­ na vysv. promÄnnÃ½ch

  tryCatch({par_est_ex3_EM_gen<-EM_general(response_with_mistakes_ex3, response_ex3[val_ind_ex3==1], val_ind_ex3, covariates_ex3, beta0_ex3, gamma_0_0_ex3, gamma_1_0_ex3)
  }, error = function(e) {
    cat("Error on iteration", i, ": ", e$message, "\n")
  })

  
  
  beta_est_ex3_EM_gen <- rbind(beta_est_ex3_EM_gen, data.frame(beta1 = par_est_ex3_EM_gen[1, 1], beta2 = par_est_ex3_EM_gen[1, 2], beta3 = par_est_ex3_EM_gen[1, 3], beta4 = par_est_ex3_EM_gen[1, 4]))
  beta_lower_est_ex3_EM_gen <- rbind(beta_lower_est_ex3_EM_gen, data.frame(beta1 = par_est_ex3_EM_gen[2, 1], beta2 = par_est_ex3_EM_gen[2, 2], beta3 = par_est_ex3_EM_gen[2, 3], beta4 = par_est_ex3_EM_gen[2, 4]))
  beta_upper_est_ex3_EM_gen <- rbind(beta_upper_est_ex3_EM_gen, data.frame(beta1 = par_est_ex3_EM_gen[3, 1], beta2 = par_est_ex3_EM_gen[3, 2], beta3 = par_est_ex3_EM_gen[3, 3], beta4 = par_est_ex3_EM_gen[3, 4]))
  
  #gamma0_est_ex3_EM_gen <- rbind(gamma0_est_ex3_EM_gen, data.frame(gamma0_1 = par_est_ex3_EM_gen[4, 1], gamma0_2 = par_est_ex3_EM_gen[4, 2], gamma0_3 = par_est_ex3_EM_gen[4, 3], gamma0_4 = par_est_ex3_EM_gen[4, 4]))
  #gamma0_lower_est_ex3_EM_gen <- rbind(gamma0_lower_est_ex3_EM_gen, data.frame(gamma0_1 = par_est_ex3_EM_gen[5, 1], gamma0_2 = par_est_ex3_EM_gen[5, 2], gamma0_3 = par_est_ex3_EM_gen[5, 3], gamma0_4 = par_est_ex3_EM_gen[5, 4]))
  #gamma0_upper_est_ex3_EM_gen <- rbind(gamma0_upper_est_ex3_EM_gen, data.frame(gamma0_1 = par_est_ex3_EM_gen[6, 1], gamma0_2 = par_est_ex3_EM_gen[6, 2], gamma0_3 = par_est_ex3_EM_gen[6, 3], gamma0_4 = par_est_ex3_EM_gen[6, 4]))
  
  
  #gamma1_est_ex3_EM_gen <- rbind(gamma1_est_ex3_EM_gen, data.frame(gamma1_1 = par_est_ex3_EM_gen[7, 1], gamma1_2 = par_est_ex3_EM_gen[7, 2], gamma1_3 = par_est_ex3_EM_gen[7, 3], gamma1_4 = par_est_ex3_EM_gen[7, 4]))
  #gamma1_lower_est_ex3_EM_gen <- rbind(gamma1_lower_est_ex3_EM_gen, data.frame(gamma1_1 = par_est_ex3_EM_gen[8, 1], gamma1_2 = par_est_ex3_EM_gen[8, 2], gamma1_3 = par_est_ex3_EM_gen[8, 3], gamma1_4 = par_est_ex3_EM_gen[8, 4]))
  #gamma1_upper_est_ex3_EM_gen <- rbind(gamma1_upper_est_ex3_EM_gen, data.frame(gamma1_1 = par_est_ex3_EM_gen[9, 1], gamma1_2 = par_est_ex3_EM_gen[9, 2], gamma1_3 = par_est_ex3_EM_gen[9, 3], gamma1_4 = par_est_ex3_EM_gen[9, 4]))


## numerickÃ¡ maximalizace ÄÃ¡st
  ## num algoritmus s sens a spec zÃ¡visejÃ­cÃ­ na vysv. promÄnnÃ½ch

  tryCatch({par_est_ex3_num_gen<-suppressWarnings(numeric_max(response_with_mistakes_ex3, response_ex3[val_ind_ex3==1], val_ind_ex3, covariates_ex3, beta0_ex3, gamma_0_0_ex3, gamma_1_0_ex3))
  }, error = function(e) {
    cat("Error on iteration", i, ": ", e$message, "\n")
  })

  
  
  beta_est_ex3_num_gen <- rbind(beta_est_ex3_num_gen, data.frame(beta1 = par_est_ex3_num_gen[1, 1], beta2 = par_est_ex3_num_gen[1, 2], beta3 = par_est_ex3_num_gen[1, 3], beta4 = par_est_ex3_num_gen[1, 4]))
  beta_lower_est_ex3_num_gen <- rbind(beta_lower_est_ex3_num_gen, data.frame(beta1 = par_est_ex3_num_gen[2, 1], beta2 = par_est_ex3_num_gen[2, 2], beta3 = par_est_ex3_num_gen[2, 3], beta4 = par_est_ex3_num_gen[2, 4]))
  beta_upper_est_ex3_num_gen <- rbind(beta_upper_est_ex3_num_gen, data.frame(beta1 = par_est_ex3_num_gen[3, 1], beta2 = par_est_ex3_num_gen[3, 2], beta3 = par_est_ex3_num_gen[3, 3], beta4 = par_est_ex3_num_gen[3, 4]))
  
  #gamma0_est_ex3_num_gen <- rbind(gamma0_est_ex3_num_gen, data.frame(gamma0_1 = par_est_ex3_num_gen[4, 1], gamma0_2 = par_est_ex3_num_gen[4, 2], gamma0_3 = par_est_ex3_num_gen[4, 3], gamma0_4 = par_est_ex3_num_gen[4, 4]))
  #gamma0_lower_est_ex3_num_gen <- rbind(gamma0_lower_est_ex3_num_gen, data.frame(gamma0_1 = par_est_ex3_num_gen[5, 1], gamma0_2 = par_est_ex3_num_gen[5, 2], gamma0_3 = par_est_ex3_num_gen[5, 3], gamma0_4 = par_est_ex3_num_gen[5, 4]))
  #gamma0_upper_est_ex3_num_gen <- rbind(gamma0_upper_est_ex3_num_gen, data.frame(gamma0_1 = par_est_ex3_num_gen[6, 1], gamma0_2 = par_est_ex3_num_gen[6, 2], gamma0_3 = par_est_ex3_num_gen[6, 3], gamma0_4 = par_est_ex3_num_gen[6, 4]))
  
  
  #gamma1_est_ex3_num_gen <- rbind(gamma1_est_ex3_num_gen, data.frame(gamma1_1 = par_est_ex3_num_gen[7, 1], gamma1_2 = par_est_ex3_num_gen[7, 2], gamma1_3 = par_est_ex3_num_gen[7, 3], gamma1_4 = par_est_ex3_num_gen[7, 4]))
  #gamma1_lower_est_ex3_num_gen <- rbind(gamma1_lower_est_ex3_num_gen, data.frame(gamma1_1 = par_est_ex3_num_gen[8, 1], gamma1_2 = par_est_ex3_num_gen[8, 2], gamma1_3 = par_est_ex3_num_gen[8, 3], gamma1_4 = par_est_ex3_num_gen[8, 4]))
  #gamma1_upper_est_ex3_num_gen <- rbind(gamma1_upper_est_ex3_num_gen, data.frame(gamma1_1 = par_est_ex3_num_gen[9, 1], gamma1_2 = par_est_ex3_num_gen[9, 2], gamma1_3 = par_est_ex3_num_gen[9, 3], gamma1_4 = par_est_ex3_num_gen[9, 4]))
  


## bÄÅ¾nÃ¡ logistickÃ¡ regrese na validaÄnÃ­ skupinu, princip je udÄlat si z toho jen logistickou regresei na validaÄnÃ­ skupinÄ, proto je vlastnÄ zbyteÄnÄ navÃ­c poÄÃ­tat sens a spec, proto zakomentovanÃ©
  ## pouze validaÄnÃ­ vÃ½bÄr pÅi sens a spec zÃ¡visejÃ­cÃ­ na vysv. promÄnnÃ½ch
  indexes1_ex3<-which(response_ex3==1 & val_ind_ex3==1)
  indexes0_ex3<-which(response_ex3==0 & val_ind_ex3==1)
  
  glm_beta<-summary(glm(response_ex3[val_ind_ex3==1]~.-1, family = binomial(link = "logit"), data = covariates_ex3[val_ind_ex3==1,]))$coefficients

  #glm_gamma0<-summary(glm(response_with_mistakes_ex3[indexes1_ex3]~.-1, family = binomial(link = "logit"), data = covariates_ex3[indexes1_ex3,]))$coefficients
  #glm_gamma1<-summary(glm(response_with_mistakes_ex3[indexes0_ex3]~.-1, family = binomial(link = "logit"), data = covariates_ex3[indexes0_ex3,]))$coefficients
  
  
  
  par_est_ex3_beta<-glm_beta[,"Estimate"]
  par_est_ex3_beta_lower<-par_est_ex3_beta-glm_beta[,"Std. Error"]*qnorm(0.975)
  par_est_ex3_beta_upper<-par_est_ex3_beta+glm_beta[,"Std. Error"]*qnorm(0.975)
  #par_est_ex3_gamma0<-glm_gamma0[,"Estimate"]
  #par_est_ex3_gamma0_lower<-par_est_ex3_gamma0-glm_gamma0[,"Std. Error"]*qnorm(0.975)
  #par_est_ex3_gamma0_upper<-par_est_ex3_gamma0+glm_gamma0[,"Std. Error"]*qnorm(0.975)
  #par_est_ex3_gamma1<-glm_gamma1[,"Estimate"]
  #par_est_ex3_gamma1_lower<-par_est_ex3_gamma1-glm_gamma1[,"Std. Error"]*qnorm(0.975)
  #par_est_ex3_gamma1_upper<-par_est_ex3_gamma1+glm_gamma1[,"Std. Error"]*qnorm(0.975)
  
  beta_est_ex3_val_gen<-rbind(beta_est_ex3_val_gen, data.frame(beta1= par_est_ex3_beta[1], beta2 = par_est_ex3_beta[2], beta3 = par_est_ex3_beta[3], beta4 = par_est_ex3_beta[4]))
  beta_lower_est_ex3_val_gen<-rbind(beta_lower_est_ex3_val_gen, data.frame(beta1= par_est_ex3_beta_lower[1], beta2 = par_est_ex3_beta_lower[2], beta3 = par_est_ex3_beta_lower[3], beta4 = par_est_ex3_beta_lower[4]))
  beta_upper_est_ex3_val_gen<-rbind(beta_upper_est_ex3_val_gen, data.frame(beta1= par_est_ex3_beta_upper[1], beta2 = par_est_ex3_beta_upper[2], beta3 = par_est_ex3_beta_upper[3], beta4 = par_est_ex3_beta_upper[4]))
  #gamma0_est_ex3_val_gen<-rbind(gamma0_est_ex3_val_gen, data.frame(gamma0_1= par_est_ex3_gamma0[1], gamma0_2 = par_est_ex3_gamma0[2], gamma0_3 = par_est_ex3_gamma0[3], gamma0_4 = par_est_ex3_gamma0[4]))
  #gamma0_lower_est_ex3_val_gen<-rbind(gamma0_lower_est_ex3_val_gen, data.frame(gamma0_1= par_est_ex3_gamma0_lower[1], gamma0_2 = par_est_ex3_gamma0_lower[2], gamma0_3 = par_est_ex3_gamma0_lower[3], gamma0_4 = par_est_ex3_gamma0_lower[4]))
  #gamma0_upper_est_ex3_val_gen<-rbind(gamma0_upper_est_ex3_val_gen, data.frame(gamma0_1= par_est_ex3_gamma0_upper[1], gamma0_2 = par_est_ex3_gamma0_upper[2], gamma0_3 = par_est_ex3_gamma0_upper[3], gamma0_4 = par_est_ex3_gamma0_upper[4]))
  #gamma1_est_ex3_val_gen<-rbind(gamma1_est_ex3_val_gen, data.frame(gamma1_1= par_est_ex3_gamma1[1], gamma1_2 = par_est_ex3_gamma1[2], gamma1_3 = par_est_ex3_gamma1[3], gamma1_4 = par_est_ex3_gamma1[4]))
  #gamma1_lower_est_ex3_val_gen<-rbind(gamma1_lower_est_ex3_val_gen, data.frame(gamma1_1= par_est_ex3_gamma1_lower[1], gamma1_2 = par_est_ex3_gamma1_lower[2], gamma1_3 = par_est_ex3_gamma1_lower[3], gamma1_4 = par_est_ex3_gamma1_lower[4]))
  #gamma1_upper_est_ex3_val_gen<-rbind(gamma1_upper_est_ex3_val_gen, data.frame(gamma1_1= par_est_ex3_gamma1_upper[1], gamma1_2 = par_est_ex3_gamma1_upper[2], gamma1_3 = par_est_ex3_gamma1_upper[3], gamma1_4 = par_est_ex3_gamma1_upper[4]))
}





end_time <- Sys.time()

as.numeric(end_time - start_time, units = "secs")


```


```{r}
beta_ex3

summary(beta_est_ex3_EM_gen)
round(apply(beta_est_ex3_EM_gen, 2, sd),3)

summary(beta_est_ex3_num_gen)
round(apply(beta_est_ex3_num_gen, 2, sd),3)

summary(beta_est_ex3_val_gen)
round(apply(beta_est_ex3_val_gen, 2, sd),3)



```


```{r}

## EM
print("EM")
not_na_indexes1<-which(!is.na(beta_lower_est_ex3_EM_gen[,1]) & !is.na(beta_upper_est_ex3_EM_gen[,1]))
length(not_na_indexes1) ## pokud je problÃ©m s invertovÃ¡nÃ­m pozorovanÃ© IM
not_na_indexes2<-which(!is.na(beta_lower_est_ex3_EM_gen[,2]) & !is.na(beta_upper_est_ex3_EM_gen[,2]))
length(not_na_indexes2)
not_na_indexes3<-which(!is.na(beta_lower_est_ex3_EM_gen[,3]) & !is.na(beta_upper_est_ex3_EM_gen[,3]))
length(not_na_indexes3)
not_na_indexes4<-which(!is.na(beta_lower_est_ex3_EM_gen[,4]) & !is.na(beta_upper_est_ex3_EM_gen[,4]))
length(not_na_indexes4)

print("beta1")
EM_unknown_gen_covered_beta1 <- sum(beta_ex3[1] >= beta_lower_est_ex3_EM_gen[not_na_indexes1,1] & beta_ex3[1] <= beta_upper_est_ex3_EM_gen[not_na_indexes1,1])
(EM_unknown_gen_coverage_CI_beta1 <- EM_unknown_gen_covered_beta1  / length(not_na_indexes1))
beta1_ex3_CI_lengths_EM_unknown_gen<-(beta_upper_est_ex3_EM_gen[not_na_indexes1,1]-beta_lower_est_ex3_EM_gen[not_na_indexes1,1])
round(mean(beta1_ex3_CI_lengths_EM_unknown_gen),3)
median(beta1_ex3_CI_lengths_EM_unknown_gen)

print("beta2")
EM_unknown_gen_covered_beta2 <- sum(beta_ex3[2] >= beta_lower_est_ex3_EM_gen[not_na_indexes2,2] & beta_ex3[2] <= beta_upper_est_ex3_EM_gen[not_na_indexes2,2])
(EM_unknown_gen_coverage_CI_beta2 <- EM_unknown_gen_covered_beta2  / length(not_na_indexes2))
beta2_ex3_CI_lengths_EM_unknown_gen<-(beta_upper_est_ex3_EM_gen[not_na_indexes2,2]-beta_lower_est_ex3_EM_gen[not_na_indexes2,2])
round(mean(beta2_ex3_CI_lengths_EM_unknown_gen),3)
median(beta2_ex3_CI_lengths_EM_unknown_gen)

print("beta3")
EM_unknown_gen_covered_beta3 <- sum(beta_ex3[3] >= beta_lower_est_ex3_EM_gen[not_na_indexes3,3] & beta_ex3[3] <= beta_upper_est_ex3_EM_gen[not_na_indexes3,3])
(EM_unknown_gen_coverage_CI_beta3 <- EM_unknown_gen_covered_beta3  / length(not_na_indexes3))
beta3_ex3_CI_lengths_EM_unknown_gen<-(beta_upper_est_ex3_EM_gen[not_na_indexes3,3]-beta_lower_est_ex3_EM_gen[not_na_indexes3,3])
round(mean(beta3_ex3_CI_lengths_EM_unknown_gen),3)
median(beta3_ex3_CI_lengths_EM_unknown_gen)

print("beta4")
EM_unknown_gen_covered_beta4 <- sum(beta_ex3[4] >= beta_lower_est_ex3_EM_gen[not_na_indexes4,4] & beta_ex3[4] <= beta_upper_est_ex3_EM_gen[not_na_indexes4,4])
(EM_unknown_gen_coverage_CI_beta4 <- EM_unknown_gen_covered_beta4  / length(not_na_indexes4))
beta4_ex3_CI_lengths_EM_unknown_gen<-(beta_upper_est_ex3_EM_gen[not_na_indexes4,4]-beta_lower_est_ex3_EM_gen[not_na_indexes4,4])
round(mean(beta4_ex3_CI_lengths_EM_unknown_gen),3)
median(beta4_ex3_CI_lengths_EM_unknown_gen)


## num
print("num")
not_na_indexes1<-which(!is.na(beta_lower_est_ex3_num_gen[,1]) & !is.na(beta_upper_est_ex3_num_gen[,1]))
length(not_na_indexes1) ## pokud je problÃ©m s invertovÃ¡nÃ­m pozorovanÃ© IM
not_na_indexes2<-which(!is.na(beta_lower_est_ex3_num_gen[,2]) & !is.na(beta_upper_est_ex3_num_gen[,2]))
length(not_na_indexes2)
not_na_indexes3<-which(!is.na(beta_lower_est_ex3_num_gen[,3]) & !is.na(beta_upper_est_ex3_num_gen[,3]))
length(not_na_indexes3)
not_na_indexes4<-which(!is.na(beta_lower_est_ex3_num_gen[,4]) & !is.na(beta_upper_est_ex3_num_gen[,4]))
length(not_na_indexes4)

print("beta1")
num_unknown_gen_covered_beta1 <- sum(beta_ex3[1] >= beta_lower_est_ex3_num_gen[not_na_indexes1,1] & beta_ex3[1] <= beta_upper_est_ex3_num_gen[not_na_indexes1,1])
(num_unknown_gen_coverage_CI_beta1 <- num_unknown_gen_covered_beta1  / length(not_na_indexes1))
beta1_ex3_CI_lengths_num_unknown_gen<-(beta_upper_est_ex3_num_gen[not_na_indexes1,1]-beta_lower_est_ex3_num_gen[not_na_indexes1,1])
round(mean(beta1_ex3_CI_lengths_num_unknown_gen),3)
median(beta1_ex3_CI_lengths_num_unknown_gen)

print("beta2")
num_unknown_gen_covered_beta2 <- sum(beta_ex3[2] >= beta_lower_est_ex3_num_gen[not_na_indexes2,2] & beta_ex3[2] <= beta_upper_est_ex3_num_gen[not_na_indexes2,2])
(num_unknown_gen_coverage_CI_beta2 <- num_unknown_gen_covered_beta2  / length(not_na_indexes2))
beta2_ex3_CI_lengths_num_unknown_gen<-(beta_upper_est_ex3_num_gen[not_na_indexes2,2]-beta_lower_est_ex3_num_gen[not_na_indexes2,2])
round(mean(beta2_ex3_CI_lengths_num_unknown_gen),3)
median(beta2_ex3_CI_lengths_num_unknown_gen)

print("beta3")
num_unknown_gen_covered_beta3 <- sum(beta_ex3[3] >= beta_lower_est_ex3_num_gen[not_na_indexes3,3] & beta_ex3[3] <= beta_upper_est_ex3_num_gen[not_na_indexes3,3])
(num_unknown_gen_coverage_CI_beta3 <- num_unknown_gen_covered_beta3  / length(not_na_indexes3))
beta3_ex3_CI_lengths_num_unknown_gen<-(beta_upper_est_ex3_num_gen[not_na_indexes3,3]-beta_lower_est_ex3_num_gen[not_na_indexes3,3])
round(mean(beta3_ex3_CI_lengths_num_unknown_gen),3)
median(beta3_ex3_CI_lengths_num_unknown_gen)

print("beta4")
num_unknown_gen_covered_beta4 <- sum(beta_ex3[4] >= beta_lower_est_ex3_num_gen[not_na_indexes4,4] & beta_ex3[4] <= beta_upper_est_ex3_num_gen[not_na_indexes4,4])
(num_unknown_gen_coverage_CI_beta4 <- num_unknown_gen_covered_beta4  / length(not_na_indexes4))
beta4_ex3_CI_lengths_num_unknown_gen<-(beta_upper_est_ex3_num_gen[not_na_indexes4,4]-beta_lower_est_ex3_num_gen[not_na_indexes4,4])
round(mean(beta4_ex3_CI_lengths_num_unknown_gen),3)
median(beta4_ex3_CI_lengths_num_unknown_gen)


## val
print("val")

print("beta1")
val_unknown_covered_beta1 <- sum(beta_ex3[1] >= beta_lower_est_ex3_val_gen[,1] & beta_ex3[1] <= beta_upper_est_ex3_val_gen[,1])
(val_unknown_coverage_CI_beta1 <- val_unknown_covered_beta1  / num_samples_ex3)
beta1_ex3_CI_lengths_val_unknown<-(beta_upper_est_ex3_val_gen[,1]-beta_lower_est_ex3_val_gen[,1])
round(mean(beta1_ex3_CI_lengths_val_unknown),3)
median(beta1_ex3_CI_lengths_val_unknown)

print("beta2")
val_unknown_covered_beta2 <- sum(beta_ex3[2] >= beta_lower_est_ex3_val_gen[,2] & beta_ex3[2] <= beta_upper_est_ex3_val_gen[,2])
(val_unknown_coverage_CI_beta2 <- val_unknown_covered_beta2  / num_samples_ex3)
beta2_ex3_CI_lengths_val_unknown<-(beta_upper_est_ex3_val_gen[,2]-beta_lower_est_ex3_val_gen[,2])
round(mean(beta2_ex3_CI_lengths_val_unknown),3)
median(beta2_ex3_CI_lengths_val_unknown)

print("beta3")
val_unknown_covered_beta3 <- sum(beta_ex3[3] >= beta_lower_est_ex3_val_gen[,3] & beta_ex3[3] <= beta_upper_est_ex3_val_gen[,3])
(val_unknown_coverage_CI_beta3 <- val_unknown_covered_beta3  / num_samples_ex3)
beta3_ex3_CI_lengths_val_unknown<-(beta_upper_est_ex3_val_gen[,3]-beta_lower_est_ex3_val_gen[,3])
round(mean(beta3_ex3_CI_lengths_val_unknown),3)
median(beta3_ex3_CI_lengths_val_unknown)

print("beta4")
val_unknown_covered_beta4 <- sum(beta_ex3[4] >= beta_lower_est_ex3_val_gen[,4] & beta_ex3[4] <= beta_upper_est_ex3_val_gen[,4])
(val_unknown_coverage_CI_beta4 <- val_unknown_covered_beta4  / num_samples_ex3)
beta4_ex3_CI_lengths_val_unknown<-(beta_upper_est_ex3_val_gen[,4]-beta_lower_est_ex3_val_gen[,4])
round(mean(beta4_ex3_CI_lengths_val_unknown),3)
median(beta4_ex3_CI_lengths_val_unknown)
```

grafickÃ© znÃ¡zornÄnÃ­ jÃ¡drovÃ½ch odhadÅ¯ hustot beta1 a beta2 dle jednotlivÃ½ch metod

```{r}
par(mfrow=c(2, 2))


dens1_em <- density(beta_est_ex3_EM_gen$beta1)
plot(dens1_em, main="JÃ¡drovÃ½ odhad hustoty beta1 \n w=1/4, n=500", xlab="beta1", ylab = "hustota", col="blue", xlim=c(-1.4, 0.6))
dens1_num <- density(beta_est_ex3_num_gen$beta1)
lines(dens1_num, col="green")
dens1_pure_val <- density(beta_est_ex3_val_gen$beta1)
lines(dens1_pure_val, col="red")


legend("topleft", legend=c("EM", "Num", "Val"), col=c("blue", "green","red"), lty=1, cex=0.8)

dens2_em <- density(beta_est_ex3_EM_gen$beta2)
plot(dens2_em, main="JÃ¡drovÃ½ odhad hustoty beta2 \n w=1/4, n=500", ylab = "hustota", xlab="beta2", col="blue", xlim=c(-0.9, 1.1))
dens2_num <- density(beta_est_ex3_num_gen$beta2)
lines(dens2_num, col="green")
dens2_pure_val <- density(beta_est_ex3_val_gen$beta2)
lines(dens2_pure_val, col="red")


legend("topleft", legend=c("EM", "Num", "Val"), col=c("blue", "green", "red"), lty=1, cex=0.8)

dens3_em <- density(beta_est_ex3_EM_gen$beta3)
plot(dens3_em, main="JÃ¡drovÃ½ odhad hustoty beta3 \n w=1/4, n=500", ylab = "hustota", xlab="beta3", col="blue", xlim=c(-0.8, 1.2), ylim=c(0,1.4))
dens3_num <- density(beta_est_ex3_num_gen$beta3)
lines(dens3_num, col="green")
dens3_pure_val <- density(beta_est_ex3_val_gen$beta3)
lines(dens3_pure_val, col="red")


legend("topleft", legend=c("EM", "Num", "Val"), col=c("blue", "green", "red"), lty=1, cex=0.8)

dens4_em <- density(beta_est_ex3_EM_gen$beta4)
plot(dens4_em, main="JÃ¡drovÃ½ odhad hustoty beta4 \n w=1/4, n=500", ylab = "hustota", xlab="beta4", col="blue", xlim=c(-0.4, 1.6))
dens4_num <- density(beta_est_ex3_num_gen$beta4)
lines(dens4_num, col="green")
dens4_pure_val <- density(beta_est_ex3_val_gen$beta4)
lines(dens4_pure_val, col="red")


legend("topleft", legend=c("EM", "Num", "Val"), col=c("blue", "green", "red"), lty=1, cex=0.8)

plot_gen_500_w4 <- recordPlot()

## replayPlot(plot_noval_500)
```

```{r}
pdf("pr3.8.500.pdf", width = 12, height = 8)

plot_gen_500_w8

dev.off()

pdf("pr3.4.500.pdf", width = 12, height = 8)

plot_gen_500_w4

dev.off()

pdf("pr3.8.1000.pdf", width = 12, height = 8)

plot_gen_1000_w8

dev.off()

pdf("pr3.4.1000.pdf", width = 12, height = 8)

plot_gen_1000_w4

dev.off()


```