Extrapolando com splines B e GAMs

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Do fundo da pilha – R, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Um problema que costuma surgir durante a modelagem com modelos aditivos gerados (GAMs), especialmente com séries temporais ou dados espaciais, é como extrapolar além do alcance dos dados usados ​​para treinar o modelo? O problema surge porque os GAMs usam splines para aprender com os dados usando funções básicas. Os splines em si são construídos a partir de funções básicas que normalmente são configuradas em termos dos dados usados ​​para ajustar o modelo. Se não houver funções básicas além do intervalo dos dados de entrada, o que exatamente está sendo usado se queremos extrapolar? Uma questão relacionada é a da penalidade de wiggliness; dependendo do tipo de base usado, a penalidade pode se estender por toda a linha real (-∞ – ∞) ou apenas pelo intervalo dos dados de entrada. Neste post, quero dar uma olhada prática no comportamento de extrapolação de splines em GAMs equipados com o mgcv pacote para R. Em particular, quero ilustrar o quão flexível é a base do spline B.

Muito do que discuto na postagem se baseia muito na página de ajuda do mgcv para a base da spline B – ?mgcv::b.spline – e uma recente discussão por e-mail com Alex Hayes, Dave Miller e Eric Pedersen, embora o que escrevo aqui reflita minha própria opinião sobre essa discussão.

Inicialmente, pensei em analisar isso novamente depois de ler uma nova pré-impressão sobre aproximações de baixa patente de um processo gaussiano (GP; Riutort-Mayol et al., 2020), onde, entre outras coisas, os autores comparam o comportamento do modelo exato do GP com sua versão baixa e com uma spline de regressão em chapa fina (TPRS). O TPRS é o tipo de coisa que você obteria por padrão com mgcv e s(), mas como os outros modelos eram totalmente bayesianos, o modelo TPRS foi ajustado usando brm() de brms pacote para que todos os modelos fossem comparáveis, sendo finalmente Stan. O modelo TPRS não fez um bom trabalho ao ajustar as observações de teste ao extrapolar além dos limites dos dados. Gostaria de saber se poderíamos fazer melhor com a base do spline B em mgcv como eu sabia que tinha flexibilidade extra para extrapolação curta além dos dados, mas nunca havia realmente analisado como funcionava ou qual era o respectivo comportamento.

Leia Também  Ciência de dados colaborativa: orientação de alto nível para análises éticas de pares científicos

Se você deseja recriar elementos do restante da postagem, precisará dos seguintes pacotes instalados:

## Packages
library('ggplot2')
library('tibble')
library('tidyr')
library('dplyr')
library('mgcv')
library('gratia')
library('patchwork')

## remotes::install_github("clauswilke/colorblindr")
library('colorblindr')
## remotes::install_github("clauswilke/relayer")
library('relayer')

Os dois últimos são usados ​​para plotagem e o relayer um pacote em particular é necessário, pois vou usar duas escalas de cores separadas nos gráficos. Se você não os tiver instalado, poderá instalá-los usando o Remotos pacote e o código nas linhas comentadas acima.

O conjunto de dados de exemplo usado na comparação foi postado no repositório GitHub da pré-impressão, por isso foi fácil pegá-los e começar a brincar. Para carregar os dados no R, podemos usar

load(url("https://bit.ly/gprocdata"))
ls()
[1] "f_true"

onde o link curto Bitly apenas vincula ao .Rdata arquivo armazenado no GitHub. Isso cria um objeto, f_true, na área de trabalho. Veremos a verdadeira função em um minuto. Após a pré-impressão, um conjunto de dados de observações ruidosas é simulado a partir da função true adicionando ruído gaussiano (μ = 0, σ = 0,2)

seed <- 1234
set.seed(seed)
gp_data <- tibble(truth = unname(f_true), x = seq(-1, 1, by = 0.002)) %>%
    mutate(y = truth + rnorm(length(truth), , 0.2))

A partir desse conjunto barulhento, amostramos 250 observações aleatoriamente e indicamos algumas delas como estando em um conjunto de testes que não usaremos ao ajustar os GAMs

