From b4aba3cea64aaa51cb4ce6657fad647720262015 Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Mon, 7 Jun 2021 18:05:37 +0200
Subject: [PATCH] Fix accuracy.R + add postTreatment draft

---
 reports/accuracy.R      |  4 ++--
 reports/postTreatment.R | 39 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 41 insertions(+), 2 deletions(-)
 create mode 100644 reports/postTreatment.R

diff --git a/reports/accuracy.R b/reports/accuracy.R
index 9a9d21d..63cd5aa 100644
--- a/reports/accuracy.R
+++ b/reports/accuracy.R
@@ -27,9 +27,9 @@ optimBeta <- function(N, n, p, beta, b, link, ncores)
         dat <- as.data.frame( cbind(fargs$Y,fargs$X) )
         res2 <- NULL
         tryCatch({
-          fm <- flexmix( cbind(V1, 1-V1) ~ .-V1, data=dat, k=K,
+          fm <- flexmix( cbind(V1, 1-V1) ~ ., data=dat, k=K,
             model = FLXMRglm(family = binomial(link = link)) )
-          pf <- mean(fm@posterior[["scaled"]][,1])
+          pf <- colMeans(fm@posterior[["scaled"]])
           out <- refit(fm)
           beta_b <- sapply( seq_len(K), function(i) {
             as.double( out@components[[1]][[i]][,1] )
diff --git a/reports/postTreatment.R b/reports/postTreatment.R
new file mode 100644
index 0000000..b9ad42a
--- /dev/null
+++ b/reports/postTreatment.R
@@ -0,0 +1,39 @@
+K <- list("2"=2, "5"=2, "10"=3)
+p <- list("2"=c(.5,.5), "5"=c(.5,.5), "10"=c(.3, .3, .4))
+b <- list("2"=c(-.2, .5), "5"=c(-.2, .5), "10"=c(-.2, 0, .5))
+beta <- list(
+  "2"=matrix( c(1,-2, 3,1), ncol=2 ),
+  "5"=matrix( c(1,2,-1,0,3, 2,-3,0,1,0), ncol=2 ),
+  "10"=matrix( c(1,2,-1,0,3,4,-1,-3,0,2, 2,-3,0,1,0,-1,-4,3,2,0, -1,1,3,-1,0,0,2,0,1,-2), ncol=3 ))
+
+for (n in c("5000", "10000", "100000")) {
+  for (d in c("2", "5", "10")) {
+    load(paste0("res_",n,"_",d,"_logit.RData"))
+
+    # p
+    for (m in 1:2) {
+      err <- c()
+      for (i in 1:K[[d]])
+        err <- c(err, mean(abs(morpheus:::.extractParam(mr, 1, i)[[1]] - p[[d]][i])))
+      print(paste0("p", m, " ", mean(err)))
+    }
+
+    # b
+    for (m in 1:2) {
+      err <- c()
+      for (i in 1:K[[d]])
+        err <- c(err, mean(abs(morpheus:::.extractParam(mr, 2+d, i)[[1]] - b[[d]][i])))
+      print(paste0("b", m, " ", mean(err)))
+    }
+
+    # beta
+    for (m in 1:2) {
+      for (i in 1:K[[d]]) {
+        err <- c()
+        for (j in 1:d)
+          err <- c(err, mean(abs(morpheus:::.extractParam(mr, 1+j, i)[[1]] - beta[[d]][1+j,i])))
+        print(paste0("beta", m, "_", i, " ", mean(err)))
+      }
+    }
+  }
+}
-- 
2.44.0