tylermorganwall / skpr

Generates and evaluates D, I, A, Alias, E, T, G, and custom optimal designs. Supports generation and evaluation of mixture and split/split-split/N-split plot designs. Includes parametric and Monte Carlo power evaluation functions. Provides a framework to evaluate power using functions provided in other packages or written by the user.
https://tylermorganwall.github.io/skpr/
GNU General Public License v3.0
115 stars 15 forks source link

Variance/covariance matrix wrong by wrong blocking structure reconstruction if some only one block plots in a given blocking level #80

Closed caebergs closed 8 months ago

caebergs commented 1 year ago

Hello,

In R/gen_design.R, calls to :

 blocklist = strsplit(initialrownames, ".", fixed = TRUE)

 existingblockstructure = do.call(rbind, blocklist)

 blockgroups = apply(existingblockstructure, 2, blockingstructure)

when two rows from separate harder to change level have "1", result in being associated to the same block. For example, 6.1.1.1 and 5.1.1.1 will be considered from the same block with the current blocking structure.

Replacing (separating by index within the same blocking level, ignoring the harder to change levels)

 blocklist = strsplit(initialrownames, ".", fixed = TRUE)

by something similar to (reforming the block identifier, including harder levels :

 headsgroup=function(initrownames,n) {

   blocklist=unlist(strsplit(initrownames,"\\."),recursive=FALSE);

   return(sapply(seq(1,length(blocklist)),function(n){paste0(blocklist[1:n],collapse=".")}))

 }

 existingblockstructure=do.call(rbind,lapply(initialrownames,headsgroup))

This code is obviously to be cleaned for readability.

Thank you in advance for implementing this solution or proposing an alternative solution.

Best regards,

Thierry

tylermorganwall commented 8 months ago

This has been fixed in the latest (v1.7.0) update with an overhaul of how skpr parses rownames. Rather than the previous behavior of marking the position within the block, such as:

       a  b  c
1.1.1  1  1  1
1.1.2  1  1 -1
1.2.1  1 -1  1
1.2.2  1 -1 -1
2.1.1 -1 -1  1
2.1.2 -1 -1 -1
2.2.1 -1  1 -1
2.2.2 -1  1  1
3.1.1 -1 -1  1
3.1.2 -1 -1 -1
3.2.1 -1  1  1
3.2.2 -1  1 -1
4.1.1  1  1 -1
4.1.2  1  1  1
4.2.1  1 -1 -1
4.2.2  1 -1  1
5.1.1  1  1  1
5.1.2  1  1 -1
5.2.1  1 -1 -1
5.2.2  1 -1  1
6.1.1 -1  1  1
6.1.2 -1  1 -1
6.2.1 -1 -1  1
6.2.2 -1 -1 -1

skpr now just numbers each block according to the order within the stratum in which it lies, independent of other layers.

         a  b  c
1.1.1   -1  1  1
1.1.2   -1  1 -1
1.2.3   -1 -1  1
1.2.4   -1 -1 -1
2.3.5   -1 -1 -1
2.3.6   -1 -1  1
2.4.7   -1  1  1
2.4.8   -1  1 -1
3.5.9    1 -1 -1
3.5.10   1 -1  1
3.6.11   1  1 -1
3.6.12   1  1  1
4.7.13   1  1  1
4.7.14   1  1 -1
4.8.15   1 -1  1
4.8.16   1 -1 -1
5.9.17  -1 -1 -1
5.9.18  -1 -1  1
5.10.19 -1  1 -1
5.10.20 -1  1  1
6.11.21  1 -1  1
6.11.22  1 -1 -1
6.12.23  1  1 -1
6.12.24  1  1  1