Modelando dados censurados com propensabilidade

Nada é perfeito, e os dados também não. Um tipo de “imperfeição” é dados ausentesonde alguns recursos não estão observados para alguns assuntos. (Um tópico para outro submit.) Outra é dados censuradosonde um evento cujas características queremos medir não ocorre no intervalo de observação. O exemplo em Richard McElreath’s Repensar estatística é hora de adoção de gatos em um abrigo de animais. Se corrigirmos um intervalo e observarmos os tempos de espera para os gatos que realmente fez Ser adotado, nossa estimativa acabará muito otimista: não levamos em consideração os gatos que não foram adotados durante esse intervalo e, portanto, teriam contribuído com tempos de espera de comprimento por mais tempo que o intervalo completo.

Neste submit, usamos um exemplo um pouco menos emocional que, no entanto, pode ser de interesse, especialmente para os desenvolvedores de pacotes R: tempo para a conclusão de R CMD verifycoletado de Cran e fornecido pelo parsnip Pacote AS check_times. Aqui, a parte censurada são aquelas verificações que erraram por qualquer motivo, ou seja, para o qual o cheque não foi concluído.

Por que nos preocupamos com a parte censurada? No cenário de adoção de gatos, isso é bastante óbvio: queremos poder obter uma estimativa realista para qualquer gato desconhecido, não apenas para os gatos que serão “sortudos”. Que tal check_times? Bem, se o seu envio é um daqueles que erraram, você ainda se importa com quanto tempo espera, portanto, mesmo que a porcentagem deles seja baixa (<1%), não queremos simplesmente excluí -los. Além disso, existe a possibilidade de que os falhos tenham demorado mais, se tivessem conclusão, devido a alguma diferença intrínseca entre os dois grupos. Por outro lado, se as falhas fossem aleatórias, os cheques de longo prazo teriam uma probability maior de ser atingido por um erro. Então, aqui também, exluir os dados censurados pode resultar em viés.

Como podemos modelar as durações para essa porção censurada, onde a “duração verdadeira” é desconhecida? Dando um passo para trás, como podemos modelar durações em geral? Fazendo o menor número possível de suposições, o distribuição máxima de entropia Para deslocamentos (no espaço ou no tempo), é exponencial. Assim, para as verificações que realmente concluíram, supõe -se que as durações sejam distribuídas exponencialmente.

Para os outros, tudo o que sabemos é que, em um mundo digital onde o cheque concluído, levaria pelo menos tanto tempo como a duração fornecida. Essa quantidade pode ser modelada pela função de distribuição cumulativa complementar exponencial (CCDF). Por que? Uma função de distribuição cumulativa (CDF) indica a probabilidade de que um valor menor ou igual a algum ponto de referência tenha sido atingido; por exemplo, “a probabilidade de durações <= 255 é 0,9”. Seu complemento, 1 - CDF, fornece a probabilidade de um valor exceder que esse ponto de referência.

Vamos ver isso em ação.

Os dados

O código a seguir funciona com as versões estáveis ​​atuais da probabilidade Tensorflow e Tensorflow, que são 1,14 e 0,7, respectivamente. Se você não tiver tfprobability Instalado, obtenha -o no GitHub:

Estas são as bibliotecas de que precisamos. A partir do tensorflow 1.14, chamamos tf$compat$v2$enable_v2_behavior() para correr com execução ansiosa.

Além das durações de cheques que queremos modelar, check_times relata vários recursos do pacote em questão, como número de pacotes importados, número de dependências, tamanho de código e arquivos de documentação and many others. O standing A variável indica se a verificação concluída ou errada.

df <- check_times %>% choose(-bundle)
glimpse(df)

