Info objective - #91
Conversation
|
Thank you for additional objective function. |
64f1b03 to
c89070a
Compare
|
Thanks, i updated the yaml so the new public functions are in the index. The low level helpers are already |
| objective_function_info <- function( | ||
| layout_df, | ||
| swap, | ||
| spatial_cols, | ||
| criterion = c("A", "D"), | ||
| L_matrix = NULL, | ||
| block_column = "block", | ||
| ... | ||
| ) { | ||
| criterion <- match.arg(criterion) | ||
|
|
||
| ci <- .compute_info(layout_df, swap, L_matrix, block_column) | ||
| info <- ci$info | ||
| v <- ci$v | ||
|
|
||
| ## get positive eigenvals from info matrix | ||
| eig <- eigen(info, symmetric = TRUE, only.values = TRUE)$values | ||
| max_eig <- max(eig) | ||
| if (max_eig <= 0) { | ||
| score <- 1e10 | ||
| pos_eig <- numeric(0) | ||
| } else { | ||
| pos_eig <- eig[eig > max_eig * 1e-10] | ||
| ## if design is disconnected, penalise greatly | ||
| ## otherwise A or D | ||
| if (length(pos_eig) < v - 1) { | ||
| score <- 1e10 | ||
| } else { | ||
| score <- switch(criterion, | ||
| A = sum(1 / pos_eig), | ||
| D = -sum(log(pos_eig)) | ||
| ) | ||
| if (!is.finite(score)) score <- 1e10 | ||
| } | ||
| } | ||
|
|
||
| return(list( | ||
| score = score, | ||
| info_matrix = info, | ||
| eigenvalues = sort(pos_eig, decreasing = TRUE), | ||
| criterion = criterion | ||
| )) | ||
| } |
There was a problem hiding this comment.
We should be able to take advantage of current_score_obj and swapped_items passed to objective function calls from the 2nd iteration onwards, so only partial calculation or object assigning might be required each iteration (?).
Please see example usage in objective_function_piepho.
| #' @return \eqn{(n \times n)} numeric matrix. | ||
| #' | ||
| #' @export | ||
| compute_L_projection <- function( |
There was a problem hiding this comment.
Order matters here and speed calls order under the hood.
It might be better that the users do not have to call this themselves by calling this function in the 1st objective function call and pass around via current_score_obj.
|
|
||
| ## L = I - X2 (X2t X2)^{-1} X2t | ||
| X2tX2_inv <- diag(1 / colSums(X2), nrow = b) | ||
| diag(n) - X2 %*% X2tX2_inv %*% t(X2) |
There was a problem hiding this comment.
Please add an explicit return statement.
There are other places as well.
| if (nrow(Sigma) != n) { | ||
| stop( | ||
| "Sigma matrix must have dimension equal ", | ||
| "to dimension of layout_df dataframe") | ||
| } | ||
| if (ncol(Sigma) != n) { | ||
| stop( | ||
| "Sigma matrix must have dimension equal ", | ||
| "to dimension of layout_df dataframe" | ||
| ) | ||
| } |
There was a problem hiding this comment.
These can be collapsed into 1 if.
| #' Calculate concurrence matrix | ||
| #' | ||
| #' @param layout_df Data frame representing the spatial layout of the | ||
| #' experiment. | ||
| #' @param treatment_column Column name containing the design's treatments in | ||
| #' \code{layout_df}. | ||
| #' @param block_column Column name containing the designs block in | ||
| #' \code{layout_df}. | ||
| #' | ||
| #' @return Symmetric \eqn{N N^\intercal}{N Nᵀ} matrix. | ||
| #' | ||
| #' @keywords internal | ||
| calc_concurrence_matrix <- function( |
There was a problem hiding this comment.
I think this should be exported and live in metrics.R.
Maybe add brief description and use calculate_ to be consistent.
| #' Calculate canonical efficiency factors | ||
| #' | ||
| #' @param layout_df Data frame representing the spatial layout of the | ||
| #' experiment. | ||
| #' @param treatment_column Column name containing the design's treatments in | ||
| #' \code{layout_df}. | ||
| #' @param L_matrix Precomputed projection matrix from | ||
| #' \code{compute_L_projection}. | ||
| #' @param block_column Column name containing the designs block in | ||
| #' \code{layout_df}. | ||
| #' | ||
| #' @return A list with \code{efficiency_factors}, \code{E}, \code{replication}. | ||
| #' | ||
| #' @keywords internal | ||
| calculate_efficiency_factors <- function( |
There was a problem hiding this comment.
Should be exported and live in metrics.R as well with brief description.
There was a problem hiding this comment.
The name should be objective_{something}.R, like objective_fisher.
You got the idea.
Also, I'll check the math later (wish me luck).
There was a problem hiding this comment.
I'd move metric functions like efficiency and matrices out of this file.
Can even do 1 file per function like a bunch of existing test-calculate_*s.
2 participants
Fisher information based objective function, with A and D optimality criterion, and optimises spatially using a user-supplied covariance structure. It is calculated from treatment contrasts, adjusted for nuisance fixed effects like block and intercept.
And L removes nuisance effects under a spatial correlation structure from user-supplied$n \times n$ covariance $\Sigma$ and nuisance effects $X_2$
With A-optimality$\mathrm{tr}(I^-)$ or contrast-space D-optimality (pseudo determinant since I has rank $\nu - 1$ ) $-\log \left| I \right|$ .
So basically the user specifies the correlation structure they want to use at analysis, such as$AR1 \times AR1$ , and the design is optimised to be efficient under that spatial structure. So it should be directly targeting the Cramér-Rao lower bound rather than the adjacency proxy currently used in
speed.