Rysowanie linii regresji częściowej

10

Czy istnieje sposób wykreślenia linii regresji takiego fragmentowego modelu, inny niż użycie linesdo wykreślenia każdego segmentu osobno lub użycie 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

Dziękuję Ci.

George Dontas
źródło

Odpowiedzi:

6

Jedynym sposobem, w jaki wiem, jak to łatwo zrobić, jest przewidywanie z modelu w całym zakresie sqfti wykreślanie prognoz. Nie ma ogólnego sposobu z ablinelub podobnym. Możesz także spojrzeć na segmentowany pakiet, który będzie pasował do tych modeli i zapewni infrastrukturę drukowania.

Robi to za pomocą prognoz i podstawowej grafiki. Po pierwsze, niektóre dane pozorne:

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)

Dopasuj model:

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

Wygeneruj niektóre dane, aby przewidzieć i przewidzieć:

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

Narysuj linie regresji:

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)

Możesz zakodować to w prostej funkcji - potrzebujesz tylko kroków z dwóch poprzednich fragmentów kodu - których możesz użyć zamiast 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)
}

Następnie:

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)

Poprzez segmentowany pakiet

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

Na podstawie tych danych nie szacuje punktu przerwania mean(sqft), ale metody ploti linesw tym pakiecie mogą pomóc Ci zaimplementować coś bardziej ogólnego niż myablinewykonanie tego zadania bezpośrednio z dopasowanego lm()modelu.

Edycja: Jeśli chcesz, aby segmentacja oszacowała lokalizację punktu przerwania, ustaw 'psi'argument na NA:

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

Następnie segmentedwypróbuje K = 10kwantyle sqft, z Kustawieniem seg.control()i domyślną wartością 10. Zobacz ?seg.controlwięcej.

Gavin Simpson
źródło
@Gavin (+1) Zdecydowanie bardziej kompletna odpowiedź niż moja; Po prostu to lubię.
chl
@Gavin Sekcja „Przez pakiet podzielony na segmenty” nie działała dla moich danych. Po uruchomieniu segmentedpolecenia otrzymałem komunikat „Nie oszacowano punktu przerwania” .
George Dontas
@ gd047: Przepraszamy, w kodzie pojawił się błąd. Musisz podać argument seq.Zz jednostronną formułą zmiennej (zmiennych), które mają segmentowaną relację z odpowiedzią. Zmodyfikowałem swoją odpowiedź, aby dołączyć seq.Z = ~ sqfti dodałem notatkę o tym, czy segmentedwybrać psidla ciebie wartości .
Gavin Simpson
@ gd047 Chciałbym usunąć moją odpowiedź, ponieważ ta odpowiada na twoje pierwotne pytanie w lepszy sposób. Czy mógłbyś zaakceptować ten zamiast mojego?
chl
@chl Oczywiście, mimo że nadal pojawia się błąd: Błąd w if (model) objF model <- mf: warunek ma długość> 1 i zostanie użyty tylko pierwszy elementmoremil<-mfa:zarsolummintjasnotjantmirprmitzablmizaslosoljadozaljanzarerejatjaon:W.zarnjansolmmisszasolmi:janjafa(moremil)objotfa
George Dontas