[R-sig-ME] PREDICTIONS from a PIECEWISE LINEAR (mixed) MODEL: THEY AIN'T LINEAR BETWEEN BREAK POINTS!!

[Ben Bolker]

On 11-03-19 07:11 AM, Federico Bonofiglio wrote:
> By the way, ..it was to hard to resist...
> I did play a little around....
> I f it doesn't bother u I took the permission to post u the fakies now.
> 
> it really seems that something weird happens to the piecewise model when
> I introduce an interaciton...
> 
> though from statistical theory this should be allowed..
> I don't know if I simply implement a wrong syntax in R... however it
> sounds quite coherent to me......
> 

  Nothing too surprising happened to me (see my version below).  I did
have some convergence issues which I managed to get around (perhaps) by
switching the optimization method (although it may just be that the
convergence criteria are less fussy for the alternative method).  The
predicted values certainly have more noise in them when they incorporate
the effects of other variables, but that's not a surprise.



set.seed(1001)
x<-rnorm(1000)
y<-exp(-x)+rnorm(1000)
plot(x,y)
abline(v=-1,col=2,lty=2)
d <- data.frame(x,y)
mod1 <- lm(y~x+x*(x>-1),data=d)
mod1A <- lm(y~x*(x>-1),data=d)  ## identical
mod1B <- lm(y~x+x:(x>-1),data=d)  ## identical
mod1C <- lm(y~I(x-1):I(x>-1),data=d)  ## if we wanted to force
continuity at I=1

newx <- seq(-3,3,length=101)
yy<-predict(mod,newdata=data.frame(x=newx))
lines(newx,yy,col=2)


#--lme

#grouping factor, unbalanced

g<-as.character(c(1:200))
## this seems a little bit odd.  Why sample from
##  half the groups at random?  Why not sample from
##  half as many groups to start with?
d$id<-sample(g,size=1000,replace=TRUE,
           prob=sample(0:1,200,rep=TRUE))


library(nlme)
mod2<-lme(y~x+x*(x>-1),random=~x|id,
          data=d,
          control=lmeControl(maxIter=1000,msMaxIter=1000,
            opt="optim"))
## convergence problem: this way doesn't complain (although?
##  may not really be better)

summary(mod2)

newframe<-data.frame(x=newx,
                     id="fake")


#predictions

yy2<-predict(mod2,level=0, newdata=newframe)

lines(newx,yy2,col="blue",lwd=2)

## looks fine to me?

# add variable in the model

d$z <- rgamma(1000,4,6)

mod3<-lme(y~x+x*(x>-1)+z,
          random=~x|id,
          data=d,
          control=lmeControl(maxIter=1000,msMaxIter=1000,
            opt="optim"))


summary(mod3)


#new id

newframe2<-data.frame(x=newx,  #fictious id
                      id="fake",
                      z=0)


#predict

yy3<-predict(mod3,level=0, newdata=newframe2)
lines(newx,yy3,col="green",lwd=2)

## still the same

# ADD INTERACTION  z:x

mod4<-lme(y~x+x*(x>-1)+
          z+z:x+
          z:x*(x>-1),
          random=~x|id,
          data=d,
          control=lmeControl(maxIter=1000,msMaxIter=1000,
            opt="optim"))


#predict

yy4<-predict(mod4,level=0, newdata=newframe2)

lines(newx,yy4,col="violet",lwd=2)


yy4B<-predict(mod4,level=0)

lines(sort(d$x),yy4B[order(d$x)],col="cyan",lwd=2)





> I excuse for the annoying if I did cause any
> 
> 
> Big thanks
> 
> Freddy
> 
> 
> 
> 
> 
> x<-rnorm(1000)
> 
> y<-exp(-x)+rnorm(1000)
> 
> plot(x,y)
> abline(v=-1,col=2,lty=2)
> 
> 
> mod<-lm(y~x+x*(x>-1))
> 
> summary(mod)
> 
> yy<-predict(mod)
> 
> lines(x[order(x)],yy[order(x)],col=2,lwd=2)
> 
> 
> #--lme
> 
> #grouping factor, unbalanced
> 
> g<-as.character(c(1:200))
> id<-sample(g,size=1000,replace=T,
> prob=sample(0:1,200,rep=T)) 
> 
> table(id)   #unbalanced
> 
> 
> 
> mod2<-lme(y~x+x*(x>-1),random=~x|id,
> data=data.frame(x,y,id))
> 
> summary(mod2)
> 
> 
> newframe<-data.frame(  #fictious id
> id="fictious",
> x)
> 
> newframe[1:5,]
> 
> #predictions
> 
> yy2<-predict(mod2,level=0, newdata=newframe)
> 
> 
> lines(x[order(x)],yy2[order(x)],col="blue",lwd=2)
> 
> 
> 
> # add variable in the model
> 
> z<-rgamma(1000,4,6)
> 
> mod3<-lme(y~x+x*(x>-1)+z
> ,random=~x|id,
> data=data.frame(x,y,z,id))
> 
> summary(mod3)
> 
> 
> #new id
> 
> newframe2<-data.frame(  #fictious id
> id="fictious",
> x,
> z)
> 
> 
> #predict
> 
> yy3<-predict(mod3,level=0, newdata=newframe2)
> 
> 
> lines(x[order(x)],yy3[order(x)],col="green",lwd=2)
> 
> 
> 
> # ADD INTERACTION  z:x
> 
> 
> 
> mod4<-lme(y~x+x*(x>-1)+
> 
> z+
>    z:x+
>          z:x*(x>-1)
> 
>    ,random=~x|id,
> data=data.frame(x,y,z,id))
> 
> 
> 
> #predict
> 
> yy4<-predict(mod4,level=0, newdata=newframe2)
> 
> 
> lines(x[order(x)],yy4[order(x)],col="violet",lwd=2)  #something bizarre
>                                                     #starts to happen
>                                                     #in the predicted values
> 
> # they begin to jiggle around the straight line
> 
> 
> 
> 
>