Observations: 13,626
Variables: 24
$ authors         1, 1, 1, 1, 5, 3, 2, 1, 4, 6, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1,…
$ imports         0, 6, 0, 0, 3, 1, 0, 4, 0, 7, 0, 0, 0, 0, 3, 2, 14, 2, 2, 0…
$ suggests        2, 4, 0, 0, 2, 0, 2, 2, 0, 0, 2, 8, 0, 0, 2, 0, 1, 3, 0, 0,…
$ relies upon         3, 1, 6, 1, 1, 1, 5, 0, 1, 1, 6, 5, 0, 0, 0, 1, 1, 5, 0, 2,…
$ Roxygen         0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0,…
$ gh              0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0,…
$ rforge          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ descr           217, 313, 269, 63, 223, 1031, 135, 344, 204, 335, 104, 163,…
$ r_count         2, 20, 8, 0, 10, 10, 16, 3, 6, 14, 16, 4, 1, 1, 11, 5, 7, 1…
$ r_size          0.029053, 0.046336, 0.078374, 0.000000, 0.019080, 0.032607,…
$ ns_import       3, 15, 6, 0, 4, 5, 0, 4, 2, 10, 5, 6, 1, 0, 2, 2, 1, 11, 0,…
$ ns_export       0, 19, 0, 0, 10, 0, 0, 2, 0, 9, 3, 4, 0, 1, 10, 0, 16, 0, 2…
$ s3_methods      3, 0, 11, 0, 0, 0, 0, 2, 0, 23, 0, 0, 2, 5, 0, 4, 0, 0, 0, …
$ s4_methods      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ doc_count       0, 3, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
$ doc_size        0.000000, 0.019757, 0.038281, 0.000000, 0.007874, 0.000000,…
$ src_count       0, 0, 0, 0, 0, 0, 0, 2, 0, 5, 3, 0, 0, 0, 0, 0, 0, 54, 0, 0…
$ src_size        0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,…
$ data_count      2, 0, 0, 3, 3, 1, 10, 0, 4, 2, 2, 146, 0, 0, 0, 0, 0, 10, 0…
$ data_size       0.025292, 0.000000, 0.000000, 4.885864, 4.595504, 0.006500,…
$ testthat_count  0, 8, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, 0, 0,…
$ testthat_size   0.000000, 0.002496, 0.000000, 0.000000, 0.000000, 0.000000,…
$ check_time      49, 101, 292, 21, 103, 46, 78, 91, 47, 196, 200, 169, 45, 2…
$ standing          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

Destas 13.626 observações, apenas 103 são censuradas:

0     1 
103 13523 

Para uma melhor legibilidade, trabalharemos com um subconjunto das colunas. Nós usamos surv_reg Para nos ajudar a encontrar um subconjunto útil e interessante de preditores:

survreg_fit <-
  surv_reg(dist = "exponential") %>% 
  set_engine("survreg") %>% 
  match(Surv(check_time, standing) ~ ., 
      knowledge = df)
tidy(survreg_fit) 

# A tibble: 23 x 7
   time period             estimate std.error statistic  p.worth conf.low conf.excessive
                                         
 1 (Intercept)     3.86      0.0219     176.     0.             NA        NA
 2 authors         0.0139    0.00580      2.40   1.65e- 2       NA        NA
 3 imports         0.0606    0.00290     20.9    7.49e-97       NA        NA
 4 suggests        0.0332    0.00358      9.28   1.73e-20       NA        NA
 5 relies upon         0.118     0.00617     19.1    5.66e-81       NA        NA
 6 Roxygen         0.0702    0.0209       3.36   7.87e- 4       NA        NA
 7 gh              0.00898   0.0217       0.414  6.79e- 1       NA        NA
 8 rforge          0.0232    0.0662       0.351  7.26e- 1       NA        NA
 9 descr           0.000138  0.0000337    4.10   4.18e- 5       NA        NA
10 r_count         0.00209   0.000525     3.98   7.03e- 5       NA        NA
11 r_size          0.481     0.0819       5.87   4.28e- 9       NA        NA
12 ns_import       0.00352   0.000896     3.93   8.48e- 5       NA        NA
13 ns_export      -0.00161   0.000308    -5.24   1.57e- 7       NA        NA
14 s3_methods      0.000449  0.000421     1.06   2.87e- 1       NA        NA
15 s4_methods     -0.00154   0.00206     -0.745  4.56e- 1       NA        NA
16 doc_count       0.0739    0.0117       6.33   2.44e-10       NA        NA
17 doc_size        2.86      0.517        5.54   3.08e- 8       NA        NA
18 src_count       0.0122    0.00127      9.58   9.96e-22       NA        NA
19 src_size       -0.0242    0.0181      -1.34   1.82e- 1       NA        NA
20 data_count      0.0000415 0.000980     0.0423 9.66e- 1       NA        NA
21 data_size       0.0217    0.0135       1.61   1.08e- 1       NA        NA
22 testthat_count -0.000128  0.00127     -0.101  9.20e- 1       NA        NA
23 testthat_size   0.0108    0.0139       0.774  4.39e- 1       NA        NA

