Nanostring-Biostats / GeoDiff

GeoDiff, an R package for count generating models for analyzing Geomx RNA data. Note that this version of the package is still under development, undergoing submission process to Bioconductor 3.14 release and still needs to complete NanoString internal verification process.
MIT License
7 stars 6 forks source link

fitNBthDE error with prior_type="equal" #10

Closed maddygriz closed 2 years ago

maddygriz commented 2 years ago

This error occurs when prior_type = "equal". It doesn't happen every time but I don't see a pattern for when it does.

Error in t(X) %*% preci1con : non-conformable arguments

The Traceback doesn't seem helpful but I'm adding it anyway. Traceback:

  1. .local(object, ...)
  2. fitNBthDE(form = form, annot = annot[ROIs_high, ], object = countmat[, ROIs_high], probenum = probenum, features_high = features_high, features_all = features_all, sizefact_start = sizefact_start, sizefact_BG = sizefact_BG, threshold_mean = threshold_mean, ...
  3. fitNBthDE(form = form, annot = annot[ROIs_high, ], object = countmat[, ROIs_high], probenum = probenum, features_high = features_high, features_all = features_all, sizefact_start = sizefact_start, sizefact_BG = sizefact_BG, threshold_mean = threshold_mean, ...
  4. .local(object, ...)
  5. fitNBthDE(form = ~region, split = FALSE, object = kidney, preci2 = 10000, prior_type = "equal", preci1con = 1/25, sizescalebythreshold = TRUE)
  6. fitNBthDE(form = ~region, split = FALSE, object = kidney, preci2 = 10000, prior_type = "equal", preci1con = 1/25, sizescalebythreshold = TRUE)
maddygriz commented 2 years ago

I am getting the same error with fitPoisthNorm()

zhiiiyang commented 2 years ago

fitNBthDE passed

> data(demoData)
> demoData <- demoData[, c(1:5, 33:37)]
> demoData <- fitPoisBG(demoData, size_scale = "sum")
Iteration = 1, squared error = 1.199470e+05
Iteration = 2, squared error = 0.000000e+00
Model converged.
> demoData <- aggreprobe(demoData, use = "cor")
> demoData <- BGScoreTest(demoData)
> demoData$slidename <- substr(demoData[["slide name"]], 12, 17)
> thmean <- 1 * mean(fData(demoData)$featfact, na.rm = TRUE)
> demo_pos <- demoData[which(!fData(demoData)$CodeClass == "Negative"), ]
> demo_neg <- demoData[which(fData(demoData)$CodeClass == "Negative"), ]
> sc1_scores <- fData(demo_pos)[, "scores"]
> names(sc1_scores) <- fData(demo_pos)[, "TargetName"]
> features_high <- ((sc1_scores > quantile(sc1_scores, probs = 0.4)) &
+                     (sc1_scores < quantile(sc1_scores, probs = 0.95))) |>
+   which() |>
+   names()
> set.seed(123)
> demoData <- fitNBth(demoData,
+                     features_high = features_high,
+                     sizefact_BG = demo_neg$sizefact,
+                     threshold_start = thmean,
+                     iterations = 5,
+                     start_para = c(200, 1),
+                     lower_sizefact = 0,
+                     lower_threshold = 100,
+                     tol = 1e-8)
Iteration = 1, squared error = 0.00211581666061465
Iteration = 2, squared error = 0.000635277560607171
Iteration = 3, squared error = 0.000136530888681976
Iteration = 4, squared error = 3.81046696630021e-05
Iteration = 5, squared error = 1.15557488250469e-05
Model converged.
> ROIs_high <- sampleNames(demoData)[which(demoData$sizefact_fitNBth * thmean > 2)]
> features_all <- rownames(demo_pos)
> 
> pData(demoData)$group <- c(rep(1, 5), rep(2, 5))
> 
> NBthDEmod1 <- fitNBthDE(
+   form = ~group,
+   split = FALSE,
+   object = demoData,
+   ROIs_high = ROIs_high,
+   features_high = features_high,
+   features_all = features_all,
+   sizefact_start = demoData[, ROIs_high][["sizefact_fitNBth"]],
+   sizefact_BG = demoData[, ROIs_high][["sizefact"]],
+   preci2 = 10000,
+   prior_type = "equal",
+   covrob = FALSE,
+   preci1con = 1/25,
+   sizescalebythreshold = TRUE
+ )
Iteration = 1, squared error = 4.081296e-05
Iteration = 2, squared error = 1.933278e-05
zhiiiyang commented 2 years ago

fitPoisthNorm passed

> library(Biobase)
> library(dplyr)
> data(demoData)
> demoData <- fitPoisBG(demoData, size_scale = "sum")
Iteration = 1, squared error = 9.471314e+06
Iteration = 2, squared error = 0.000000e+00
Model converged.
> demoData <- aggreprobe(demoData, use = "cor")
> demoData <- BGScoreTest(demoData)
> thmean <- 1 * mean(fData(demoData)$featfact, na.rm = TRUE)
> demo_pos <- demoData[which(!fData(demoData)$CodeClass == "Negative"), ]
> demo_neg <- demoData[which(fData(demoData)$CodeClass == "Negative"), ]
> sc1_scores <- fData(demo_pos)[, "scores"]
> names(sc1_scores) <- fData(demo_pos)[, "TargetName"]
> features_high <- ((sc1_scores > quantile(sc1_scores, probs = 0.4)) &
+                     (sc1_scores < quantile(sc1_scores, probs = 0.95))) |>
+   which() |>
+   names()
> set.seed(123)
> features_high <- sample(features_high, 100)
> demoData <- fitNBth(demoData,
+                     features_high = features_high,
+                     sizefact_BG = demo_neg$sizefact,
+                     threshold_start = thmean,
+                     iterations = 5,
+                     start_para = c(200, 1),
+                     lower_sizefact = 0,
+                     lower_threshold = 100,
+                     tol = 1e-8)
Iteration = 1, squared error = 0.00127298898031945
Iteration = 2, squared error = 3.85529407814789e-05
Iteration = 3, squared error = 2.75554613662818e-06
Iteration = 4, squared error = 3.04023597732423e-07
Iteration = 5, squared error = 6.68035435765034e-08
Model converged.
> ROIs_high <- sampleNames(demoData)[which((quantile(fData(demoData)[["para"]][, 1],
+                                                    probs = 0.90, na.rm = TRUE) -
+                                             notes(demoData)[["threshold"]]) * demoData$sizefact_fitNBth > 2)]
> features_all <- rownames(demo_pos)
> thmean <- mean(fData(demo_neg)[["featfact"]])
> demoData <- fitPoisthNorm(
+   object = demoData,
+   split = FALSE,
+   ROIs_high = ROIs_high,
+   features_high = features_high,
+   features_all = features_all,
+   sizefact_start = demoData[, ROIs_high][["sizefact_fitNBth"]],
+   sizefact_BG = demoData[, ROIs_high][["sizefact"]],
+   threshold_mean = thmean,
+   preci2 = 10000,
+   prior_type = "equal",
+   covrob = FALSE,
+   preci1con = 1 / 25
+ )
probe finished
Iteration = 1, squared error = 5.48724768110968e-05
probe finished
Iteration = 2, squared error = 5.77221071476496e-05
Model converged.