set.seed(seed)
r_samp <- sample_n(gp_data, size = 250) %>%
    arrange(x) %>%
    mutate(data_set = case_when(x < -0.8 ~ "test",
                                x > 0.8 ~ "test",
                                x > -0.45 & x < -0.36 ~ "test",
                                x > -0.05 & x < 0.05 ~ "test",
                                x > 0.45 & x < 0.6 ~ "test",
                                TRUE ~ "train"))

Finalmente, visualizamos a verdadeira função e as observações barulhentas que obtivemos dela

Leia Também  Crie e converta rabiscos | R-bloggers

ggplot(r_samp, aes(x = x, y = y, colour = data_set)) +
    geom_line(aes(y = truth, colour = NULL), show.legend = FALSE, alpha = 0.5) +
    geom_point() +
    scale_colour_brewer(palette = "Set1", name = "Data set")

A verdadeira função e as observações barulhentas extraídas dela. Os pontos azuis são as observações de treinamento que usaremos para ajustar os modelos, enquanto os pontos vermelhos são observações de teste usadas para investigar como os modelos interpolam e extrapolam.

Os pontos vermelhos são as observações do teste e serão usados ​​para analisar o comportamento dos splines sob condições de interpolação e extrapolação.

Splines de placa fina

Primeiro, veremos como as estrias de chapa fina se comportam sob extrapolação, recriando o comportamento da pré-impressão. Começo ajustando dois GAMs nos quais usamos 50 funções básicas (k = 50) a partir da base TPRS (bs = “tp”) O argumento m controla a ordem da penalidade derivada; o padrão é m = 2, para uma segunda penalidade derivada (penalizando a curvatura do spline). Para o segundo modelo, usamos m = 1, indicando uma penalidade na primeira derivada do TPRS, que penaliza desvios de uma função plana. Observe que filtramos a amostra de dados ruidosos para incluir apenas as observações do treinamento.

m_tprs2 <- gam(y ~ s(x, k = 50, bs = "tp", m = 2),
               data = filter(r_samp, data_set == "train"), method = "REML")
## first order penalty
m_tprs1 <- gam(y ~ s(x, k = 50, bs = "tp", m = 1),
               data = filter(r_samp, data_set == "train"), method = "REML")

Não vou me preocupar em analisar o diagnóstico de modelo nesta postagem e, em vez disso, pular para a análise de como esses dois modelos se comportam quando previmos além dos limites dos dados de treinamento.

Em seguida, defino algumas novas observações para prever a partir dos dois modelos

new_data <- tibble(x = seq(-1.5, 1.5, by = 0.002))

Lembre-se de que os dados de treinamento cobriram o intervalo de -0,8 a 0,8; portanto, estamos extrapolando bastante proporcionalmente do suporte dos dados de treinamento. Agora podemos prever a partir dos dois modelos

p_tprs2 <- as_tibble(predict(m_tprs2, new_data, se.fit = TRUE)) %>%
    rename(fit_tprs_2 = fit, se_tprs_2 = se.fit)
p_tprs1 <- as_tibble(predict(m_tprs1, new_data, se.fit = TRUE)) %>%
    rename(fit_tprs_1 = fit, se_tprs_1 = se.fit)

Observe que nomeamos as duas colunas de dados com algumas informações necessárias para a plotagem, portanto, os sublinhados são importantes.

Em seguida, organizamos alguns dados para obter as previsões em um formato organizado e adequado para plotagem

crit <- qnorm((1 - 0.89) / 2, lower.tail = FALSE)
new_data_tprs <- bind_cols(new_data, p_tprs2, p_tprs1) %>%
    pivot_longer(fit_tprs_2:se_tprs_1, names_sep = '_',
                 names_to = c('variable', 'spline', 'order')) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    mutate(upr_ci = fit + (crit * se), lwr_ci = fit - (crit * se))

A idéia básica aqui é que lançamos os dados em uma versão longa e fina muito geral e extraímos variáveis ​​indicando o tipo de valor (fit = montado e se = erro padrão), o tipo de spline e a ordem da penalidade, dividindo o sublinhado em cada uma das colunas de entrada. Em seguida, convertemos o quadro de dados longo e fino em uma versão um pouco mais ampla, na qual temos acesso ao fit e se variáveis, antes de calcular um intervalo credível de 89% nos valores previstos.

Agora podemos plotar os dados mais os valores previstos

