demonstration-of-parts-of-the-package.Rmd
library(mostlytidyMMM)These functions define adstock and saturation functions ($ f(x) = A(1-e^{-kx} $) for use in recipes (ie tidymodels pre-model transformation pipelines).
#not run!
suppressMessages(suppressWarnings(library(tidyverse)))
suppressMessages(suppressWarnings(library(recipes)))
suppressMessages(suppressWarnings(library(tune)))
suppressMessages(suppressWarnings(library(dials)))
suppressMessages(suppressWarnings(library(mostlytidyMMM)))
#create two datasets
mktdata<-rbind(tibble(prod='brand',store='store1',
sales=c(100,100,100,100,100),tv=c(10,100,0,0,100),
search=c(0,10,20,50,50)) ,
tibble(prod='brand',store='store2',
sales=c(10,10,10,10,10),tv=c(0,0,0,0,0),
search=c(0,2,2,0,0) ) ) |>
group_by(prod,store)
mktdata2<-tibble(prod='brand',store='all',
sales=100,tv=1000,search=1000) |> group_by(prod,store)
#creating a recipe with a saturaiton step in it
rec_obj <-
recipe(sales ~ ., data = mktdata) |>
step_saturation(c(tv,search),asymptote=200,saturation_speed=.003 ) |>
prep(training = mktdata)
#showing what a tunable recipe looks liek:
rec_obj2<- recipe(sales ~ ., data = mktdata) |>
step_saturation(c(tv,search),asymptote=tune(),saturation_speed=tune() )
extract_parameter_set_dials(rec_obj2)
#> Collection of 2 parameters for tuning
#>
#> identifier type object
#> asymptote asymptote nparam[+]
#> saturation_speed saturation_speed nparam[+]
#custom print shows values of saturation hyperparameters:
print(rec_obj)
#> Recipe
#>
#> Inputs:
#>
#> role #variables
#> outcome 1
#> predictor 4
#>
#> Training data contained 10 data points and no missing data.
#>
#> Operations:
#>
#> Saturation (asymptote= 200 saturation_speed= 0.003 Transformation ontv, search [trained]
#or shows that they are set to be tuned:
print(rec_obj2)
#> Recipe
#>
#> Inputs:
#>
#> role #variables
#> outcome 1
#> predictor 4
#>
#> Operations:
#>
#> Saturation (asymptote= tune saturation_speed= tune Transformation onc(tv, search)
#producing the transformed dataset:
bake(rec_obj,mktdata)
#> # A tibble: 10 x 5
#> # Groups: prod, store [2]
#> prod store tv search sales
#> <fct> <fct> <dbl> <dbl> <dbl>
#> 1 brand store1 5.91 0 100
#> 2 brand store1 51.8 5.91 100
#> 3 brand store1 0 11.6 100
#> 4 brand store1 0 27.9 100
#> 5 brand store1 51.8 27.9 100
#> 6 brand store2 0 0 10
#> 7 brand store2 0 1.20 10
#> 8 brand store2 0 1.20 10
#> 9 brand store2 0 0 10
#> 10 brand store2 0 0 10
suppressMessages(suppressWarnings(library(tidyverse)))
suppressMessages(suppressWarnings(library(recipes)))
suppressMessages(suppressWarnings(library(tune)))
suppressMessages(suppressWarnings(library(dials)))
suppressMessages(suppressWarnings(library(mostlytidyMMM)))
#create two datasets:
mktdata<-rbind(tibble(prod='brand',store='store1',
sales=c(100.,100.,100.,100.,100.),
tv=c(10.,100.,0.,0.,100),
search=c(0,10,20,50.,50.),
week=c(1,2,3,4,5)) ,
tibble(prod='brand',store='store2',
sales=c(10.,10,10,10,10),tv=c(0.,0,0,0,0),
search=c(0.,2,2,0,0) ,
week=c(1,2,3,4,5)) ) |>
group_by(prod,store)
mktdata2<-tibble(prod='brand',store='all',sales=100,tv=1000,search=1000,week=6) |> group_by(prod,store)
#build a recipe with two different adstock steps -- could be done in one step
rec_obj <- recipe(sales ~ ., data = mktdata) |> step_adstock(tv,retention=.5,groups=c('prod','store')) |>
step_adstock(search,retention=.5,groups=c('prod','store'),time_id='week') |>
prep(training = mktdata)
#showing off the custom print function
print(rec_obj)
#> Recipe
#>
#> Inputs:
#>
#> role #variables
#> outcome 1
#> predictor 5
#>
#> Training data contained 10 data points and no missing data.
#>
#> Operations:
#>
#> Adstock Transformation with retention 0.5 ontv [trained]
#> Adstock Transformation with retention 0.5 onsearch [trained]
#see the final transformed output
bake(rec_obj,mktdata)
#> # A tibble: 10 x 6
#> # Groups: prod, store [2]
#> prod store tv search week sales
#> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
#> 1 brand store1 10 0 1 100
#> 2 brand store1 105 10 2 100
#> 3 brand store1 52.5 25 3 100
#> 4 brand store1 26.2 62.5 4 100
#> 5 brand store1 113. 81.2 5 100
#> 6 brand store2 0 0 1 10
#> 7 brand store2 0 2 2 10
#> 8 brand store2 0 3 3 10
#> 9 brand store2 0 1.5 4 10
#> 10 brand store2 0 0.75 5 10
#note what happens to store = 'all' in mktdata2:
bake(rec_obj,mktdata2)
#> Warning: There were 2 warnings in `summarise()`.
#> The first warning was:
#> i In argument: `min = min(delta, na.rm = T)`.
#> i In group 1: `prod = "brand"`, `store = "all"`.
#> Caused by warning in `min()`:
#> ! no non-missing arguments to min; returning Inf
#> i Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
#> There were 2 warnings in `summarise()`.
#> The first warning was:
#> i In argument: `min = min(delta, na.rm = T)`.
#> i In group 1: `prod = "brand"`, `store = "all"`.
#> Caused by warning in `min()`:
#> ! no non-missing arguments to min; returning Inf
#> i Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
#> # A tibble: 1 x 6
#> # Groups: prod, store [1]
#> prod store tv search week sales
#> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
#> 1 brand NA 1000 1000 6 100
#check to make sure that per-variable application of the step_adstock and step-saturation do not break the grouping structure even though
#the output tibble from bake is ungrouped (irritatingly)
rec_both_steps<-recipe(sales~.,data=mktdata) |>
step_adstock(tv,retention=.1) |> step_saturation(tv,asymptote=1000,saturation_speed=.001) |>
step_adstock(search,retention=.1,time_id='week') |> prep()
bake(rec_both_steps,mktdata)
#> # A tibble: 10 x 6
#> # Groups: prod, store [2]
#> prod store tv search week sales
#> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
#> 1 brand store1 9.95 0 1 100
#> 2 brand store1 96.1 10 2 100
#> 3 brand store1 10.0 21 3 100
#> 4 brand store1 1.01 52.1 4 100
#> 5 brand store1 95.3 55.2 5 100
#> 6 brand store2 0 0 1 10
#> 7 brand store2 0 2 2 10
#> 8 brand store2 0 2.2 3 10
#> 9 brand store2 0 0.22 4 10
#> 10 brand store2 0 0.022 5 10
#And adstock steps throw an error if time series isn't regular -- last point is two weeks
mktdata3<-rbind(tibble(prod='brand',store='store1',
sales=c(100.,100.,100.,100.,100.),
tv=c(10.,100.,0.,0.,100),
search=c(0,10,20,50.,50.),
week=c(1,2,3,4,6)) ,
tibble(prod='brand',store='store2',
sales=c(10.,10,10,10,10),tv=c(0.,0,0,0,0),
search=c(0.,2,2,0,0) ,
week=c(1,2,3,4,6)) ) |>
group_by(prod,store)
#not run to allow vignette to knit -- uncomment to see error
# rec1<-recipe(mktdata3) |> step_adstock(search,retention=.1,groups=c('prod','store'),
# time_id='week')|>prep()
#bake(rec1,mktdata3)It is not uncommon for a Marketing Mix Model to have dozens of variables. The functions allow roles and media transformations (ie the adstock and saturation steps) to be applied based on the configuration tables.
suppressMessages(suppressWarnings(library(tidyverse)))
suppressMessages(suppressWarnings(library(mostlytidyMMM)))
suppressMessages(suppressWarnings(library(recipes)))
suppressMessages(suppressWarnings(library(dials)))
suppressMessages(suppressWarnings(library(tune)))
#get control spreadsheet
control_file<-system.file('example model control.xlsx',package='mostlytidyMMM')
var_controls<-readxl::read_xlsx(control_file,'variables')
transform_controls<-readxl::read_xlsx(control_file,'role controls')
data_set=rename_columns_per_controls(read.csv(system.file("example2.csv", package = "mostlytidyMMM")),variable_controls=var_controls)
recipe0<-recipe(head(data_set,n=1) )
recipe1<-recipe0 |> bulk_update_role() |> bulk_add_role()
add_steps_media(recipe1)
#> Recipe
#>
#> Inputs:
#>
#> role #variables
#> group 2
#> outcome 1
#> postprocess 3
#> predictor 7
#> price 1
#> programmatic 2
#> time 1
#> time_id 1
#> trend 1
#> tv 4
#>
#> Operations:
#>
#> Adstock Transformation with retention tune on Adstock Transformation with retention programmatic_retention onhas_role("programmatic")
#> Saturation (asymptote= tune saturation_speed= tune Transformation on Saturation (asymptote= programmatic_asymptote saturation_speed= programmatic_saturation_speed Transformation onhas_role("programmatic")
#> Adstock Transformation with retention tune on Adstock Transformation with retention tv_retention onhas_role("tv")
#> Saturation (asymptote= tune saturation_speed= tune Transformation on Saturation (asymptote= tv_asymptote saturation_speed= tv_saturation_speed Transformation onhas_role("tv")
recipe1
#> Recipe
#>
#> Inputs:
#>
#> role #variables
#> group 2
#> outcome 1
#> postprocess 3
#> predictor 7
#> price 1
#> programmatic 2
#> time 1
#> time_id 1
#> trend 1
#> tv 4The ulam function generates a stan program and then runs sampling from it. It is very flexible, giving a user access to most of what stan can do. create_ulam_list() is a function that translates a lmer()-style formula into the requisite inputs for ulam.
suppressMessages(suppressWarnings(library(recipes)))
suppressMessages(suppressWarnings(library(tune)))
suppressMessages(suppressWarnings(library(dials)))
suppressMessages(suppressWarnings(library(tidyverse)))
suppressMessages(suppressWarnings(library(mostlytidyMMM)))
#get the control file:
control_file<-system.file('example model control.xlsx',package='mostlytidyMMM')
#get each relevant table of the control file:
var_controls<-readxl::read_xlsx(control_file,'variables')
transform_controls<-readxl::read_xlsx(control_file,'role controls')
workflow_controls<-readxl::read_xlsx(control_file,"workflow") |> select(-desc)
#pull data, add fourier transform columns
data1<-read.csv(system.file('example2.csv',package='mostlytidyMMM'))|>rename_columns_per_controls(variable_controls=var_controls)|>
rename_columns_per_controls()|> mutate(week=as.Date(week,"%m/%d/%Y"))|>
add_fourier_vars(vc=var_controls) |> add_groups_and_sort(vc=var_controls)
#append several more columns
lotsa_vars<-paste0('x',1:60)
new_x<-replicate(60,runif(nrow(data1)))
names(new_x)<-lotsa_vars
new_x<-as_tibble(new_x)
#> Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
#> `.name_repair` is omitted as of tibble 2.0.0.
#> i Using compatibility `.name_repair`.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#bring the data together:
data_for_lotsa_vars<-cbind(data1 ,new_x)
#make recipes
recipea<-recipe(head(data_for_lotsa_vars,n=1) )
recipeb<-recipea |> bulk_update_role() |> bulk_add_role()
recipeb<-recipeb |> add_steps_media() |> step_select(-has_role('postprocess'))
recipec <-recipeb |> step_mutate(week=as.numeric(week)-19247.65) |>
update_role(c(sin1,sin2,sin3,cos1,cos2,cos3),new_role='time') |>
add_role(c(sin1,sin2,sin3,cos1,cos2,cos3),new_role='predictor')
recipec<-recipec |> update_role(starts_with('V'),new_role='predictor')
# build a model formula in the lmer style from the recipe:
(formula_with_lotsa_vars<-create_formula(recipec))
#> [1] "sales ~ price + TV1 + TV2 + ProgVideo1 + trend + V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 + V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 + V21 + V22 + V23 + V24 + V25 + V26 + V27 + V28 + V29 + V30 + V31 + V32 + V33 + V34 + V35 + V36 + V37 + V38 + V39 + V40 + V41 + V42 + V43 + V44 + V45 + V46 + V47 + V48 + V49 + V50 + V51 + V52 + V53 + V54 + V55 + V56 + V57 + V58 + V59 + V60 + (sin1|store) + (cos1|store) + (1|store) + (TV1|store)"
#build a list suitable for input to rethinking::ulam from the formula and recipe:
(create_ulam_list(model_formula=formula_with_lotsa_vars))
#> [[1]]
#> sales ~ normal(big_model, big_sigma)
#> <environment: 0x0000000024b1ef00>
#>
#> [[2]]
#> big_model <- big_model_5 + a0 + store_int[store_id] + b_sin1_interact_store[store_id] +
#> b_cos1_interact_store[store_id] + b_TV1_interact_store[store_id]
#>
#> [[3]]
#> big_model_5 <- big_model_4 + b_V56 * V56 + b_V57 * V57 + b_V58 *
#> V58 + b_V59 * V59 + b_V60 * V60
#>
#> [[4]]
#> big_model_4 <- big_model_3 + b_V41 * V41 + b_V42 * V42 + b_V43 *
#> V43 + b_V44 * V44 + b_V45 * V45 + b_V46 * V46 + b_V47 * V47 +
#> b_V48 * V48 + b_V49 * V49 + b_V50 * V50 + b_V51 * V51 + b_V52 *
#> V52 + b_V53 * V53 + b_V54 * V54 + b_V55 * V55
#>
#> [[5]]
#> big_model_3 <- big_model_2 + b_V26 * V26 + b_V27 * V27 + b_V28 *
#> V28 + b_V29 * V29 + b_V30 * V30 + b_V31 * V31 + b_V32 * V32 +
#> b_V33 * V33 + b_V34 * V34 + b_V35 * V35 + b_V36 * V36 + b_V37 *
#> V37 + b_V38 * V38 + b_V39 * V39 + b_V40 * V40
#>
#> [[6]]
#> big_model_2 <- big_model_1 + b_V11 * V11 + b_V12 * V12 + b_V13 *
#> V13 + b_V14 * V14 + b_V15 * V15 + b_V16 * V16 + b_V17 * V17 +
#> b_V18 * V18 + b_V19 * V19 + b_V20 * V20 + b_V21 * V21 + b_V22 *
#> V22 + b_V23 * V23 + b_V24 * V24 + b_V25 * V25
#>
#> [[7]]
#> big_model_1 <- b_price * price + b_TV1 * TV1 + b_TV2 * TV2 +
#> b_ProgVideo1 * ProgVideo1 + b_trend * trend + b_V1 * V1 +
#> b_V2 * V2 + b_V3 * V3 + b_V4 * V4 + b_V5 * V5 + b_V6 * V6 +
#> b_V7 * V7 + b_V8 * V8 + b_V9 * V9 + b_V10 * V10
#>
#> $b_sin1_interact_store
#> b_sin1_interact_store[store_id] ~ normal(0, slope_sigma)
#> <environment: 0x0000000022f7a0b8>
#>
#> $b_cos1_interact_store
#> b_cos1_interact_store[store_id] ~ normal(0, slope_sigma)
#> <environment: 0x0000000022f7a0b8>
#>
#> $b_TV1_interact_store
#> b_TV1_interact_store[store_id] ~ normal(0, slope_sigma)
#> <environment: 0x0000000022f7a0b8>
#>
#> $store
#> store_int[store_id] ~ normal(65, int_sigma)
#> <environment: 0x0000000022f86eb8>
#>
#> $b_trend
#> b_trend ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V1
#> b_V1 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V2
#> b_V2 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V3
#> b_V3 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V4
#> b_V4 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V5
#> b_V5 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V6
#> b_V6 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V7
#> b_V7 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V8
#> b_V8 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V9
#> b_V9 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V10
#> b_V10 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V11
#> b_V11 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V12
#> b_V12 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V13
#> b_V13 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V14
#> b_V14 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V15
#> b_V15 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V16
#> b_V16 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V17
#> b_V17 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V18
#> b_V18 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V19
#> b_V19 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V20
#> b_V20 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V21
#> b_V21 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V22
#> b_V22 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V23
#> b_V23 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V24
#> b_V24 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V25
#> b_V25 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V26
#> b_V26 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V27
#> b_V27 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V28
#> b_V28 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V29
#> b_V29 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V30
#> b_V30 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V31
#> b_V31 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V32
#> b_V32 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V33
#> b_V33 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V34
#> b_V34 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V35
#> b_V35 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V36
#> b_V36 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V37
#> b_V37 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V38
#> b_V38 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V39
#> b_V39 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V40
#> b_V40 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V41
#> b_V41 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V42
#> b_V42 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V43
#> b_V43 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V44
#> b_V44 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V45
#> b_V45 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V46
#> b_V46 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V47
#> b_V47 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V48
#> b_V48 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V49
#> b_V49 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V50
#> b_V50 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V51
#> b_V51 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V52
#> b_V52 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V53
#> b_V53 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V54
#> b_V54 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V55
#> b_V55 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V56
#> b_V56 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V57
#> b_V57 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V58
#> b_V58 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V59
#> b_V59 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_V60
#> b_V60 ~ normal(0, 10)
#> <environment: 0x0000000022f73040>
#>
#> $b_price
#> b_price ~ normal(-1, 10)
#> <environment: 0x0000000023227d00>
#>
#> $b_TV1
#> b_TV1 ~ normal(1, 10)
#> <environment: 0x00000000231163e0>
#>
#> $b_TV2
#> b_TV2 ~ normal(1, 10)
#> <environment: 0x00000000230f7128>
#>
#> $b_ProgVideo1
#> b_ProgVideo1 ~ normal(1, 10)
#> <environment: 0x00000000230f2b78>
#>
#> [[77]]
#> a0 ~ normal(50, 50)
#> <environment: 0x0000000024b1ef00>
#>
#> [[78]]
#> big_sigma ~ half_cauchy(0, 100)
#> <environment: 0x0000000024b1ef00>
#>
#> [[79]]
#> int_sigma ~ half_cauchy(0, 10)
#> <environment: 0x0000000024b1ef00>
#>
#> [[80]]
#> slope_sigma ~ half_cauchy(0, 10)
#> <environment: 0x0000000024b1ef00>