Disagreements

The purpose of this script is to test out different “disagreements” measure.

1 Functions

1.1 Data Generating Process

### Data Generating Process ###
dgp <- function(N, periods_type, s_epsilon, s_alpha, s_beta, p_num, rho_1, rho_2, y_initial, seed){
  x = seed
  
  ## Generate the graph ##
  p_num_value <- if(p_num == "low"){
    N * .2
  } else if(p_num == "high"){
    N * .8
  }
  
  p = p_num_value/N
  
  set.seed(6*x + 7)
  G <- erdos.renyi.game(N, p, type=c("gnp"), directed = FALSE, loops = F) %>%
    as_adjacency_matrix(sparse = F)
  
  ## Specify the number of time periods ##
  periods <- if(periods_type == "squared"){
    N^2 + N + 1
  } else if(periods_type == "plus2"){
    N + 2
  }
  
  d_max = max(rowSums(G))
  
  set.seed(x)
  rho_1_values <- if(grepl("unif", rho_1)){
    range <- str_extract_all(pattern = "-?\\.?[0-9]", rho_1)
    runif(N, as.numeric(range[[1]][1]), as.numeric(range[[1]][2]))
  } else if(rho_1  == "seq(1/N, 1, 1/N)"){
    seq(1/N, 1, 1/N)
  } else if (is.numeric(rho_1)){
    rep(rho_1, N)
  } 
  
  
  set.seed(x + 7)
  rho_2_values <- if(rho_2 == "standard normal"){
    rnorm(N, 0, 1)
  } else if(rho_2 == "unif(-1, 1)"){
    runif(N, -1, 1)
  } else if(rho_2 == "norm(0.1, 1)"){
    rnorm(N, 0.1, 1)
  } else if(rho_2  == "seq(1/N, 1, 1/N)"){
    seq(1/N, 1, 1/N)
  } else if(rho_2 == "zero"){
    rep(0, N)
  } else if(rho_2 == "depends on rho_1"){
    (1-rho_1)/d_max - 0.01
  }else if (is.numeric(rho_2)){
    rho_2
  } 
  


  ## determine alpha, beta, and epsilon ##
  set.seed(3*x + 7)
  alpha = rnorm(N, mean = 0, sd = s_alpha)
  
  set.seed(4*x + 7)
  beta = rnorm(N, mean = 0, sd = s_beta)
  
  set.seed(5*x + 7)
  epsilon = rnorm(N, mean = 0, sd = s_epsilon)
  
  
  ## Generate Y ##
  
  set.seed(2*x + 7)
  y_initial_values <- if(y_initial == "standard normal"){
    rnorm(N, 0, 1)
  } else if(y_initial == "zero"){
    rep(0, N)
  }
  
  Y <- matrix(nrow = N, ncol = periods + 1)
  colnames(Y) <- 0:periods
  
  Y[, "0"] <- y_initial_values
  
  for(t1 in 1:periods){
    
    t = t1 - 1
    
    y_t <- Y[, paste(t)]
    
    y_t1 <- alpha + beta + rho_1_values * y_t + (rho_2_values * G %*% y_t )+ epsilon
    
    Y[, paste(t1)] <- y_t1
    
  }
  
  comparison_matrix = (diag(N) + G)
  comparison_matrix2 = (rho_1_values * diag(N) + rho_2_values * G)

  return(list(Y = Y, N = N, periods = periods, periods_type = periods_type, rankY = rankMatrix(Y)[1],
              s_epsilon = s_epsilon, 
              s_alpha = s_alpha,
              s_beta = s_beta, 
              p_num = p_num, 
              rho_1 = rho_1, 
              rho_2 = rho_2_values, 
              d_max = d_max,
              y_initial = y_initial, seed = seed, 
              comparison_matrix = comparison_matrix, 
              comparison_matrix2 = comparison_matrix2, 
              G=G))
  
}

1.2 g(\(\theta\)) function

# g of theta  
g_theta <- function(matrix, theta){
  ifelse(abs(matrix) > theta, 1, 0)
}

1.3 original disagreements function

