Search code examples
rgraphstability

How to plot Chow-test results in R?

I have run a Chow-test in R using the following code:

mydata <- read.csv(file="chow.csv", header=TRUE, sep=",")
sctest(fuel~pred, data=mydata, type="Chow", point=44)

I am able to successfully obtain the F Statistic and p-value, but I would like to see the graph, i.e. the splitting regression where you have a curve divided into two (at point 44), like the image attached.

Could anybody tell me if R can plot this graph with the information given? Thank you. enter image description here

Solution

  • I don't think that you can plot this with the results from sctest but it should be easy to plot it from your data. Simply create each of the relevant linear regressions and then plot the line segments.

    Since you do not provide any data, I will illustrate with the NHANES data from the NHANES package.

    library(NHANES)
data(NHANES)
plot(Height ~ Age, data=NHANES, pch=20, col="#00000022")

    Underlying data

    There seems to be a structural change in Height as a function of Age at Age=16. You want to compare the linear model on the full data with the two linear models, Age<=16 and Age > 16. So just generate them and plot them.

    ## Full data
LM0 = lm(Height ~ Age, data=NHANES)
abline(LM0, col="red", type=2, lwd=3)

## Age <= 16
LM1 = lm(Height ~ Age, data=NHANES[NHANES$Age <= 16,])
y1 = predict(LM1, newdata=data.frame(Age=c(0,16)))

## Age > 16
LM2 = lm(Height ~ Age, data=NHANES[NHANES$Age > 16,])
y2 = predict(LM2, newdata=data.frame(Age=c(16,80)))

lines(x=c(0,16), y=y1, col="green", lty=3, lwd=3)
lines(x=c(16,80), y=y2, col="green", lty=3, lwd=3)

    Plot with regression lines

    Edit

    Based on data provided in the comments, here is another example. At least with the few points given, dividing this data into two models does not seem obvious, but I picked a division point that emphasized different models.

    mydata =  structure(list(fuel = c(346591L, 344841L, 369015L, 346062L, 
    337494L, 355830L, 380516L, 391917L, 356761L, 359326L, 348208L, 
    358639L, 351659L, 323356L, 364712L), pred = c(346293.4, 344855.5, 
    355014.2, 342456.1, 347042.1, 371000.1, 376011.6, 391521.4, 
    364751.5, 358650, 356547.3, 355159.8, 359598.2, 349122.6, 357969.6)), 
    .Names = c("fuel", "pred"), row.names = c(NA, 15L), class = "data.frame")

plot(pred ~ fuel, data=mydata)
LM0 = lm(pred ~ fuel, data=mydata)
abline(LM0, lty=2)

LM1 = lm(pred ~ fuel, data=mydata[mydata$fuel <= 347500,])
y1  = predict(LM1, newdata=data.frame(fuel=c(320000,347500)))
lines(c(320000,350000), y1, lty=3)

LM2 = lm(pred ~ fuel, data=mydata[mydata$fuel >  347500,])
y2  = predict(LM2, newdata=data.frame(fuel=c(347500,395000)))
lines(c(350000,395000), y2, lty=3)

    Second example.