Como suavizar dados e forçar a monotonicidade

14

Eu tenho alguns dados que gostaria de suavizar para que os pontos suavizados diminuam monotonicamente. Meus dados diminuem acentuadamente e depois começam a se espalhar. Aqui está um exemplo usando R

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))
ggplot(df, aes(x=x, y=y))+geom_line()

plotagem de dados para suavizar

O que é uma boa técnica de suavização que eu poderia usar? Além disso, seria bom poder forçar o 1º ponto suavizado a ficar próximo ao ponto observado.

Ben
fonte
1
Percebo que seus valores de exemplo são inteiros. Seus valores reais são importantes? Se assim fosse, então (enquanto isso não é garantia de monotonicidade, para dados como estes geralmente vão dar-lhe de qualquer maneira), algo como isso poderia ser útil:plot(y~x,data=df); f=fitted( glm( y~ns(x,df=4), data=df,family=quasipoisson)); lines(df$x,f)
Glen_b -Reinstate Monica
Q semelhante com a resposta: stats.stackexchange.com/questions/206073/…
kjetil b halvorsen

Respostas:

18

Você pode fazer isso usando splines penalizados com restrições de monotonicidade através das funções mono.con()e pcls()no pacote mgcv . Há um pouco de brincadeira por fazer, porque essas funções não são tão amigáveis ​​quanto as do usuário gam(), mas as etapas são mostradas abaixo, com base principalmente no exemplo de ?pcls, modificado para se adequar aos dados de amostra que você forneceu:

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))

## Set up the size of the basis functions/number of knots
k <- 5
## This fits the unconstrained model but gets us smoothness parameters that
## that we will need later
unc <- gam(y ~ s(x, k = k, bs = "cr"), data = df)

## This creates the cubic spline basis functions of `x`
## It returns an object containing the penalty matrix for the spline
## among other things; see ?smooth.construct for description of each
## element in the returned object
sm <- smoothCon(s(x, k = k, bs = "cr"), df, knots = NULL)[[1]]

## This gets the constraint matrix and constraint vector that imposes
## linear constraints to enforce montonicity on a cubic regression spline
## the key thing you need to change is `up`.
## `up = TRUE` == increasing function
## `up = FALSE` == decreasing function (as per your example)
## `xp` is a vector of knot locations that we get back from smoothCon
F <- mono.con(sm$xp, up = FALSE)   # get constraints: up = FALSE == Decreasing constraint!

Agora precisamos preencher o objeto que é passado para pcls()conter detalhes do modelo restrito penalizado que queremos ajustar

## Fill in G, the object pcsl needs to fit; this is just what `pcls` says it needs:
## X is the model matrix (of the basis functions)
## C is the identifiability constraints - no constraints needed here
##   for the single smooth
## sp are the smoothness parameters from the unconstrained GAM
## p/xp are the knot locations again, but negated for a decreasing function
## y is the response data
## w are weights and this is fancy code for a vector of 1s of length(y)
G <- list(X = sm$X, C = matrix(0,0,0), sp = unc$sp,
          p = -sm$xp, # note the - here! This is for decreasing fits!
      y = df$y,
          w = df$y*0+1)
G$Ain <- F$A    # the monotonicity constraint matrix
G$bin <- F$b    # the monotonicity constraint vector, both from mono.con
G$S <- sm$S     # the penalty matrix for the cubic spline
G$off <- 0      # location of offsets in the penalty matrix

Agora podemos finalmente fazer o encaixe

## Do the constrained fit 
p <- pcls(G)  # fit spline (using s.p. from unconstrained fit)

pcontém um vetor de coeficientes para as funções básicas correspondentes ao spline. Para visualizar o spline ajustado, podemos prever a partir do modelo em 100 locais no intervalo de x. Fazemos 100 valores para obter uma boa linha suave no gráfico.

## predict at 100 locations over range of x - get a smooth line on the plot
newx <- with(df, data.frame(x = seq(min(x), max(x), length = 100)))