# difference
disagreements <- function(original,estimate, N, theta){
  return(sum(abs(original - g_theta(estimate, theta)) == 1)/N^2)
}

1.4 sarah disagreements function

disagreements2 <- function(original, estimate, N){
  return(sum(abs(original - estimate))/N^2)
}

1.5 Ridge

### Ridge ###
ridge <- function(output, lambda, theta){
  
  Y<- output$Y
  
  
  N <- output$N
  periods <- output$periods
  
  F_Y <- c()
  Y_v <- c()
  
  for(i in 1:periods){
    t <- Y[, paste(i - 1)]
    
    f_t <- matrix(NA, ncol = N*N, nrow = N)
    
    for(j in 1:N){
      f_t[j, ]<- c(rep(rep(0, N), j-1), t(t), rep(rep(0, N), N-j))
    }
    
    F_Y <- rbind(F_Y, f_t)
  
    Y_v <- rbind(Y_v, t(t(Y[, paste(i)])))
  
  }
  
  
  A = try(solve(t(F_Y) %*% F_Y + lambda * diag(N^2)) %*% t(F_Y) %*% Y_v, silent = T)
  A <- t(matrix(A, N, N))
  
  error <- ifelse(sum(ifelse(class(A) == "try-error", 1, 0)) >0, "error", "no error")
  error_message <- ifelse(error  == "error", paste(attributes(A)$condition), "")

  dis <- disagreements(original = output$comparison_matrix, estimate = A, N = output$N, theta = theta)
  dis2 <- disagreements2(output$comparison_matrix2, A, output$N)
  diagnostics <- rbind(N= output$N, periods = output$periods_type, rankY = output$rankY, 
                       s_alpha = output$s_alpha,s_epsilon =  output$s_epsilon, 
                       s_beta = output$s_beta, p_num = output$p_num, rho_1 = output$rho_1, 
                       rho_2 = output$rho_2, y_initial = output$y_initial, error = error, 
                       error_message = error_message, seed = output$seed, lambda = lambda, theta = theta,
                       disagreements = round(dis, 3), 
                       disagreements2 = round(dis2, 3))

  return(list(A = A, diagnostics = diagnostics, comparison_matrix = output$comparison_matrix, comparison_matrix2= output$comparison_matrix2))
  
}

2 Output

We show the output for \(N = 15\). We set \(\lambda = 0.00001\) (because we previously found that \(\lambda = 0.00001\) vs \(\lambda = 0.1\) doesn’t make a difference).

We vary the following

variable	options
periods_type	squared, plus 2
\(s_\alpha\)	0, 0.2
\(s_\beta\)	0, 0.2
\(s_\epsilon\)	0, 0.2
\(\rho_1\)	0.2, 0.6
\(y_0\)	zero, standard normal
\(p_{numerator}\)	low (0.2), high (0.8)
\(\theta\)	0.1, 0.01, 0.0001, 0.0001

We base \(\rho_2\) off the following:

\(\rho_1 + d_{max}\rho_2 < 1 \implies \rho_2 < \frac{1 - \rho_1}{d_{max}}\), so we pick

\(\rho_2 = \frac{1 - \rho_1}{d_{max}} - 0.01\)

N = 15
periods_type =c("squared", "plus2")
s_alpha = c(0, 0.2)
s_beta = c(0, 0.2)
s_epsilon = c(0, 0.2)
rho_1 = c(0.2, 0.6)
y_initial = c("standard normal", "zero")
p_num = c("low", "high")
theta = c(0.1, 0.001)
lambda = 0.00001


tests <- expand.grid(N = N, 
                     periods_type = periods_type, 
                     s_alpha = s_alpha, 
                     s_beta = s_beta,
                     s_epsilon = s_epsilon, 
                     rho_1 = rho_1,
                     y_initial = y_initial, 
                     p_num = p_num,
                     lambda = lambda,
                     theta= theta,
            stringsAsFactors = F)

