Model liniowy - prognozowanie

Wersja pdf

Załadujmy ten sam model co poprzednio:

library(tidyverse)
devtools::install_github("kassambara/datarium")
data("marketing", package = "datarium")
model <- lm(sales ~ youtube + facebook + newspaper, data = marketing)

Stwórzmy nowe dane testowe:

new<-data.frame( youtube = c(12, 19, 24), facebook=c(40,50,60), newspaper=c(25,55,85) )
new
##   youtube facebook newspaper
## 1      12       40        25
## 2      19       50        55
## 3      24       60        85

Wykonajmy prognozowanie:

pred<-predict(model, newdata = new)
pred
##        1        2        3 
## 11.59111 13.76563 15.84863

lub ręcznie:

x.new<-cbind(rep(1,3),as.matrix(new))
pred2<-x.new %*% model$coefficients
pred2
##          [,1]
## [1,] 11.59111
## [2,] 13.76563
## [3,] 15.84863

Wersja pełna:

predf<-predict(model, new, se.fit = TRUE)
predf
## $fit
##        1        2        3 
## 11.59111 13.76563 15.84863 
## 
## $se.fit
##         1         2         3 
## 0.3068124 0.3271400 0.4201326 
## 
## $df
## [1] 196
## 
## $residual.scale
## [1] 2.022612

se.fit to odchylenie standardowe względem średniejfit.

Przedziały ufności

Przedziały ufności możemy otrzymać następująco:

predict(model, new, se.fit = TRUE, interval = "confidence")
## $fit
##        fit      lwr      upr
## 1 11.59111 10.98603 12.19618
## 2 13.76563 13.12047 14.41080
## 3 15.84863 15.02007 16.67719
## 
## $se.fit
##         1         2         3 
## 0.3068124 0.3271400 0.4201326 
## 
## $df
## [1] 196
## 
## $residual.scale
## [1] 2.022612

lub ręcznie za pomocą odpowiednich kwantyli rozkładu \(t\):

kt <- c(-1, 1) * qt(0.05 / 2, predf$df, lower.tail = FALSE)
kt
## [1] -1.972141  1.972141
pu <- predf$fit + outer(predf$se.fit, kt)
pu
##       [,1]     [,2]
## 1 10.98603 12.19618
## 2 13.12047 14.41080
## 3 15.02007 16.67719

outer - tutaj mnożenie wyraz za wyrazem.

Innymi słowymi, z prawodpodobieństwem 95% średnie wyniki sprzedaży są zwarte w odpowiednim przedziale.

Wiarygodność przewidywań

Przedziały ufności możemy otrzymać następująco:

predict(model, new, se.fit = TRUE, interval = "prediction")
## $fit
##        fit       lwr      upr
## 1 11.59111  7.556598 15.62562
## 2 13.76563  9.724919 17.80635
## 3 15.84863 11.774611 19.92265
## 
## $se.fit
##         1         2         3 
## 0.3068124 0.3271400 0.4201326 
## 
## $df
## [1] 196
## 
## $residual.scale
## [1] 2.022612

lub ręcznie:

se.PI <- sqrt(predf$se.fit ^ 2 + predf$residual.scale ^ 2)
wp <- predf$fit + outer(se.PI, kt)
wp
##        [,1]     [,2]
## 1  7.556598 15.62562
## 2  9.724919 17.80635
## 3 11.774611 19.92265

Interpretacja graficzna

Wygenerujmy więcej nowych danych:

n <- 20    
yt <- runif(n, min = 15, max = 150)+rnorm(4,2,0.5)
fb <- runif(n, min = 30, max = 70)
np <- runif(n, min=50, max=144)
new2<-data.frame( youtube = yt, facebook=fb, newspaper=np )
pred2p<-predict(model, new2, interval = "prediction")
head(pred2p)
##        fit       lwr      upr
## 1 18.27201 14.114844 22.42918
## 2 13.92545  9.888016 17.96288
## 3 17.88640 13.818078 21.95472
## 4 15.30945 11.240940 19.37796
## 5 18.39750 14.367768 22.42723
## 6 11.74739  7.607411 15.88737
pred2c<-predict(model, new2, interval = "confidence")
head(pred2c)
##        fit      lwr      upr
## 1 18.27201 17.10116 19.44287
## 2 13.92545 13.30118 14.54972
## 3 17.88640 17.08633 18.68647
## 4 15.30945 14.50844 16.11046
## 5 18.39750 17.82514 18.96986
## 6 11.74739 10.63911 12.85566
plot(pred2c[,1], ylab="", ylim=c(min(pred2p[,2]), max(pred2p[,3])), pch=20)
lines(pred2c[,2], col="blue", lty=2)
lines(pred2p[,2], col="red", lty=3)
lines(pred2c[,3], col="blue", lty=2)
lines(pred2p[,3], col="red", lty=3)