Tracer une ligne de régression par morceaux

10

Existe-t-il un moyen de tracer la ligne de régression d'un modèle par morceaux comme celui-ci, autre que d'utiliser linespour tracer chaque segment séparément ou en utilisant geom_smooth(aes(group=Ind), method="lm", fill=FALSE)?

m.sqft <- mean(sqft)
model <- lm(price~sqft+I((sqft-m.sqft)*Ind))
# sqft, price: continuous variables, Ind: if sqft>mean(sqft) then 1 else 0

plot(sqft,price)
abline(reg = model)
Warning message:
In abline(reg = model) :
  only using the first two of 3regression coefficients

Je vous remercie.

George Dontas
la source

Réponses:

6

La seule façon dont je sais comment le faire facilement est de prédire à partir du modèle sur toute la plage sqftet de tracer les prédictions. Il n'y a pas de manière générale avec ablineou similaire. Vous pouvez également jeter un œil au package segmenté qui s'adaptera à ces modèles et fournira l'infrastructure de traçage pour vous.

Faire cela via des prédictions et des graphiques de base. Tout d'abord, quelques données factices:

set.seed(1)
sqft <- runif(100)
sqft <- ifelse((tmp <- sqft > mean(sqft)), 1, 0) + rnorm(100, sd = 0.5)
price <- 2 + 2.5 * sqft
price <- ifelse(tmp, price, 0) + rnorm(100, sd = 0.6)
DF <- data.frame(sqft = sqft, price = price,
                 Ind = ifelse(sqft > mean(sqft), 1, 0))
rm(price, sqft)
plot(price ~ sqft, data = DF)

Monter le modèle:

mod <- lm(price~sqft+I((sqft-mean(sqft))*Ind), data = DF)

Générez des données pour prévoir et prévoir:

m.sqft <- with(DF, mean(sqft))
pDF <- with(DF, data.frame(sqft = seq(min(sqft), max(sqft), length = 200)))
pDF <- within(pDF, Ind <- ifelse(sqft > m.sqft, 1, 0))
pDF <- within(pDF, price <- predict(mod, newdata = pDF))

Tracez les lignes de régression:

ylim <- range(pDF$price, DF$price)
xlim <- range(pDF$sqft, DF$sqft)
plot(price ~ sqft, data = DF, ylim = ylim, xlim = xlim)
lines(price ~ sqft, data = pDF, subset = Ind > 0, col = "red", lwd = 2)
lines(price ~ sqft, data = pDF, subset = Ind < 1, col = "red", lwd = 2)

Vous pouvez coder cela en une fonction simple - vous n'avez besoin que des étapes des deux blocs de code précédents - que vous pouvez utiliser à la place de abline:

myabline <- function(model, data, ...) {
    m.sqft <- with(data, mean(sqft))
    pDF <- with(data, data.frame(sqft = seq(min(sqft), max(sqft),
                                            length = 200)))
    pDF <- within(pDF, Ind <- ifelse(sqft > m.sqft, 1, 0))
    pDF <- within(pDF, price <- predict(mod, newdata = pDF))
    lines(price ~ sqft, data = pDF, subset = Ind > 0, ...)
    lines(price ~ sqft, data = pDF, subset = Ind < 1, ...)
    invisible(model)
}

Alors:

ylim <- range(pDF$price, DF$price)
xlim <- range(pDF$sqft, DF$sqft)
plot(price ~ sqft, data = DF, ylim = ylim, xlim = xlim)
myabline(mod, DF, col = "red", lwd = 2)

Via le package segmenté

require(segmented)
mod2 <- lm(price ~ sqft, data = DF)
mod.s <- segmented(mod2, seg.Z = ~ sqft, psi = 0.5,
                   control = seg.control(stop.if.error = FALSE))
plot(price ~ sqft, data = DF)
plot(mod.s, add = TRUE)
lines(mod.s, col = "red")

Avec ces données, il n'évalue pas le point d'arrêt à mean(sqft), mais les méthodes plotet linesdans ce package peuvent vous aider à implémenter quelque chose de plus générique que myablinede faire ce travail directement pour vous à partir du lm()modèle ajusté .

Modifier: si vous souhaitez segmenté pour estimer l'emplacement du point d'arrêt, définissez l' 'psi'argument sur NA:

mod.s <- segmented(mod2, seg.Z = ~ sqft, psi = NA,
                   control = seg.control(stop.if.error = FALSE))

Puis segmentedessaiera les K = 10quantiles de sqft, avecK réglage seg.control()et qui sont par défaut 10. Voir ?seg.controlpour en savoir plus.

Gavin Simpson
la source
@Gavin (+1) Réponse beaucoup plus complète que la mienne; Je l'aime juste.
chl
@Gavin La section "Via le package segmenté" n'a pas fonctionné pour mes données. J'ai obtenu un "Aucun point d'arrêt estimé" après avoir exécuté la segmentedcommande.
George Dontas
@ gd047: Excuses, il y a eu une erreur dans le code que j'ai montré. Vous devez fournir à l'argument seq.Zune formule unilatérale des variables qui ont une relation segmentée avec la réponse. J'ai modifié ma réponse pour inclure seq.Z = ~ sqftet ajouter une note sur le segmentedchoix des valeurs de psipour vous.
Gavin Simpson
@ gd047 Je voudrais supprimer ma réponse car celle-ci répond mieux à votre question d'origine. Cela vous dérangerait d'accepter celui-ci au lieu du mien?
chl
moel<-mF:unergumentjesnotjenterpretunebleuneslogjecuneljenunejetjeon:Wunernjengmessunege:jenjeF(moel)objF