Parece que se escolhermos importsAssim, relies uponAssim, r_sizeAssim, doc_sizeAssim, ns_import e ns_export Acabamos com uma mistura de preditores (comparativamente) poderosos de diferentes espaços semânticos e de diferentes escalas.

Antes de podar o quadro de dados, economizamos a variável de destino. Em nosso modelo e configuração de treinamento, é conveniente ter dados censurados e sem censura armazenados separadamente, então aqui criamos dois Matrizes alvo em vez de uma:

# verify occasions for failed checks
# _c stands for censored
check_time_c <- df %>%
  filter(standing == 0) %>%
  choose(check_time) %>%
  as.matrix()

# verify occasions for profitable checks 
check_time_nc <- df %>%
  filter(standing == 1) %>%
  choose(check_time) %>%
  as.matrix()

Agora, podemos aumentar o zoom das variáveis ​​de interesse, configurando um DataFrame para os dados censurados e outro para os dados sem censura cada. Todos os preditores são normalizados para evitar transbordamento durante a amostragem. Adicionamos uma coluna de 1s para uso como interceptação.

df <- df %>% choose(standing,
                    relies upon,
                    imports,
                    doc_size,
                    r_size,
                    ns_import,
                    ns_export) %>%
  mutate_at(.vars = 2:7, .funs = operate(x) (x - min(x))/(max(x)-min(x))) %>%
  add_column(intercept = rep(1, nrow(df)), .earlier than = 1)

# dataframe of predictors for censored knowledge  
df_c <- df %>% filter(standing == 0) %>% choose(-standing)
# dataframe of predictors for non-censored knowledge 
df_nc <- df %>% filter(standing == 1) %>% choose(-standing)

É isso para os preparativos. Mas é claro que estamos curiosos. Os tempos de verificação parecem diferentes? Os preditores – os que escolhemos – parecem diferentes?

Comparando alguns percentis significativos para ambas as lessons, vemos que as durações para verificações incompletas são mais altas do que as para verificações concluídas, além do percentil 100%. Não é de surpreender que, dada a enorme diferença no tamanho da amostra, a duração máxima é maior para verificações concluídas. Caso contrário, porém, não parece que as verificações de pacote erradas “iriam demorar mais”?

Que tal os preditores? Não vemos nenhuma diferença para relies upono número de dependências de pacotes (além de, novamente, o máximo mais alto alcançado para pacotes cuja verificação concluída):

Mas para todos os outros, vemos o mesmo padrão relatado acima para check_time. O número de pacotes importados é maior para dados censurados em todos os percentis, além do máximo:

O mesmo para ns_exporto número estimado de funções ou métodos exportados:

Bem como para ns_importo número estimado de funções ou métodos importados:

Mesmo padrão para r_sizeo tamanho do disco de arquivos no R diretório:

E finalmente, vemos isso para doc_size também, onde doc_size é do tamanho de .Rmd e .Rnw arquivos:

Dada a nossa tarefa em questão – o modelo verifique as durações, levando em consideração os dados sem censura, bem como os dados censurados – não nos debruçaremos mais sobre as diferenças entre os dois grupos; No entanto, achamos interessante relacionar esses números.

Então agora, de volta ao trabalho. Precisamos criar um modelo.

O modelo

Conforme explicado na introdução, para verificação concluída, a duração é modelada usando um PDF exponencial. Isso é tão simples quanto adicionar tfd_exponencial () para a função do modelo, tfd_joint_distribution_sequencial (). Para a parte censurada, precisamos do CCDF exponencial. Este não é, a partir de hoje, facilmente adicionado ao modelo. O que podemos fazer é calcular seu valor e adicioná -lo à probabilidade “principal” do modelo. Veremos isso abaixo ao discutir a amostragem; Por enquanto, significa que a definição do modelo termina direta, pois cobre apenas os dados não censurados. É feito apenas do referido PDF e anteriores exponenciais para os parâmetros de regressão.