Para gerar valores previstos, usamos Predict.matrix(), que gera uma matriz que, quando múltiplos por coeficientes, pproduzem valores previstos a partir do modelo ajustado:

fv <- Predict.matrix(sm, newx) %*% p
newx <- transform(newx, yhat = fv[,1])

plot(y ~ x, data = df, pch = 16)
lines(yhat ~ x, data = newx, col = "red")

Isso produz:

insira a descrição da imagem aqui

Vou deixar que você obtenha os dados em um formato organizado para plotar com o ggplot ...

Você pode forçar um ajuste mais próximo (para responder parcialmente à sua pergunta sobre como ajustar o mais suave no primeiro ponto de dados) aumentando a dimensão da função base de x. Por exemplo, definir kigual a 8( k <- 8) e executar novamente o código acima, obtemos

insira a descrição da imagem aqui

Você não pode aumentar kmuito esses dados e precisa ter cuidado com o excesso de ajuste; tudo o que você pcls()está fazendo é resolver o problema dos mínimos quadrados penalizados, dadas as restrições e as funções básicas fornecidas, não está executando a seleção de suavidade para você - não que eu saiba ...)

Se você deseja interpolação, veja a função R básica ?splinefunque possui splines Hermite e splines cúbicos com restrições de monotonicidade. Nesse caso, você não pode usá-lo, pois os dados não são estritamente monotônicos.

Restabelecer Monica - G. Simpson
fonte
Obrigado. Tenho certeza de que sua solução é apropriada, mas é tão complexa e ofuscada que não consigo usá-la. splinefunfoi o meu pensamento inicial também (eu estou interpolando), mas spline(x=df$x, y=df$y, n=nrow(df), method="monoH.FC")e spline(x=df$x, y=df$y, n=nrow(df), method="hyman")ambos geram erros #
1955 Ben Ben
1
Se você apenas tentar, tenho certeza que você pode usá-lo; Eu tenho pouca ideia do que está acontecendo aqui, mas resolvi e indiquei os lugares que você precisaria para mudar as coisas. Supondo que você conheça R, é claro . A maioria dos detalhes é implementacional que você pode ignorar se tudo o que você deseja fazer se encaixa em um spline monotonicamente restrito. Gostaria que eu anotasse um pouco mais o código para destacar mais sobre o que cada etapa está fazendo? A referência em ?mono.contem mais detalhes sobre o método.
Reintegrar Monica - G. Simpson
Quanto ao porquê splinefunestá levantando um erro; Acabei de perceber que você pode ajustar um spline monotônico que interpola dados que não são monotônicos. A observação em x = 6é maior yque as observações em x = 5. Você apenas terá que ignorar essa parte da resposta :-)
Restabelece Monica - G. Simpson
Entendi. E não há necessidade - eu sou um usuário R bastante experiente. Eu só gosto de entender a matemática por trás do que eu uso e essa solução parece ter muita coisa acontecendo por baixo do capô. Obrigado novamente por sua ajuda.
19416 Ben
Eu adicionei algumas notas para explicar o que cada coisa é ou faz; o ponto principal a ser observado é que as restrições de monotonicidade estão sendo impostas por um conjunto específico de restrições de desigualdade que mono.conretorna para um spline cúbico. ?pclstem exemplos de splines de placa fina e modelos aditivos que são menos amigáveis ​​que o anterior, mas que podem expor um pouco mais da matemática se você estiver familiarizado com a matemática para esses tipos de spline (eu também não sou familiar).
Reintegrar Monica - G. Simpson
13

O recente pacote de fraude de Natalya Pya e baseado no artigo "Modelos de aditivos com restrição de forma" de Pya & Wood (2015) pode tornar muito mais fácil parte do processo mencionado na excelente resposta de Gavin.

library(scam)
con <- scam(y ~ s(x, k = k, bs = "mpd"), data = df)
plot(con)

Existem várias funções bs que você pode usar - no exemplo acima, usei mpd para "P-spline decrescente monotônico", mas também possui funções que reforçam a convexidade ou concavidade, separadamente ou ao lado das restrições monotônicas.

srepho
fonte