Предсказать сегментированный ЛМ за пределами пакета - PullRequest
Новая
       64

Предсказать сегментированный ЛМ за пределами пакета

0 голосов
/ 08 июня 2018

У меня есть массив выходных данных от сотен сегментированных линейных моделей (созданных с использованием сегментированного пакета в R).Я хочу иметь возможность использовать эти выходные данные для новых данных, используя функцию прогнозирования.Чтобы было ясно, у меня нет сегментированных объектов линейной модели в моей рабочей области;Я просто сохранил и заново импортировал соответствующие результаты (например, коэффициенты и точки останова).По этой причине я не могу просто использовать функциюgnast.segmented из сегментированного пакета.

Ниже приведен пример игрушки, основанный на этой ссылке , которая кажется многообещающей, но не соответствует выходным данным.функции Предсказание. Сегментированная. 

library(segmented) 
set.seed(12) 

xx <- 1:100 
zz <- runif(100) 
yy <- 2 + 1.5*pmax(xx-35,0) - 1.5*pmax(xx-70,0) + 
               15*pmax(zz-0.5,0) + rnorm(100,0,2) 
dati <- data.frame(x=xx,y=yy,z=zz) 

out.lm<-lm(y~x,data=dati) 
o<-## S3 method for class 'lm': 
     segmented(out.lm,seg.Z=~x,psi=list(x=c(30,60)), 
           control=seg.control(display=FALSE)) 

# Note that coefficients with U in the name are differences in slopes, not slopes. 
# Compare:
slope(o)
coef(o)[2] + coef(o)[3]
coef(o)[2] + coef(o)[3] + coef(o)[4]

# prediction 
pred <- data.frame(x = 1:100) 
pred$dummy1 <- pmax(pred$x - o$psi[1,2], 0) 
pred$dummy2 <- pmax(pred$x - o$psi[2,2], 0) 
pred$dummy3 <- I(pred$x > o$psi[1,2]) * (coef(o)[2] + coef(o)[3])
pred$dummy4 <- I(pred$x > o$psi[2,2]) * (coef(o)[2] + coef(o)[3] + coef(o)[4]) 
names(pred)[-1]<- names(model.frame(o))[-c(1,2)] 

# compute the prediction, using standard predict function 
# computing confidence intervals further 
# suppose that the breakpoints are fixed 
pred <- data.frame(pred, predict(o, newdata= pred, 
                       interval="confidence")) 

# Try prediction using the predict.segment version to compare
test <- predict.segmented(o)
plot(pred$fit, test, ylim = c(0, 100))
abline(0,1, col = "red")
# At least one segment not being predicted correctly?

Могу ли я использовать базовую функцию r предиката () (не функцию segmented.predict ()) с коэффициентами и точками останова, сохраненными из сегментированных линейных моделей?

ОБНОВЛЕНИЕ Я понял, что в коде выше есть проблемы (не используйте его).С помощью некоторого обратного инжиниринга функции segmented.predict () я создал матрицу проектирования и использую ее для прогнозирования значений, а не непосредственно с помощью функции предиката ().Я пока не считаю это полным ответом на первоначальный вопрос, потому что предикат () также может дать доверительные интервалы для предсказания, и я еще не реализовал это - вопрос еще открыт для кого-то, чтобы добавить доверительные интервалы. 

library(segmented)


## Define function for making matrix of dummy variables (this is based on code from predict.segmented())
dummy.matrix <- function(x.values, x_names, psi.est = TRUE, nameU, nameV, diffSlope, est.psi) {
  # This function creates a model matrix with dummy variables for a segmented lm with two breakpoints.
  # Inputs:
  # x.values: the x values of the segmented lm
  # x_names: the name of the column of x values
  # psi.est: this is legacy from the predict.segmented function, leave it set to 'TRUE'
  # obj: the segmented lm object
  # nameU: names (class character) of 3rd and 4th coef, which are "U1.x" "U2.x" for lm with two breaks. Example: names(c(obj$coef[3], obj$coef[4]))
  # nameV: names (class character) of 5th and 6th coef, which are "psi1.x" "psi2.x" for lm with two breaks. Example: names(c(obj$coef[5], obj$coef[6]))
  # diffSlope: the coefficients (class numeric) with the slope differences; called U1.x and U2.x for lm with two breaks. Example: c(o$coef[3], o$coef[4])
  # est.psi: the estimated break points (class numeric); these are the estimated breakpoints from segmented.lm. Example: c(obj$psi[1,2], obj$psi[2,2])
  #
  n <- length(x.values)
  k <- length(est.psi)
  PSI <- matrix(rep(est.psi, rep(n, k)), ncol = k)
  newZ <- matrix(x.values, nrow = n, ncol = k, byrow = FALSE)
  dummy1 <- pmax(newZ - PSI, 0)
  if (psi.est) {
    V <- ifelse(newZ > PSI, -1, 0)
    dummy2 <- if (k == 1) 
      V * diffSlope
    else V %*% diag(diffSlope)
    newd <- cbind(x.values, dummy1, dummy2)
    colnames(newd) <- c(x_names, nameU, nameV)
  } else {
    newd <- cbind(x.values, dummy1)
    colnames(newd) <- c(x_names, nameU)
  }
  # if (!x_names %in% names(coef(obj.seg))) 
  #   newd <- newd[, -1, drop = FALSE]
  return(newd)
}

## Test dummy matrix function----------------------------------------------
set.seed(12)
xx<-1:100
zz<-runif(100)
yy<-2+1.5*pmax(xx-35,0)-1.5*pmax(xx-70,0)+15*pmax(zz-.5,0)+rnorm(100,0,2)
dati<-data.frame(x=xx,y=yy,z=zz)
out.lm<-lm(y~x,data=dati)

#1 segmented variable, 2 breakpoints: you have to specify starting values (vector) for psi:
o<-segmented(out.lm,seg.Z=~x,psi=c(30,60),
             control=seg.control(display=FALSE))
slope(o)
plot.segmented(o)
summary(o)


# Test dummy matrix fn with the same dataset
newdata <- dati
nameU1 <- c("U1.x", "U2.x")
nameV1 <- c("psi1.x", "psi2.x")
diffSlope1 <- c(o$coef[3], o$coef[4])
est.psi1 <- c(o$psi[1,2], o$psi[2,2])

test <- dummy.matrix(x.values = newdata$x, x_names = "x", psi.est = TRUE, 
                     nameU = nameU1, nameV = nameV1, diffSlope = diffSlope1, est.psi = est.psi1)


# Predict response variable using matrix multiplication
col1 <- matrix(1, nrow = dim(test)[1])
test <- cbind(col1, test) # Now test is the same as model.matrix(o)
predY <- coef(o) %*% t(test)
plot(predY[1,])
lines(predict.segmented(o), col = "blue") # good, predict.segmented gives same answer
...