Quanto ao último, usamos os anteriores gaussianos de 0 centros para todos os parâmetros. Desvios padrão de 1 acabaram funcionando bem. Como os anteriores são todos iguais, em vez de listar um monte de tfd_normalS, podemos criá -los todos de uma vez como

tfd_sample_distribution(tfd_normal(0, 1), sample_shape = 7)

O tempo médio de verificação é modelado como uma combinação afim dos seis preditores e da interceptação. Aqui está o modelo completo, instanciado usando apenas os dados sem censura:

mannequin <- operate(knowledge) {
  tfd_joint_distribution_sequential(
    record(
      tfd_sample_distribution(tfd_normal(0, 1), sample_shape = 7),
      operate(betas)
        tfd_independent(
          tfd_exponential(
            price = 1 / tf$math$exp(tf$transpose(
              tf$matmul(tf$solid(knowledge, betas$dtype), tf$transpose(betas))))),
          reinterpreted_batch_ndims = 1)))
}

m <- mannequin(df_nc %>% as.matrix())

Sempre, testamos se as amostras desse modelo têm as formas esperadas:

samples <- m %>% tfd_sample(2)
samples

((1))
tf.Tensor(
(( 1.4184642   0.17583323 -0.06547955 -0.2512014   0.1862184  -1.2662812
   1.0231884 )
 (-0.52142304 -1.0036682   2.2664437   1.29737     1.1123234   0.3810004
   0.1663677 )), form=(2, 7), dtype=float32)

((2))
tf.Tensor(
((4.4954767  7.865639   1.8388556  ... 7.914391   2.8485563  3.859719  )
 (1.549662   0.77833986 0.10015647 ... 0.40323067 3.42171    0.69368565)), form=(2, 13523), dtype=float32)

Parece bom: temos uma lista do comprimento dois, um elemento para cada distribuição no modelo. Para ambos os tensores, a dimensão 1 reflete o tamanho do lote (que definimos arbitrariamente como 2 neste teste), enquanto a dimensão 2 é 7 para o número de anteriores normais e 13523 para o número de durações previstas.

Qual é a probabilidade dessas amostras?

m %>% tfd_log_prob(samples)

tf.Tensor((-32464.521   -7693.4023), form=(2,), dtype=float32)

Aqui também, a forma está correta e os valores parecem razoáveis.

A próxima coisa a fazer é definir o alvo que queremos otimizar.

Alvo de otimização

Resumo, a coisa a maximizar é a probabilidade de log dos dados – ou seja, as durações medidas – sob o modelo. Agora, aqui os dados vêm em duas partes e o alvo também. Primeiro, temos os dados não censurados, para os quais

m %>% tfd_log_prob(record(betas, tf$solid(target_nc, betas$dtype)))

calculará a probabilidade de log. Segundo, para obter a probabilidade de log para os dados censurados, escrevemos uma função personalizada que calcula o log do CCDF exponencial:

get_exponential_lccdf <- operate(betas, knowledge, goal) {
  e <-  tfd_independent(tfd_exponential(price = 1 / tf$math$exp(tf$transpose(tf$matmul(
    tf$solid(knowledge, betas$dtype), tf$transpose(betas)
  )))),
  reinterpreted_batch_ndims = 1)
  cum_prob <- e %>% tfd_cdf(tf$solid(goal, betas$dtype))
  tf$math$log(1 - cum_prob)
}

Ambas as peças são combinadas em uma pequena função de invólucro que nos permite comparar o treinamento, incluindo e excluindo os dados censurados. Não faremos isso nesta postagem, mas você pode estar interessado em fazê -lo com seus próprios dados, especialmente se a proporção de peças censuradas e sem censura for um pouco menos desequilibrada.

get_log_prob <-
  operate(target_nc,
           censored_data = NULL,
           target_c = NULL) {
    log_prob <- operate(betas) {
      log_prob <-
        m %>% tfd_log_prob(record(betas, tf$solid(target_nc, betas$dtype)))
      potential <-
        if (!is.null(censored_data) && !is.null(target_c))
          get_exponential_lccdf(betas, censored_data, target_c)
      else
        0
      log_prob + potential
    }
    log_prob
  }

