I have x, y data where x is position (along a transect) and y is a continuous variable (e.g. temperature).
I want to fit a linear model, a 4 parameter logistic model and a piecewise constant function (aka step-function) with one breakpoint to these data and compare models to determine whether y changes gradually, more abruptly than linear or abruptly in a step-like fashion.
(NOTE: plotting the data reveals that the logistic model is most appropriate in this case but I will ultimately have to loop across many different dependent variables without this information ahead of time).
I am using nls() with manual formula (y~int+slope) for the linear model and SSfpl for the 4-parameter logistic.
l<-nls(y~int+x*slope,data=data.frame(cbind(x,y)),start=list(int=0,slope=1)) # linear model fl<-nls(y~SSfpl(x,A,B,xmid,scal),data=data.frame(cbind(x,y))) # four parameter logistic
My problem is fitting a piecewise constant model. There seem to existing functions for a piecewise linear model, but not many for a step-function. I found:
which gives me the index of an optimal breakpoint in x:
But I'm not sure how to compare the output otherwise to the linear or logistic models. For instance, I think I'll have to do model comparison using AIC but not sure how to calculate this based on the class of the object returned by jumpoints(). So I have manually tried to generate the model using the optimal breakpoint:
Now I can get AIC scores for all the models. But I've noted that when I plot the fit of my manual step-function (pcm) to that generated in cumSeg (pc), it looks like the fit of my manual model isn't as good...
par(mfrow=c(1,2)) plot(x,y) lines(x,predict(cpm) # versus plot(pc)
Wondering why this is happening?
And also looking for alternative suggestions for addressing the overall problem. On the model fitting side of things: the main problem is how do I tell whether or not a variable is stepped versus otherwise changing more or less linearly? (I can have a logistic model that fits better than a linear model but has similar slopes at xmid...so that's not really an answer).
On the model comparison end of things: I'm using AICc() function from the MuMIn package to get the AIC values for the different models and the akaike.weights() function from the qpcR package to do the delta AIC "test"....but maybe there are other ideas?
Thanks! Much appreciated!