Dear all,
I'm trying to understand how the strucchange package is working and I have been looking at the examples given for the Fstats() function.
The first example (Nile), shows one peak in the F-stats and one breakpoint is estimated, that can be plotted using the following code

## Nile data with one breakpoint: the annual flows drop in 1898
## because the first Ashwan dam was built
data("Nile")
plot(Nile)
## test the null hypothesis that the annual flow remains constant
## over the years
fs.nile <- Fstats(Nile ~ 1)
plot(fs.nile, alpha=0.05)
sctest(fs.nile)
## visualize the breakpoint implied by the argmax of the F statistics
plot(Nile)
lines(breakpoints(fs.nile))


In the second example, the authors state the presence of "at least" two breakpoints. When plotting the F-statistics using the following code, we see indeed two peaks in the F-statistics, that coincides with the dates given by the authors: c.a 1973 and 1983 but when trying to add those breakpoints to the time series, only one is taken into account


## UK Seatbelt data: a SARIMA(1,0,0)(1,0,0)_12 model
## (fitted by OLS) is used and reveals (at least) two
## breakpoints - one in 1973 associated with the oil crisis and
## one in 1983 due to the introduction of compulsory
## wearing of seatbelts in the UK.
data("UKDriverDeaths")
seatbelt <- log10(UKDriverDeaths)
seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12))
colnames(seatbelt) <- c("y", "ylag1", "ylag12")
seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12))
plot(seatbelt[,"y"], ylab = expression(log[10](casualties)))
## compute F statistics for potential breakpoints between
## 1971(6) (corresponds to from = 0.1) and 1983(6) (corresponds to
## to = 0.9 = 1 - from, the default)
## compute F statistics
fs <- Fstats(y ~ ylag1 + ylag12, data = seatbelt, from = 0.1)
## this gives the same result
fs <- Fstats(y ~ ylag1 + ylag12, data = seatbelt, from = c(1971, 6),
to = c(1983, 6))
## plot the F statistics
plot(fs, alpha = 0.05)

###TRYING to plot the 2 breakpoints:
lines(breakpoints(fs,breaks=2))

...only one breakpoint is plotted. When looking at the results from

> breakpoints(fs,breaks=2)

Optimal 2-segment partition:

Call: breakpoints.Fstats(obj = fs, breaks = 2)
Breakpoints at observation number:
46
Corresponding to breakdates:
1973(10)

We see that even though the F-statistics seem to show the existence of 2 breakpoints, only one is detected by the breakpoints() function. Does anyone know how this is possible? I'm totally new to strucchange so it might well be something obvious I'm missing here!


OTHER SIDE QUESTION: can strucchange be used if the y variable is binary???

Thanks heaps if anyone can help
Geraldine

From: Mabille, Geraldine
Sent: 29. mai 2012 11:30
To: 'r-help@r-project.org'
Subject: GAM interactions, by example

Dear all,
I'm using the mgcv library by Simon Wood to fit gam models with interactions and I have been reading (and running) the "factor 'by' variable example"   given on the gam.models help page (see below, output from the two first models b, and b1).
The example explains that both b and b1 fits are similar: "note that the preceding fit (here b) is the same as (b1)...."
I agree with the idea that it "looks" the same but when I look at the results from both models (summary b and summary b1) I see that the results look in fact quite different (edf, and also deviance explained for example???)
Are those two models (b and b1) really testing the same things??? If yes, why are the results so different between models???
Thanks a lot if anyone can help with that...
Geraldine


dat <- gamSim(4)

## fit model...
b <- gam(y ~ fac+s(x2,by=fac)+s(x0),data=dat)
plot(b,pages=1)
summary(b)

Family: gaussian
Link function: identity

Formula:
y ~ fac + s(x2, by = fac) + s(x0)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.1784     0.1985   5.937 6.59e-09 ***
fac2         -1.2148     0.2807  -4.329 1.92e-05 ***
fac3          2.2012     0.2436   9.034  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
             edf Ref.df      F  p-value
s(x2):fac1 5.364  6.472  2.285   0.0312 *
s(x2):fac2 4.523  5.547 11.396 4.59e-11 ***
s(x2):fac3 8.024  8.741 43.456  < 2e-16 ***
s(x0)      1.000  1.000  0.237   0.6269
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.634   Deviance explained = 65.3%
GCV score = 4.0288  Scale est. = 3.8082    n = 400

## note that the preceding fit is the same as....
b1<-gam(y ~ s(x2,by=as.numeric(fac==1))+s(x2,by=as.numeric(fac==2))+
            s(x2,by=as.numeric(fac==3))+s(x0)-1,data=dat)
## ... the `-1' is because the intercept is confounded with the
## *uncentred* smooths here.
plot(b1,pages=1)
summary(b1)

Family: gaussian
Link function: identity

Formula:
y ~ s(x2, by = as.numeric(fac == 1)) + s(x2, by = as.numeric(fac ==
    2)) + s(x2, by = as.numeric(fac == 3)) + s(x0) - 1

Approximate significance of smooth terms:
                             edf Ref.df       F  p-value
s(x2):as.numeric(fac == 1) 6.341  7.447   6.214 3.38e-07 ***
s(x2):as.numeric(fac == 2) 3.393  3.961  14.727 4.07e-11 ***
s(x2):as.numeric(fac == 3) 9.015  9.737 104.760  < 2e-16 ***
s(x0)                      1.000  1.000   0.266    0.606
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.631   Deviance explained =   75%
GCV score = 4.0345  Scale est. = 3.8353    n = 400


	[[alternative HTML version deleted]]