ggplot(mapping = aes(x = x, y = y)) +
    geom_ribbon(data = new_data_tprs,
                mapping = aes(ymin = lwr_ci, ymax = upr_ci, x = x,
                              fill = order),
                inherit.aes = FALSE, alpha = 0.2) +
    geom_point(data = r_samp, aes(colour = data_set)) +
    geom_line(data = new_data_tprs, aes(y = fit, x = x, colour2 = order),
              size = 1) %>%
    rename_geom_aes(new_aes = c("colour" = "colour2")) +
    scale_colour_brewer(palette = "Set1", aesthetics = "colour",
                        name = "Data set") +
    scale_colour_OkabeIto(aesthetics = "colour2", name = "Penalty") +
    scale_fill_OkabeIto(name = "Penalty") +
    coord_cartesian(ylim = c(-2, 2)) +
    labs(title = "Extrapolating with thin plate splines",
         subtitle = "How behaviour varies with derivative penalties of different order")

Médias preditivas posteriores para os dois modelos de spline de regressão de placa fina mostrando o comportamento de interpolação e extrapolação com penalidades de primeira e segunda derivadas.

Com a penalidade padrão da segunda derivada, vemos que, sob extrapolação, o spline exibe comportamento linear. Para o primeiro modelo de penalização deriavtive, o comportamento é prever um valor constante. Os intervalos credíveis também são irrealisticamente estreitos no caso do modelo TPRS com a primeira penalidade derivada. Também não é um trabalho particularmente bom estimar qualquer amostra de teste fora da faixa de x nos dados de treinamento. Os modelos se saem melhor ao interpolar, exceto pela seção ao redor x = 0,5.

B splines

ESTÁ BEM. E as estrias B? Com o construtor B spline em mgcv temos muito controle sobre como estabelecemos a base e a penalidade de wiggliness. Veremos mais dessas opções posteriormente, mas primeiro, veremos o comportamento padrão em que a penalidade opera apenas no intervalo das observações de treinamento.

m_bs_default <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, 2)),
                    data = filter(r_samp, data_set == "train"), method = "REML")
Warning in smooth.construct.bs.smooth.spec(object, dk$data, dk$knots): there is
*no* information about some basis coefficients

Aqui, solicitamos um spline B cúbico com uma penalidade de segunda ordem – esse é o spline B cúbico comum ou de jardim, onde a penalidade de wigglines cobre o intervalo de x nos dados de treinamento. Ignore o aviso; isso ocorre apenas porque temos muitas funções e algumas não são suportadas por nenhum dos dados devido aos buracos causados ​​pelas observações do teste.

Se queremos que a penalidade se estenda um pouco além da faixa de x, precisamos passar um conjunto de pontos finais sobre os quais os nós serão definidos. Precisamos especificar os dois pontos extremos extremos que abrangem a região que queremos prever e dois nós internos que cobrem o intervalo dos dados, além de um pouco. Nós especificamos esses nós abaixo

knots <- list(x = c(-2, -0.9, 0.9, 2))

e depois passar knots ao knots argumento ao ajustar o modelo

m_bs_extrap <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, 2)), method = "REML",
                   data = filter(r_samp, data_set == "train"), knots = knots)
Warning in smooth.construct.bs.smooth.spec(object, dk$data, dk$knots): there is
*no* information about some basis coefficients

A única diferença aqui é como especificamos que queremos que a penalidade se estenda para além dos limites das observações de treinamento. Você receberá outro aviso aqui. Isso sempre acontece quando você define nós externos além do alcance dos dados; é inofensivo.

Podemos visualizar as diferenças nas bases usando basis() de gratia pacote

bs_default <- basis(s(x, k = 50, bs = "bs", m = c(3, 2)), knots = knots,
                    data = filter(new_data, x >= -0.8 & x <= 0.8))
bs_extrap <- basis(s(x, k = 50, bs = "bs", m = c(3, 2)), knots = knots,
                   data = new_data)
lims <- lims(x = c(-1.5, 1.5))
vlines <- geom_vline(data = tibble(x = c(-0.8, 0.8)),
                     aes(xintercept = x), lty = "dashed")
(draw(bs_default) + lims + vlines) / (draw(bs_extrap) + lims + vlines) +
    plot_annotation(tag_levels = 'A')