cl <- parallel::makeCluster(3)
doParallel::registerDoParallel(cl)
tbl <- foreach(i = 1:nrow(tests), .combine = rbind, .packages=c('tidyverse', 'igraph', 'Matrix')) %dopar% {
  
  output <- dgp(N = tests[i, "N"], 
                periods_type = tests[i, "periods_type"], 
                p_num = tests[i, "p_num"], 
                s_epsilon = tests[i, "s_epsilon"],  
                s_alpha =  tests[i, "s_alpha"], 
                s_beta = tests[i, "s_beta"], 
                rho_1 = tests[i, "rho_1"],  
                rho_2 = "depends on rho_1",
                y_initial = tests[i, "y_initial"], 
                seed = 1
  )
  
  out <- ridge(output = output, lambda = tests[i, "lambda"], theta = tests[i, "theta"])
  
  final <- t(out$diagnostics)
  
  return(final)
  
}
parallel::stopCluster(cl)

2.1 Trying Two Theta Values

2.1.1 Output

tbl <- as.data.frame(tbl) %>%
  mutate(across(c(N, rankY, s_alpha, s_epsilon, s_beta, rho_1, rho_2, lambda, theta, disagreements, disagreements2), as.numeric),
         across(c(rho_2), function(x) round(x, 3)))

options(DT.options = list(pageLength = 25))
tbl %>% 
  datatable() 

2.1.2 Plot Original Disagreement vs New Disagreement

plot(tbl$disagreements, tbl$disagreements2)

2.1.3 Comparing thetaas

Is 0.1 or 0.001 more successful at minimizing disagreements?

There are 60 cases where 0.001 has lower disagreement values and 60 cases where 0.1 has lower disagreement rates.

Note that this follows high (0.8) and low (0.2) \(p_{numerator}\). In all cases where \(p_{numerator}\) is low, 0.1 does better (or equal). In all cases where \(p_{numerator}\) is high, 0.001 has lower disagreement values ( or equal).

tbl_summ <- tbl %>%
  select(N, periods, rankY, s_alpha, s_epsilon, s_beta, p_num, rho_1, rho_2, y_initial, lambda, theta, disagreements) %>%
  pivot_wider(names_from = theta, values_from = disagreements, names_prefix = "dis1_") %>%
  mutate(best_theta = case_when(dis1_0.1 > dis1_0.001 ~ "0.001", 
                                dis1_0.1 < dis1_0.001 ~ "0.1", 
                                TRUE ~ "equal"), 
         min = as.numeric(case_when(dis1_0.1 > dis1_0.001 ~ dis1_0.001, 
                         dis1_0.1 < dis1_0.001 ~ dis1_0.1,
                         TRUE ~ dis1_0.001))) 

tbl_summ %>%
  group_by(best_theta, p_num) %>%
  count()

## # A tibble: 4 × 3
## # Groups:   best_theta, p_num [4]
##   best_theta p_num     n
##   <chr>      <chr> <int>
## 1 0.001      high     60
## 2 0.1        low      60
## 3 equal      high      4
## 4 equal      low       4

2.1.4 Distribution of Minimum Disagreement Values

Let’s say that we choose \(\theta\) such that disagreements are minimized. What does the distribution look like?

plot(density(tbl_summ$min), main = "Distribution of Minimum Disagreement Values")

2.2 Four Theta Values

Now, let’s add \(\theta \in \{0.0001, 1\}\).

N = 15
periods_type =c("squared", "plus2")
s_alpha = c(0, 0.2)
s_beta = c(0, 0.2)
s_epsilon = c(0, 0.2)
rho_1 = c(0.2, 0.6)
y_initial = c("standard normal", "zero")
p_num = c("low", "high")
theta = c(1, 0.0001)
lambda = 0.00001


tests2 <- expand.grid(N = N, 
                     periods_type = periods_type, 
                     s_alpha = s_alpha, 
                     s_beta = s_beta,
                     s_epsilon = s_epsilon, 
                     rho_1 = rho_1,
                     y_initial = y_initial, 
                     p_num = p_num,
                     lambda = lambda,
                     theta= theta,
            stringsAsFactors = F)

