Dual use of probz accessor

[dmckean]

probz's accessor includes a warning that it is modifying the value.

## TODO Dangerous to have accessor do something more than return the value of it
## slot.  Further
## probz(object) <- probz(object)
## will not behave as expected
#' @rdname probz-method
#' @aliases probz,MixtureModel-method
setMethod("probz", "MixtureModel", function(object) {
  ## because first iteration not saved
  object@probz/(iter(object)-1)
})

I think that probz is written this way because it carries unscaled probabilities that are added to during model fitting. The CPP functions update_probz_batch and compute_probz support this. R functions combine_batch, combine_multibatch_pooled and .computeProbz also have code that indicate that the probabilities must be rescaled to add up to 1.

Coercion methods do not seem to account for this rescaling. Should these, and other functions that act as constructors, copy the probz matrix directly using a slot accessor?

As currently defined, the matrix rows in a CopyNumberModel's probz slot is scaled to add up to 1. As a result, the rowSums of probz(cn.model) or probCopyNumber(cn.model) add up to 1/(cn.model@mcmc.parms@iter - 1)

[rscharpf]

Reviewing the unit tests for the copy number models (test_copynumber_models.R), there is a function called probCopyNumber that provides probabilities that are scaled to sum to 1. I think this is what should be exposed to the user rather than probz. On February 5, 2018 at 5:01:30 PM, Dave McKean (notifications@github.com<mailto:notifications@github.com>) wrote: probz's accessor includes a warning<https://github.com/scristia/CNPBayes/blob/991b84bbf4093b7d49676d3cba4c72cfea078db0/R/methods-MixtureModel.R#L216-L225> that it is modifying the value. ## TODO Dangerous to have accessor do something more than return the value of it ## slot. Further ## probz(object) <- probz(object) ## will not behave as expected #' @Rdname probz-method #' @Aliases probz,MixtureModel-method setMethod("probz", "MixtureModel", function(object) { ## because first iteration not saved object@probz/(iter(object)-1) }) I think that probz is written this way because it carries unscaled probabilities that are added to during model fitting. The CPP functions update_probz_batch and compute_probz support this. R functions combine_batch, combine_multibatch_pooled and .computeProbz also have code that indicate that the probabilities must be rescaled to add up to 1. Coercion methods<https://github.com/scristia/CNPBayes/blob/0ff9ecd54c69ced186f9993adde812fb56bc2e81/R/coerce-methods.R#L35> do not seem to account for this rescaling. Should these, another functions that act as constructors copy the probz matrix directly using a slot accessor? As currently defined, the matrix rows in a CopyNumberModel's probz slot is scaled to add up to 1. As a result, the rowSums of probz(cn.model) or probCopyNumber(cn.model) add up to 1/(cn.model@mcmc.parms@iter - 1) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#10>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AAyQw8APb6teN7t4rcOoyqICvGmPsF5Dks5tR3ougaJpZM4R6KsP>. WARNING: E-mail sent over the Internet is not secure. Information sent by e-mail may not remain confidential. DISCLAIMER: This e-mail is intended only for the individual to whom it is addressed. It may be used only in accordance with applicable laws. If you received this e-mail by mistake, notify the sender and destroy the e-mail.

[dmckean]

I am not observing this problem in models that are created from a new constructor.

I think that we have a problem where coercing a model will rescale the probz matrix to 1, leading to the probz accessor in the subsequent model to incorrectly scale the values. This is true regardless of if the output is coming from probz() or probCopyNumber().

probCopyNumber() calls .prob_copynumber(), I don't see it doing any scaling to one other than in case 2:

Case 1. If the model is does not have manyToOneMapping, .prob_copynumber() returns the probz matrix via the probz() accessor.
Case 2. If the model mapping reduces to a single state, it returns a Nx1 matrix of ones.
Case 3. Otherwise, return a matrix where probz columns (as returned by probz()) have been added together for states collapsed by the copy number mapping.

I'm observing this problem in probz() matrices for any pooled model. This problem persists for any pooled copy number model, except when Case 2 applies. This is happening because these models are being created via coercion methods that rely on the the probz() accessor.

All of the MixtureModel-derived coercion methods copy over the probz matrix like so

    probz=probz(from),

I think the solution would be to either ditch the accessor function within coercion methods (perhaps probz=from@probz, or to rescale them, as done in combineModels() in diagnostics.R probz=probz(from) * (iter(from) - 1).

These are the defined coercion methods.

setAs("BatchModel", "MultiBatchModel", function(from, to))
# used to create MultiBatchModelExample, which you can see is improperly scaled by running 

summary(rowSums(probz(MultiBatchModelExample)))
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#0.003356 0.003356 0.003356 0.003356 0.003356 0.003356 

setAs("MarginalModel", "SingleBatchModel", function(from, to))
# same applies here, try
summary(rowSums(probz(SingleBatchModelExample)))
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#   0.25    0.25    0.25    0.25    0.25    0.25 

setAs("SingleBatchPooled", "SingleBatchCopyNumber", function(from, to)) - used in SingleBatchCopyNumber() in copynumber-models.R.

setAs("MultiBatchPooled", "MultiBatchCopyNumberPooled", function(from, to) - used in MultiBatchCopyNumberPooled() in copynumber-models.R

setAs("SingleBatchModel", "SingleBatchPooled", function(from, to) - used in SingleBatchPooled(), .SingleBatchPooled() in methods-SingleBatchPooled.R

setAs("MultiBatchModel", "MultiBatchPooled", function(from, to) - used in MultiBatchPooled() in methods-MultiBatchPooled.R (2 uses)

[rscharpf]

Thanks for these details. I think the best way to proceed would be to add a unit test for each of these assertions. They should all fail. When they stop failing, we will have some assurance that they've been addressed and we can avoid repeating past mistakes in the future