Bases de spline B cúbicas com nós cobrindo o intervalo de observações de treinamento (A) e com nós externos cobrindo o intervalo de dados de treinamento mais a região onde queremos extrapolar. Usar os nós externos tem o efeito de estender a penalidade de peruca sobre a região que queremos prever. As linhas tracejadas são desenhadas em x = -0,8 e x = 0,8, os limites das observações de treinamento.

Tecnicamente, as funções básicas no painel superior se estenderiam um pouco para a região de previsão, mas basis() ainda não consegue lidar com o uso de um conjunto de dados para configurar a base e outro no qual avaliar. Como temos funções básicas que se estendem ao longo do intervalo de previsão, a penalidade de wiggliness também pode ser aplicada nessa região.

Agora, prevemos a partir de ambos os modelos como antes e repetimos a disputa de dados

p_bs_default <- as_tibble(predict(m_bs_default, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_default = fit, se_bs_default = se.fit)
p_bs_extrap <- as_tibble(predict(m_bs_extrap, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_extrap = fit, se_bs_extrap = se.fit)

new_data_bs_eg <- bind_cols(new_data, p_bs_default, p_bs_extrap) %>%
    pivot_longer(fit_bs_default:se_bs_extrap, names_sep = '_',
                 names_to = c('variable', 'spline', 'penalty')) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    mutate(upr_ci = fit + (crit * se), lwr_ci = fit - (crit * se))

A única diferença aqui é que eu codifiquei nos nomes das variáveis ​​se usamos a penalidade padrão ou a que foi estendida além dos limites dos dados. Traçamos os ajustes com

ggplot(mapping = aes(x = x, y = y)) +
    geom_ribbon(data = new_data_bs_eg,
                mapping = aes(ymin = lwr_ci, ymax = upr_ci, x = x, fill = penalty),
                inherit.aes = FALSE, alpha = 0.2) +
    geom_point(data = r_samp, aes(colour = data_set)) +
    geom_line(data = new_data_bs_eg, aes(y = fit, x = x, colour2 = penalty),
              size = 1) %>%
    rename_geom_aes(new_aes = c("colour" = "colour2")) +
    scale_colour_brewer(palette = "Set1", aesthetics = "colour", name = "Data set") +
    scale_colour_OkabeIto(aesthetics = "colour2", name = "Penalty") +
    scale_fill_OkabeIto(name = "Penalty") +
    coord_cartesian(ylim = c(-2, 2)) +
    labs(title = "Extrapolating with B splines",
         subtitle = "How behaviour varies when the penalty extends beyond the data")

Meios preditivos posteriores para os dois modelos de spline B que mostram o comportamento de interpolação e extrapolação quando a penalidade cobre o intervalo de dados e quando se estende além desse intervalo.

Como esses dois modelos usaram penalidades de segunda derivada, ambos extrapolam linearmente além do alcance das observações de treinamento. É importante ressaltar, porém, que temos um comportamento muito diferente dos intervalos credíveis, especialmente no final mais baixo da x, em que o amplo intervalo é uma melhor representação da incerteza que temos nas previsões extrapoladas. Esse é um comportamento melhor, pois pelo menos estamos sendo honestos sobre a incerteza ao extrapolar.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

Comparando bases diferentes

Até agora, tão desinteressante. Antes de chegarmos às coisas boas e demonstrar outros recursos da base do spline B em mgcv, vamos comparar rapidamente os modelos de spline TPRS e B com um processo gaussiano suave projetado para corresponder de perto à função de geração de dados. Observe que este GP é instalado usando mgcv onde precisamos especificar a escala de comprimento e, como tal, não deve ser diretamente comparável com os modelos GP exatos ou com classificação baixa de Riutort-Mayol et al. (2020).

No mgcv um GP pode ser ajustado usando bs = “gp”. Quando fazemos isso, o significado do m argumento muda. Aqui, estamos solicitando uma função de covariância de Matérn com ν = 3/2 e escala de comprimento de 0,15. Esses valores foram escolhidos para corresponder aos da função verdadeira.

m_gp <- gam(y ~ s(x, k = 50, bs = "gp", m = c(3, 0.15)),
            data = filter(r_samp, data_set == "train"), method = "REML")

Novamente, temos algumas disputas a fazer para juntar tudo isso em um objeto que podemos traçar facilmente

p_bs <- as_tibble(predict(m_bs_extrap, new_data, se.fit = TRUE)) %>%
    rename(fit_bs = fit, se_bs = se.fit)
p_tprs <- as_tibble(predict(m_tprs2, new_data, se.fit = TRUE)) %>%
    rename(fit_tprs = fit, se_tprs = se.fit)
p_gp <- as_tibble(predict(m_gp, new_data, se.fit = TRUE)) %>%
    rename(fit_gp = fit, se_gp = se.fit)

new_data_bases <- bind_cols(new_data, p_tprs, p_bs, p_gp) %>%
    pivot_longer(fit_tprs:se_gp, names_sep = '_',
                 names_to = c('variable', 'spline')) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    mutate(upr_ci = fit + (2 * se), lwr_ci = fit - (2 * se))

E finalmente traçamos usando

ggplot(mapping = aes(x = x, y = y)) +
    geom_ribbon(data = new_data_bases,
                mapping = aes(ymin = lwr_ci, ymax = upr_ci, x = x, fill = spline),
                inherit.aes = FALSE, alpha = 0.2) +
    geom_point(data = r_samp, aes(colour = data_set)) +
    geom_line(data = new_data_bases, aes(y = fit, x = x, colour2 = spline),
              size = 1) %>%
    rename_geom_aes(new_aes = c("colour" = "colour2")) +
    scale_colour_brewer(palette = "Set1", aesthetics = "colour", name = "Data set") +
    scale_colour_OkabeIto(aesthetics = "colour2", name = "Basis") +
    scale_fill_OkabeIto(name = "Basis") +
    coord_cartesian(ylim = c(-2, 2)) +
    labs(title = "Extrapolating with splines",
         subtitle = "How behaviour varies with different basis types")
Warning: Ignoring unknown aesthetics: colour2

Médias preditivas posteriores para três GAMs; um spline de placa fina com penalidade da 2ª derivada, um spline B com penalidade da 2ª derivada estendido no intervalo para previsão e um processo gaussiano com uma função de covariância de Matérn (ν = 3/2) com escala de comprimento = 0,15

Claramente, o GP se aproxima dos dados do teste ao extrapolar, mas essa não é uma comparação justa, pois eu disse ao modelo qual era a escala de comprimento correta! Poderíamos tentar estimar isso a partir dos dados, ajustando modelos sobre uma grade de valores prováveis ​​para o parâmetro de escala de comprimento e usando o modelo com a menor pontuação REML, mas não mostrarei como fazer isso aqui; Eu tenho um código de exemplo nos suplementos para Simpson (2018) mostrando como fazer isso, você está interessado.

Mais com splines B

Não estamos restritos a usar a penalidade da segunda derivada com splines B; podemos usar penalidades de terceira, segunda, primeira ou até zero ordem com splines B cúbicos. Como o comportamento deles varia ao interpolar e extrapolar?

Por conveniência, ajustarei os três modelos a um formato comum, mesmo que já tenhamos visto e ajustado o primeiro modelo com a penalidade da segunda derivada. Observe como especificamos a ordem da penalidade derivada passando um segundo valor ao argumento m; m = 1 é uma primeira penalidade derivada, m = 0 uma penalidade derivada de zero, etc.

m_bs_2 <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, 2)), method = "REML",
              data = filter(r_samp, data_set == "train"), knots = knots)
