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

Can generate but not evaluate design - undefined columns selected #54

Closed comatrion closed 4 years ago

comatrion commented 4 years ago

Hi

The script below is auto-generated by the skprGUI. I believe it generates a valid split-plot design (X3 is HTC) but when I attempt to evaluate the design I get an error.

The error

Error in `[.data.frame`(run_matrix_processed, , col_name) : 
  undefined columns selected
Calls: eval_design ... calculate_degrees_of_freedom -> unique -> [ -> [.data.frame
Execution halted

The script

# This is the R code used to generate these results in skpr.
# Copy this into an R script and rerun to reproduce these results.

library(skpr)

# Generating candidate set:
candidateset = expand.grid(X1 = seq(-1,1, length.out = 3), 
                           X2 = seq(-1,1, length.out = 3), 
                           X3 = seq(-1,1, length.out = 3))

# Generating design for hard-to-change factors: 
design_htc = gen_design(candidateset = candidateset, 
                        model = ~X3, 
                        trials = 4)

# Generating design:
design = gen_design(candidateset = candidateset, 
                    model = ~X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3 + I(X1^2) + I(X2^2), 
                    trials = 48, 
                    splitplotdesign = design_htc, 
                    splitplotsizes = 12, 
                    splitcolumns = TRUE)

# Evaluating Design:
eval_design(design = design, 
            model = ~X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3 + I(X1^2) + I(X2^2), 
            alpha = 0.05, 
            blocking = TRUE)

## How to analyze this experiment when the data have been collected:
## (to run, remove one # from this section) 
## First, assign the results to a column in the data frame. Each 
## entry in the vector corresponds to the result from that run in the design. 

#design$Y = results 

## Now analyze the blocked linear model with lmer (from the lme4 package) and lmerTest:

#library(lmerTest)
#lme4::lmer(formula = Y ~ (1|Block1) + X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3 +
 I(X1^2) + I(X2^2), data = design))

Env: R 3.6.1, skpr 0.62.0

comatrion commented 4 years ago

I think the error might possibly be occurring in the call to unique here https://github.com/tylermorganwall/skpr/blob/72350dd2d9e82eb89252c932ad0a706ac4ec4f06/R/calculate_degrees_of_freedom.R#L41-L45 and is because at this point splitterms[j] contains I(X1^2) which is not in the expected format of an interaction term e.g. X1:X2.

comatrion commented 4 years ago

If there was an error it seems to have been resolved - the script runs without error in 0.64.1.

Loading required package: shiny
     parameter            type     power
1  (Intercept)    effect.power 0.1539629
2           X1    effect.power 0.9999705
3           X2    effect.power 0.9999794
4           X3    effect.power 0.2063510
5      I(X1^2)    effect.power 0.7677688
6      I(X2^2)    effect.power 0.7370548
7        X1:X2    effect.power 0.9998038
8        X1:X3    effect.power 0.9999705
9        X2:X3    effect.power 0.9999794
10 (Intercept) parameter.power 0.1539629
11          X1 parameter.power 0.9999705
12          X2 parameter.power 0.9999794
13          X3 parameter.power 0.2063510
14     I(X1^2) parameter.power 0.7677688
15     I(X2^2) parameter.power 0.7370548
16       X1:X2 parameter.power 0.9998038
17       X1:X3 parameter.power 0.9999705
18       X2:X3 parameter.power 0.9999794
============Evaluation Info=============
• Alpha = 0.05 • Trials = 48 • Blocked = TRUE
• Evaluating Model = ~X1 + X2 + X3 + X1:X2 + X1:X3 + X2:X3 + I(X1^2) + I(X2^2)
• Anticipated Coefficients = c(1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000)
• Number of Blocks = Level 1: 4
• Variance Ratios  = Level 1: 1