log_prob <-
  get_log_prob(
    check_time_nc %>% tf$transpose(),
    df_c %>% as.matrix(),
    check_time_c %>% tf$transpose()
  )

Amostragem

Com o modelo e o destino definidos, estamos prontos para fazer amostragem.

n_chains <- 4
n_burnin <- 1000
n_steps <- 1000

# hold observe of some diagnostic output, acceptance and step dimension
trace_fn <- operate(state, pkr) {
  record(
    pkr$inner_results$is_accepted,
    pkr$inner_results$accepted_results$step_size
  )
}

# get form of preliminary values 
# to start out sampling with out producing NaNs, we'll feed the algorithm
# tf$zeros_like(initial_betas)
# as a substitute 
initial_betas <- (m %>% tfd_sample(n_chains))((1))

Para o número de etapas do LeapFrog e o tamanho da etapa, a experimentação mostrou que uma combinação de 64 / 0.1 produziu resultados razoáveis:

hmc <- mcmc_hamiltonian_monte_carlo(
  target_log_prob_fn = log_prob,
  num_leapfrog_steps = 64,
  step_size = 0.1
) %>%
  mcmc_simple_step_size_adaptation(target_accept_prob = 0.8,
                                   num_adaptation_steps = n_burnin)

run_mcmc <- operate(kernel) {
  kernel %>% mcmc_sample_chain(
    num_results = n_steps,
    num_burnin_steps = n_burnin,
    current_state = tf$ones_like(initial_betas),
    trace_fn = trace_fn
  )
}

# vital for efficiency: run HMC in graph mode
run_mcmc <- tf_function(run_mcmc)

res <- hmc %>% run_mcmc()
samples <- res$all_states

Resultados

Antes de inspecionarmos as correntes, aqui está uma rápida olhada na proporção de etapas aceitas e no tamanho médio da etapa por parâmetro:

0.995

0.004953894

Também armazenamos tamanhos de amostra eficazes e o rhat métricas para adição posterior à sinopse.

effective_sample_size <- mcmc_effective_sample_size(samples) %>%
  as.matrix() %>%
  apply(2, imply)
potential_scale_reduction <- mcmc_potential_scale_reduction(samples) %>%
  as.numeric()

Nós então convertemos o samples Tensor a uma matriz R para uso no pós -processamento.

# 2-item record, the place every merchandise has dim (1000, 4)
samples <- as.array(samples) %>% array_branch(margin = 3)

Quão bem a amostra funcionou? As correntes se misturam bem, mas para alguns parâmetros, a autocorrelação ainda é bastante alta.

prep_tibble <- operate(samples) {
  as_tibble(samples,
            .name_repair = ~ c("chain_1", "chain_2", "chain_3", "chain_4")) %>%
    add_column(pattern = 1:n_steps) %>%
    collect(key = "chain", worth = "worth",-pattern)
}

plot_trace <- operate(samples) {
  prep_tibble(samples) %>%
    ggplot(aes(x = pattern, y = worth, colour = chain)) +
    geom_line() +
    theme_light() +
    theme(
      legend.place = "none",
      axis.title = element_blank(),
      axis.textual content = element_blank(),
      axis.ticks = element_blank()
    )
}

plot_traces <- operate(samples) {
  plots <- purrr::map(samples, plot_trace)
  do.name(grid.prepare, plots)
}

plot_traces(samples)

Figura 1: Gráficos de rastreamento para os 7 parâmetros.

Agora, para uma sinopse das estatísticas de parâmetros posteriores, incluindo os indicadores usuais de amostragem por parâmetro Tamanho eficaz da amostra e rhat.

all_samples <- map(samples, as.vector)

means <- map_dbl(all_samples, imply)

sds <- map_dbl(all_samples, sd)

hpdis <- map(all_samples, ~ hdi(.x) %>% t() %>% as_tibble())

abstract <- tibble(
  imply = means,
  sd = sds,
  hpdi = hpdis
) %>% unnest() %>%
  add_column(param = colnames(df_c), .after = FALSE) %>%
  add_column(
    n_effective = effective_sample_size,
    rhat = potential_scale_reduction
  )

abstract

# A tibble: 7 x 7
  param       imply     sd  decrease higher n_effective  rhat
                     