m_bs_1 <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, 1)), method = "REML",
              data = filter(r_samp, data_set == "train"), knots = knots)
m_bs_0 <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, )), method = "REML",
              data = filter(r_samp, data_set == "train"), knots = knots)

Mais uma vez, repetimos a disputa de dados para obter algo que possamos traçar

p_bs_2 <- as_tibble(predict(m_bs_2, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_2 = fit, se_bs_2 = se.fit)
p_bs_1 <- as_tibble(predict(m_bs_1, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_1 = fit, se_bs_1 = se.fit)
p_bs_0 <- as_tibble(predict(m_bs_0, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_0 = fit, se_bs_0 = se.fit)

new_data_order <- bind_cols(new_data, p_bs_2, p_bs_1, p_bs_0) %>%
    pivot_longer(fit_bs_2:se_bs_0, names_sep = '_',
                 names_to = c('variable', 'spline', 'order')) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    mutate(upr_ci = fit + (2 * se), lwr_ci = fit - (2 * se))

Observe novamente como estou definindo os nomes das colunas que contêm valores ajustados e seus erros padrão para facilitar a retirada desses dados durante o pivot_longer() degrau.

Traçamos os valores previstos com

ggplot(mapping = aes(x = x, y = y)) +
    geom_ribbon(data = new_data_order,
                mapping = aes(ymin = lwr_ci, ymax = upr_ci, x = x, fill = order),
                inherit.aes = FALSE, alpha = 0.2) +
    geom_point(data = r_samp, aes(colour = data_set)) +
    geom_line(data = new_data_order, aes(y = fit, x = x, colour2 = order),
              size = 1) %>%
    rename_geom_aes(new_aes = c("colour" = "colour2")) +
    scale_colour_brewer(palette = "Set1", aesthetics = "colour", name = "Data set") +
    scale_colour_OkabeIto(aesthetics = "colour2", name = "Penalty") +
    scale_fill_OkabeIto(name = "Penalty") +
    coord_cartesian(ylim = c(-2, 2)) +
    labs(title = "Extrapolating with B splines",
         subtitle = "How behaviour varies with penalties of different order")

Meios preditivos posteriores para três GAMs usando splines B com diferentes ordens de penalidade derivativa, todos cobrindo a região onde queremos prever as amostras de teste; uma spline B com penalidade de 2ª derivada, uma spline B com penalidade de 1ª derivada e uma spline B com penalidade de derivada zero.

O enredo mostra as diferentes penalidades, levando a uma ampla gama de comportamentos. O spline com a penalidade de ordem zero interpola mal, aparentemente indo para a média geral dos dados durante cada uma das seções de teste dentro do intervalo de x. Ao extrapolar, vemos novamente esse comportamento de “reversão média”, o que significa que ele se sai bem ao extrapolar para grandes valores de x, mas se sai extremamente mal na extremidade inferior do x. Os intervalos credíveis para este modelo também são irrealisticamente estreitos, como os do modelo TPRS com a penalidade da 1ª derivada que vimos anteriormente.

O modelo com a primeira penalidade derivada tem um comportamento razoável; extrapola como uma função amplamente plana, continuando a partir dos valores mínimo e máximo de x, como o TPRS se encaixa com uma primeira penalidade derivada que vimos acima, mas os intervalos credíveis são muito mais realistas para o spline B do que para o TPRS. Observe também que os intervalos do spline B com a primeira penalidade derivada não explodem tão rapidamente quanto os do spline B se ajustam à penalidade do segundo derivado.

Penalidades múltiplas

Um truque final que a base da spline B mgcv A vantagem é que você pode combinar várias penalidades em um único spline. Poderíamos aplicar splines B cúbicos com uma, duas, três ou até quatro penalidades. As penalidades adicionais são especificadas passando mais valores para m: m = c(3, 2, 1) seria um spline B cúbico com penalidade de uma segunda derivada e uma primeira derivada, enquanto m = c(3, 2, 1, 0) você obteria um spline cúbico com todas as três penalidades. Você pode misturar e combinar o quanto quiser, com algumas exceções:

  • você pode ter apenas uma penalidade para cada pedido; portanto, não é possível penalizar uma derivada de maneira mais rigorosa, adicionando mais de uma penalidade a ela; m = c(3, 2, 2, 1) por exemplo não é permitido e

  • você só pode ter valores para m[i] (Onde i > 1) que existem para a ordem dada do spline B, ou seja, onde m[i] ≤ m[1].

No código abaixo, eu encaixo dois modelos adicionais com misturas de penalidades e os comparo com a penalidade de segunda derivada padrão (ajustada anteriormente). Em cada caso, estou novamente usando o knots argumento para estender as penalidades ao longo do intervalo que podemos querer prever.

m_bs_21 <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, 2, 1)), method = "REML",
                data = filter(r_samp, data_set == "train"), knots = knots)
