Comment lisser les données et forcer la monotonie

14

J'ai quelques données que je voudrais lisser pour que les points lissés diminuent de façon monotone. Mes données diminuent fortement puis commencent à se stabiliser. Voici un exemple utilisant 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()

tracé des données à lisser

Quelle bonne technique de lissage pourrais-je utiliser? De plus, ce serait bien si je pouvais forcer le 1er point lissé à être proche de mon point observé.

Ben
la source
1
Je remarque que vos valeurs d'exemple sont entières. Vos vraies valeurs comptent-elles? S'ils l'étaient, alors (bien que ce ne soit pas une garantie de monotonie, pour des données comme celles-ci, ils le donneront généralement de toute façon), quelque chose comme ça pourrait être utile: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 similaire avec réponse: stats.stackexchange.com/questions/206073/…
kjetil b halvorsen

Réponses:

18

Vous pouvez le faire en utilisant des splines pénalisées avec des contraintes de monotonie via les fonctions mono.con()et pcls()dans le package mgcv . Il y a un peu à faire car ces fonctions ne sont pas aussi conviviales que gam(), mais les étapes sont indiquées ci-dessous, basées principalement sur l'exemple de ?pcls, modifiées pour s'adapter aux exemples de données que vous avez donnés:

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!

Maintenant, nous devons remplir l'objet qui est passé à pcls()contenant les détails du modèle contraint pénalisé que nous voulons adapter

## 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

Maintenant, nous pouvons enfin faire le montage

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

pcontient un vecteur de coefficients pour les fonctions de base correspondant à la spline. Pour visualiser la spline ajustée, nous pouvons prédire à partir du modèle à 100 emplacements sur la plage de x. Nous faisons 100 valeurs afin d'obtenir une belle ligne lisse sur l'intrigue.

## 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)))

Pour générer des valeurs prédites que nous utilisons Predict.matrix(), qui génère une matrice telle que lorsque plusieurs par coefficients pproduisent des valeurs prédites à partir du modèle ajusté:

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")

Cela produit:

entrez la description de l'image ici

Je vous laisse le soin de mettre les données sous une forme ordonnée pour tracer avec ggplot ...

Vous pouvez forcer un ajustement plus étroit (pour répondre partiellement à votre question sur l'ajustement plus fluide du premier point de données) en augmentant la dimension de la fonction de base de x. Par exemple, en définissant kégal à 8( k <- 8) et en réexécutant le code ci-dessus, nous obtenons

entrez la description de l'image ici

Vous ne pouvez pas pousser kbeaucoup plus haut pour ces données, et vous devez faire attention au sur-ajustement; tout ce qu'il pcls()y a à faire est de résoudre le problème des moindres carrés pénalisés étant donné les contraintes et les fonctions de base fournies, il n'effectue pas de sélection de lissage pour vous - pas que je sache ...)

Si vous voulez une interpolation, alors voyez la fonction R de base ?splinefunqui a des splines Hermite et des splines cubiques avec des contraintes de monotonie. Dans ce cas, vous ne pouvez toutefois pas l'utiliser, car les données ne sont pas strictement monotones.

Réintégrer Monica - G. Simpson
la source
Merci. Je suis sûr que votre solution est appropriée, mais elle est si complexe et obscurcie que je ne peux tout simplement pas l'utiliser. splinefunétait ma pensée initiale aussi (j'interpole) mais spline(x=df$x, y=df$y, n=nrow(df), method="monoH.FC")et les spline(x=df$x, y=df$y, n=nrow(df), method="hyman")deux soulèvent des erreurs
Ben
1
Si vous essayez, je suis sûr que vous pouvez l'utiliser; Je n'ai aucune idée de ce qui se passe sous le capot ici, mais j'ai travaillé, et j'ai indiqué les endroits où vous auriez besoin de changer les choses. En supposant que vous connaissiez certains R bien sûr . La plupart des détails sont de mise en œuvre que vous pouvez ignorer si tout ce que vous voulez faire s’adapte à une spline contrainte monotone. Souhaitez-vous que j'annote un peu plus le code pour souligner davantage ce que fait chaque étape? La référence en ?mono.cona plus de détails sur la méthode.
Reinstate Monica - G. Simpson
Quant à savoir pourquoi splinefunsoulève une erreur; Je viens de réaliser que vous pouvez adapter une spline monotone qui interpole des données qui ne sont pas elles-mêmes monotones. L'observation à x = 6est plus grande yque les observations à x = 5. Vous n'aurez qu'à ignorer cette partie de la réponse :-)
Réintégrer Monica - G. Simpson
Je l'ai. Et pas besoin - je suis un utilisateur R assez expérimenté. J'aime juste comprendre les mathématiques derrière ce que j'utilise et cette solution semble avoir beaucoup de choses sous le capot. Merci encore pour votre aide.
Ben
J'ai ajouté quelques notes pour expliquer ce que chaque chose est ou fait; le point principal à noter est que les contraintes de monotonie sont imposées par un ensemble spécifique de contraintes d'inégalité qui mono.conrevient pour une spline cubique. ?pclscontient des exemples de splines à plaques minces et de modèles additifs qui sont moins conviviaux que ceux ci-dessus, mais qui pourraient exposer un peu plus les mathématiques si vous êtes familier avec les mathématiques de ces types de splines (je ne suis pas très familier moi-même).
Reinstate Monica - G. Simpson
13

Le récent paquet d'arnaques de Natalya Pya et basé sur le document "Shape constrained additive models" de Pya & Wood (2015) peut faciliter une partie du processus mentionné dans l'excellente réponse de Gavin.

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

Il existe un certain nombre de fonctions bs que vous pouvez utiliser - dans ce qui précède, j'ai utilisé mpd pour "monotonique décroissante P-spline", mais il a également des fonctions qui imposent la convexité ou la concavité séparément ou à côté des contraintes monotones.

srepho
la source