1 intercept  4.05  0.0158  4.02   4.08       508.   1.17
2 relies upon    1.34  0.0732  1.18   1.47      1000    1.00
3 imports    2.89  0.121   2.65   3.12      1000    1.00
4 doc_size   6.18  0.394   5.40   6.94       177.   1.01
5 r_size     2.93  0.266   2.42   3.46       289.   1.00
6 ns_import  1.54  0.274   0.987  2.06       387.   1.00
7 ns_export -0.237 0.675  -1.53   1.10        66.8  1.01

Figura 2: Meios posteriores e HPDIs.

A partir das parcelas de diagnóstico e rastreamento, o modelo parece funcionar razoavelmente bem, mas como não há uma métrica de erro direto envolvido, é difícil saber se as previsões reais chegariam a um intervalo apropriado.

Para garantir que sim, inspecionamos previsões de nosso modelo e também de surv_reg. Desta vez, também dividimos os dados em conjuntos de treinamento e teste. Aqui estão primeiro as previsões de surv_reg:

train_test_split <- initial_split(check_times, strata = "standing")
check_time_train <- coaching(train_test_split)
check_time_test <- testing(train_test_split)

survreg_fit <-
  surv_reg(dist = "exponential") %>% 
  set_engine("survreg") %>% 
  match(Surv(check_time, standing) ~ relies upon + imports + doc_size + r_size + 
        ns_import + ns_export, 
      knowledge = check_time_train)
survreg_fit(sr_fit)

# A tibble: 7 x 7
  time period         estimate std.error statistic  p.worth conf.low conf.excessive
                                    
1 (Intercept)  4.05      0.0174     234.    0.             NA        NA
2 relies upon      0.108     0.00701     15.4   3.40e-53       NA        NA
3 imports      0.0660    0.00327     20.2   1.09e-90       NA        NA
4 doc_size     7.76      0.543       14.3   2.24e-46       NA        NA
5 r_size       0.812     0.0889       9.13  6.94e-20       NA        NA
6 ns_import    0.00501   0.00103      4.85  1.22e- 6       NA        NA
7 ns_export   -0.000212  0.000375    -0.566 5.71e- 1       NA        NA

survreg_pred <- 
  predict(survreg_fit, check_time_test) %>% 
  bind_cols(check_time_test %>% choose(check_time, standing))  

ggplot(survreg_pred, aes(x = check_time, y = .pred, colour = issue(standing))) +
  geom_point() + 
  coord_cartesian(ylim = c(0, 1400))

Figura 3: Previsões do conjunto de testes de Surv_reg. Um outlier (do valor 160421) é excluído by way of coord_carttesian () para evitar distorcer o gráfico.

Para o modelo MCMC, reiniciamos apenas o conjunto de treinamento e obtemos o resumo do parâmetro. O código é análogo ao acima e não é mostrado aqui.

Agora podemos prever no conjunto de testes, para simplificar apenas o uso dos meios posteriores:

df <- check_time_test %>% choose(
                    relies upon,
                    imports,
                    doc_size,
                    r_size,
                    ns_import,
                    ns_export) %>%
  add_column(intercept = rep(1, nrow(check_time_test)), .earlier than = 1)

mcmc_pred <- df %>% as.matrix() %*% abstract$imply %>% exp() %>% as.numeric()
mcmc_pred <- check_time_test %>% choose(check_time, standing) %>%
  add_column(.pred = mcmc_pred)

ggplot(mcmc_pred, aes(x = check_time, y = .pred, colour = issue(standing))) +
  geom_point() + 
  coord_cartesian(ylim = c(0, 1400)) 

Figura 4: Previsões do conjunto de testes do modelo MCMC. Não há outliers, apenas usando a mesma escala acima para comparação.

Isso parece bom!

Envolver

Mostramos como modelar dados censurados – ou melhor, um subtipo frequente dele envolvendo durações – usando tfprobability. O check_times dados de parsnip Foi uma escolha divertida, mas essa técnica de modelagem pode ser ainda mais útil quando a censura é mais substancial. Espero que seu submit tenha fornecido algumas orientações sobre como lidar com dados censurados em seu próprio trabalho. Obrigado pela leitura!