vignettes/jointfp_model.R

In mfp: Multivariable Fractional Polynomials

library(mfp)


x <- xorig <- runif(1000, min=-0.5, max=1.5)
f <- rbinom(1000, 1, 0.5)*1
xf <- x
xf[f==1] <- 0

y <- 5 + 3*x^2+rnorm(1000, 0, 1)

fit <- mfp(y~fp(xf)+fp(x)+f)




p<-predict(fit, type="response", newdata=data.frame(x=x))
p
pt <- predict(fit, type="terms", terms="x", seq=list(seq(-0.4, 0.4, 0.01)), ref=list(0))
pt


y2 <- 5 + (x>0)*3*x^2 + (x<0)*10*x^3 + (x<0)*3*x+rnorm(1000, 0, 1)
dataset <- data.frame(x=x, y=y, y2=y2)

plot(x, y)
plot(x, y2)

# simple mfp
fit <- mfp(y~fp(x), dataset)
summary(fit)

p<-predict(fit, type="terms", terms="x", seq=list(seq(-0.5, 1.5, 0.01)), ref=list(0))

plot(p$x$variable, p$x$contrast, type="l")


# joint mfp

x1 <- x2 <- xorig
x1[xorig>0.5] <- 0.5
x2[xorig<0.5] <- 0.5
dataset <- data.frame(x=xorig, y=y, x1=x1, x2=x2)

fitj <- mfp(y~fp(x1)+ fp(x2), dataset)
summary(fitj)
x<- seq(-0.5, 1.5, 0.01)
x2 <- x1 <- x
x1[x>0.5] <- 0.5
x2[x<0.5] <- 0.5

#lp <- 3.6209 -3.6979*(x1+1) + 3.0029*(x1+1)*log(x1+1) + 2.3824*(x2+0.1)^2
#lp <- 0.205 + 2.65663*(x2+0.1)^2 -0.18889*(x1+0.2)/0.1
#lp <- 5.36799 + 2.69344*(x2+0.1)^2 - 1.12503*(x1+0.5)
lp <- 5.09896 + 3.04023*x2^2 - 3.19599*(x1+0.5) + 3.07416 * (x1+0.5)^2

plot(xorig[1:100], y[1:100])
lines(x, lp)
abline(h=0, lty=2)



### y2
# simple mfp
fit2 <- mfp(y2~fp(x), dataset)
summary(fit2)

p<-predict(fit2, type="terms", terms="x", seq=list(seq(-0.49, 1.5, 0.01)), ref=list(0))

plot(p$x$variable, p$x$contrast, type="l")


# joint mfp

x1 <- x2 <- xorig
x2p <- function(x, joinpoint=0, lower=TRUE){
  y <- x
  if(lower) y[x>joinpoint] <- joinpoint
  else y[x<joinpoint] <- joinpoint
  return(y)
}
x1[xorig>0] <- 0
x2[xorig<0] <- 0
dataset <- data.frame(x=xorig, y=y, x1=x1, x2=x2, y2=y2)

fitj2 <- mfp(y2~fp(x1)+ fp(x2), dataset)
fitj3 <- mfp(y2~fp(x2p(x, joinpoint=0, lower=TRUE))+ fp(x2p(x, joinpoint=0, lower=FALSE)), dataset)
summary(fitj2)
summary(fitj3)


x<- seq(-0.49, 1.5, 0.01)
x2 <- x1 <- x
x1[x>0] <- 0
x2[x<0] <- 0

#lp <- 3.6209 -3.6979*(x1+1) + 3.0029*(x1+1)*log(x1+1) + 2.3824*(x2+0.1)^2
#lp <- 0.205 + 2.65663*(x2+0.1)^2 -0.18889*(x1+0.2)/0.1
#lp <- 5.36799 + 2.69344*(x2+0.1)^2 - 1.12503*(x1+0.5)
#lp <- 5.09896 + 3.04023*x2^2 - 3.19599*(x1+0.5) + 3.07416 * (x1+0.5)^2
#lp <- 5.02938 + 1.29664*log(x1+0.5) + 0.10722*(log(x1+0.5))^2 + 2.94951*x2^2
lp <- 2.25368 + 8.96298*(x1+0.5) - 14.71014*(x1+0.5)^3 + 2.60244*(x2+0.1)^2

plot(xorig[1:100], y2[1:100])
lines(x, lp)
abline(h=0, lty=2)

p <- predict(fitj3, newdata=dataset, type="response")
p <- predict(fitj3, type="terms", terms="x", seq=list(seq(-0.49, 0.5, 0.01)), ref=list(0))