cl <- parallel::makeCluster(3)
doParallel::registerDoParallel(cl)
tbl2 <- foreach(i = 1:nrow(tests2), .combine = rbind, .packages=c('tidyverse', 'igraph', 'Matrix')) %dopar% {
  
  output <- dgp(N = tests2[i, "N"], 
                periods_type = tests2[i, "periods_type"], 
                p_num = tests2[i, "p_num"], 
                s_epsilon = tests2[i, "s_epsilon"],  
                s_alpha =  tests2[i, "s_alpha"], 
                s_beta = tests2[i, "s_beta"], 
                rho_1 = tests2[i, "rho_1"],  
                rho_2 = "depends on rho_1",
                y_initial = tests2[i, "y_initial"], 
                seed = 1
  )
  
  out <- ridge(output = output, lambda = tests2[i, "lambda"], theta = tests2[i, "theta"])
  
  final <- t(out$diagnostics)
  
  return(final)
  
}
parallel::stopCluster(cl)

2.2.1 Output Table

tbl2 <- as.data.frame(tbl2)%>%
  mutate(across(c(N, rankY, s_alpha, s_epsilon, s_beta, rho_1, rho_2, lambda, theta, disagreements, disagreements2), as.numeric),
         across(c(rho_2), function(x) round(x, 3)))

final_df <- rbind(tbl, tbl2)
final_df %>% 
  datatable() 

2.2.2 Plotting original vs new disagreement

Is there a relationship between original and new disagreements measures?

plot(final_df$disagreements, final_df$disagreements2)

2.2.3 Comparing thetas

Which values of \(\theta\) minimize disagreements?

final_summ <- final_df %>%
  select(N, periods, rankY, s_alpha, s_epsilon, s_beta, p_num, rho_1, rho_2, y_initial, lambda, theta, disagreements) %>%
  arrange(theta) %>%
  pivot_wider(names_from = theta, values_from = disagreements, names_prefix = "dis1_")


best <- apply(X = final_summ[, c("dis1_0.0001",  "dis1_0.001", "dis1_0.1","dis1_1")], MARGIN = 1, FUN =function(x){ which(x == min(x))}) %>%
  lapply(., names)

final_summ$`min_0.0001` <- sapply(best, function(x)is.element("dis1_0.0001", x), simplify = T) %>% as.numeric
final_summ$`min_0.001` <- sapply(best, function(x)is.element("dis1_0.001", x), simplify = T) %>% as.numeric
final_summ$`min_0.1` <- sapply(best, function(x)is.element("dis1_0.1", x), simplify = T) %>% as.numeric
final_summ$`min_1` <- sapply(best, function(x)is.element("dis1_1", x), simplify = T) %>% as.numeric
final_summ$best_theta <- lapply(best, function(x){ifelse(length(x) >1, paste(x, collapse = ", "), x)}) %>% as.character

final_summ <- final_summ %>%
  mutate(min = case_when(min_1 == 1 ~ dis1_1, 
                         min_0.1 == 1 ~ dis1_0.1, 
                        min_0.001 == 1 ~ dis1_0.001, 
                        min_0.0001 == 1 ~ dis1_0.0001, 
                        TRUE ~ min_1))


final_summ %>%
  group_by(min_1, min_0.1, min_0.001, min_0.0001, p_num) %>%
  count() %>%
  arrange(p_num)

## # A tibble: 8 × 6
## # Groups:   min_1, min_0.1, min_0.001, min_0.0001, p_num [8]
##   min_1 min_0.1 min_0.001 min_0.0001 p_num     n
##   <dbl>   <dbl>     <dbl>      <dbl> <chr> <int>
## 1     0       0         0          1 high     25
## 2     0       0         1          0 high      6
## 3     0       0         1          1 high     29
## 4     1       1         1          1 high      4
## 5     0       1         0          0 low       7
## 6     1       0         0          0 low      50
## 7     1       1         0          0 low       3
## 8     1       1         1          1 low       4

2.2.4 Distribution of Minimum Disagreement Values

Let’s say that we choose \(\theta\) such that disagreements are minimized. What does the distribution look like?

plot(density(final_summ$min), main = "Distribution of Minimum Disagreement Values")