Warning in smooth.construct.bs.smooth.spec(object, dk$data, dk$knots): there is
*no* information about some basis coefficients
m_bs_210 <- gam(y ~ s(x, k = 50, bs = "bs", m = c(3, 2, 1, )), method = "REML",
                data = filter(r_samp, data_set == "train"), knots = knots)
Warning in smooth.construct.bs.smooth.spec(object, dk$data, dk$knots): there is
*no* information about some basis coefficients

Novamente, fazemos a mesma disputa, desta vez codificando as misturas de pedidos nos nomes das colunas

p_bs_21 <- as_tibble(predict(m_bs_21, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_21 = fit, se_bs_21 = se.fit)
p_bs_210 <- as_tibble(predict(m_bs_210, new_data, se.fit = TRUE)) %>%
    rename(fit_bs_210 = fit, se_bs_210 = se.fit)

new_data_multi <- bind_cols(new_data, p_bs_2, p_bs_21, p_bs_210) %>%
    pivot_longer(fit_bs_2:se_bs_210, names_sep = '_',
                 names_to = c('variable', 'spline', 'order')) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    mutate(upr_ci = fit + (2 * se), lwr_ci = fit - (2 * se),
           penalty = case_when(order == "2" ~ "2",
                               order == "21" ~ "2, 1",
                               order == "210" ~ "2, 1, 0"))

O último passo aqui usa case_when() para escrever uma formatação melhor para as penalidades, para obter uma legenda melhor na trama, que produzimos com

ggplot(mapping = aes(x = x, y = y)) +
    geom_ribbon(data = new_data_multi,
                mapping = aes(ymin = lwr_ci, ymax = upr_ci, x = x, fill = penalty),
                inherit.aes = FALSE, alpha = 0.2) +
    geom_point(data = r_samp, aes(colour = data_set)) +
    geom_line(data = new_data_multi, aes(y = fit, x = x, colour2 = penalty),
              size = 1) %>%
    rename_geom_aes(new_aes = c("colour" = "colour2")) +
    scale_colour_brewer(palette = "Set1", aesthetics = "colour", name = "Data set") +
    scale_colour_OkabeIto(aesthetics = "colour2", name = "Penalty") +
    scale_fill_OkabeIto(name = "Penalty") +
    coord_cartesian(ylim = c(-2, 2)) +
    labs(title = "Extrapolating with B splines",
         subtitle = "How behaviour changes when combining multiple penalties")

Médias preditivas posteriores para três GAMs usando splines B com misturas de penalidades derivadas, todas cobrindo a região onde queremos prever as amostras de teste; um spline B com penalidade única de 2ª derivada, um spline B com penalidade de 2ª e 1ª derivada e um spline B com 2nd1st e 0º penalidades derivativas.

Ao misturar as penalidades, misturamos alguns dos recursos de comportamento. Por exemplo, o estranho comportamento de interpolação do spline B com penalidade derivada de zero é essencialmente removido quando combinado com penalidade de segunda e primeira derivada.

Dados os dados, as previsões que essencialmente preveem funções constantes além do alcance dos dados, mas com intervalos amplos e credíveis, são provavelmente as mais realistas; em cada caso em que usamos uma spline B que incluía uma primeira penalidade derivada, pelo menos, cobriu a maior parte da observação de teste além da faixa de x.

No entanto, em nenhum dos ataques obtemos comportamento que se aproxime de ajustar as observações do teste além do treinamento de x nos dados de treinamento, mesmo ao usar um processo gaussiano que supostamente corresponde a pelo menos a forma geral da função verdadeira.

Referências

Riutort-Mayol, G., Bürkner, P.-C., Andersen, M.R., Solin, A. e Vehtari, A. (2020). O espaço prático de hilbert aproxima processos gaussianos bayesianos para programação probabilística. Disponível em: http://arxiv.org/abs/2004.11408.

Simpson, G.L. (2018). Modelagem de séries temporais paleoecológicas usando modelos aditivos generalizados. Fronteiras em Ecologia e Evolução 6, 149. doi: 10.3389 / fevo.2018.00149.

var vglnk = {key: ‘949efb41171ac6ec1bf7f206d57e90b8’};

(função (d, t) {
var s = d.createElement
s.src = ‘//cdn.viglink.com/api/vglnk.js’;
var r = d.getElementsByTagName
} (documento, ‘script’));

Para Deixe um comentário para o autor, siga o link e comente no blog: Do fundo da pilha – R.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br