# Tag Archives: smoothing

On Thursday, March 2nd, I will give the first lecture of the PhD course on advanced tools for econometrics, on nonlinearities. Slides are available online.

I will give a short graduate course for PhD students, in Rennes, on Thurday mornings, in March (2nd, 9th, 23rd and 30th). The agenda will be

1. Nonlinear Regression Models and Smoothing Techniques

2. Bootstrapping and Regression

3. Penalized Regression Models and LASSO

4. Quantile Regression and Expectiles

There will be slides available by the end of February.

# Statistics, and the Goldilocks Principle

By the end of May, in Toronto, we had that great talk at the SSC by Jeff Rosenthal, on monte carlo techniques, and Jeff mention the name of “the Goldilocks principle” (it was in the contect of MCMC, and I did mention it in my talk in London on MCMC, when I discussed the value of the rejection rate of the Hastings Metropolis algorithm, which should be not to large, and not too small…). In the story, Goldilocks, there are always three alternative, one is always too much in one extreme (too hot – for the soup – or too large – for the bed, or the chaiir), one is too much in the opposite extreme (too cold, or too small), and one is “just right“.

# Conditional dependence measures

This week, I spend some time at the Workshop on Nonparametric Curve Smoothing conference at Concordia. Yesterday afternoon, Noël Veraverbeke show an interesting graph, to illustrate conditional copulas (and the derivation of conditional dependence measures, such as Kendall’s tau, or Spearman’s rho). A long time ago, in my PhD thesis (mainly on conditional copulas) I did try to derive conditional dependence measures (in a dedicated chapter). In my PhD, I was interested to describe the dependence of a pair $(Y_1,Y_2)$ given $(Y_1,Y_2)\in\mathcal{V}$, where $\mathcal V$ is a region of interest, such has tails. So I wanted to study the behavior of $(Y_1,Y_2)$ given $\{Y_1>t,Y_2>t\}$. This has interpretation when studying large risks, but also in joint life mortality.

In the paper Noël mentioned, they want to describe the dependence of a pair $(Y_1,Y_2)$ given a covariate $X=x$. And he came up with this very nice example: consider expected lifetimes, for male and female, in various countries. You can get zipped files with data for male, female and we can use the GPD per capita as our covariate. Here is the code to visualize life expectancies,

b1=read.table("sp.dyn.le00.fe.in_Indicator_en_csv_v2.csv",header=TRUE,sep=",",skip=2)
b1b=b1[,c(1,2,55)]
b2b=b2[,c(1,2,55)]
b3b=b3[,c(1,2,55)]
names(b1b)[3]="LEF"
names(b2b)[3]="LEM"
names(b3b)[3]="GPD"
b=merge(b1b,b2b)
b=merge(b,b3b)
plot(b$LEM,b$LEF,xlab="Life Expectancy (male vs. female)")

With this graph, we cannot visualize the link with the covariate,

b$cgpd=cut(b$GPD,quantile(b$GPD,seq(0,1,by=1/6),na.rm=TRUE)) levels(b$cgpd)=as.character(1:6)
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")
plot(b$LEM,b$LEF,xlab="Life Expectancy (male vs. female)",pch=19,col=CL[as.numeric(b$cgpd)]) Here, poor countries are in red, and rich countries in blue, Clearly, life expectancy is connected to the wealth of the country, plot(b$GPD,b$LEF,xlab="(Female) Life Expectancy vs. GPD (log scale)",pch=19,col=CL[as.numeric(b$cgpd)],log="x")
plot(b$GPD,b$LEM,xlab="(Male) Life Expectancy vs. GPD (log scale)",pch=19,col=CL[as.numeric(b$cgpd)],log="x") The idea here is to consider the conditional dependence structure, given the wealth. If we want something smooth (this is actually the goal of the workshop, but I’d like to make that quickly) consider some weighted version of Kendall’s tau, based on the idea mentioned in a post on http://stackoverflow.com/ The idea is to use concordance and discordance counts, with replications of the data, based on the weights P = function(t) { r_ndx = row(t) c_ndx = col(t) sum(t * mapply(function(r, c){sum(t[(r_ndx > r) & (c_ndx > c)])}, r = r_ndx, c = c_ndx))} Q = function(t) { r_ndx = row(t) c_ndx = col(t) sum(t * mapply( function(r, c){ sum(t[(r_ndx > r) & (c_ndx < c)]) }, r = r_ndx, c = c_ndx) ) } kendall_tau_c = function(t){ t = as.matrix(t) m = min(dim(t)) n = sum(t) ks_tauc = (m*2*(P(t)-Q(t)))/((n*n)*(m-1)) } I=is.na(b$GPD)
bw=density(log(b$GPD[!I]))$bw
kendall.weight=function(x){
df=data.frame(Y1=b$LEF, Y2=b$LEM, freq=trunc(dnorm(log(b$GPD)-log(x),sd=bw)*100)) df=df[!is.na(df$freq),]
dfrep=data.frame( lapply(df, function(x){rep(x, df$freq)})) t=xtabs(~ Y1+Y2, dfrep) return(kendall_tau_c(t))} Here, I use weights using some Gaussian kernel on the logarithm of the GPD per capita (my standard deviation for the Gaussian weight being equal to the bandwidth of the Gaussian kernel of the density of the log of the GPD per capita), then, we can compute various conditional Kendall’s tau, T=exp(seq(6,11.5,length=50)) K=Vectorize(kendall.weight)(T) and plot them, plot(T,K,type="l",xlab="Conditional Kendall's tau vs. GPD (log scale)") There is more “correlation” between lifetimes of men and women in poor countries than rich country (which is also what Noël observed). Now, we can also play with time, because we have those statistics for several years. # Smoothing mortality rates This morning, I was working with Julie, a student of mine, coming from Rennes, on mortality tables. Actually, we work on genealogical datasets from a small region in Québec, and we can observe a lot of volatiliy. If I borrow one of her graph, we get something like Since we have some missing data, we wanted to use some Generalized Nonlinear Models. So let us see how to get a smooth estimator of the mortality surface. We will write some code that we can use on our data later on (the dataset we have has been obtained after signing a lot of official documents, and I guess I cannot upload it here, even partially). DEATH <- read.table( "http://freakonometrics.free.fr/Deces-France.txt", header=TRUE) EXPO <- read.table( "http://freakonometrics.free.fr/Exposures-France.txt", header=TRUE,skip=2) library(gnm) D=DEATH$Male
E=EXPO$Male A=as.numeric(as.character(DEATH$Age))
Y=DEATH$Year I=(A<100) base=data.frame(D=D,E=E,Y=Y,A=A) subbase=base[I,] subbase=subbase[!is.na(subbase$A),]

The first idea can be to use a Poisson model, where the mortality rate is a smooth function of the age and the year, something like

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s(x,t)}])$that can be estimated using

library(mgcv)
regbsp=gam(D~s(A,Y,bs="cr")+offset(log(E)),data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regbsp,newdata=data.frame(A=a,Y=y,E=1))
vX=trunc(seq(0,99,length=41))
vY=trunc(seq(1900,2005,length=41))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The mortality surface is here

It is also possible to extract the average value of the years, which is the interpretation of the $a_x$ coefficient in the Lee-Carter model,

predAx=function(a) mean(predict(regbsp,newdata=data.frame(A=a,
Y=seq(min(subbase$Y),max(subbase$Y)),E=1)))
plot(seq(0,99),Vectorize(predAx)(seq(0,99)),col="red",lwd=3,type="l")

We have the following smoothed mortality rate

Recall that the Lee-Carter model is

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}a_x+b_x\cdot k_t}])$

where parameter estimates can be obtained using

regnp=gnm(D~factor(A)+Mult(factor(A),factor(Y))+offset(log(E)),
data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regnp,newdata=data.frame(A=a,Y=y,E=1))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The (crude) mortality surface is

with the following $a_x$ coefficients.

plot(seq(1,99),coefficients(regnp)[2:100],col="red",lwd=3,type="l")

Here we have a lot of coefficients, and unfortunately, on a smaller dataset, we have much more variability. Can we smooth our Lee-Carter model ? To get something which looks like

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s_a(x)+s_b(x)\cdot s_k(t)}])$

Actually, we can, and the code is rather simple

library(splines)
knotsA=c(20,40,60,80)
knotsY=c(1920,1945,1980,2000)
regsp=gnm(D~bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+
Mult(bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3),
bs(subbase$Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))+
offset(log(E)),data=subbase, family=quasipoisson)
BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3) BpY=bs(seq(min(subbase$Y),max(subbase$Y)),knots=knotsY,Boundary.knots= range(subbase$Y),degre=3)
predmodel=function(a,y)
predict(regsp,newdata=data.frame(A=a,Y=y,E=1)) v
Z=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The mortality surface is now

and again, it is possible to extract the average mortality rate, as a function of the age, over the years,

BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3) Ax=BpA%*%coefficients(regsp)[2:8] plot(seq(0,99),Ax,col="red",lwd=3,type="l") We can then play with the smoothing parameters of the spline functions, and see the impact on the mortality surface knotsA=seq(5,95,by=5) knotsY=seq(1910,2000,by=10) regsp=gnm(D~bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+
Mult(bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3), bs(Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))
+offset(log(E)),data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regsp,newdata=data.frame(A=a,Y=y,E=1))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

We now have to use those functions our our small data sample ! That should be fun….

# Some heuristics about spline smoothing

Let us continue our discussion on smoothing techniques in regression. Assume that .$\mathbb{E}(Y\vert X=x)=h(x)$ where $h(\cdot)$ is some unkown function, but assumed to be sufficently smooth. For instance, assume that $h(\cdot)$ is continuous, that $h'(\cdot)$ exists, and is continuous, that  $h''(\cdot)$ exists and is also continuous, etc. If $h(\cdot)$ is smooth enough, Taylor’s expansion can be used. Hence, for $x\in(\alpha,\beta)$

$h(x)=h(\alpha)+\sum_{k=1}^ d \frac{(x-\alpha)^k}{k!}h^{(k)}(x_0)+\frac{1}{d!}\int_{\alpha}^x [x-t]^d h^{(d+1)}(t)dt$

which can also be writen as

$h(x)=\sum_{k=0}^ d a_k (x-\alpha)^k +\frac{1}{d!}\int_{\alpha}^x [x-t]^d h^{(d+1)}(t)dt$

for some $a_k$‘s. The first part is simply a polynomial.

The second part, is some integral. Using Riemann integral, observe that

$\frac{1}{d!}\int_{\alpha}^x [x-t]^d h^{(d+1)}(t)dt\sim \sum_{i=1}^ j b_i (x-x_i)_+^d$

for some $b_i$‘s, and some

$\alpha < x_1< x_2< \cdots < x_{j-1} < x_j < \beta$

Thus,

$h(x) \sim \sum_{k=0}^ d a_k (x-\alpha)^k +\sum_{i=1}^ j b_i (x-x_i)_+^d$

Nice! We have our linear regression model. A natural idea is then to consider a regression of $Y$ on $\boldsymbol{X}$ where

$\boldsymbol{X} = (1,X,X^2,\cdots,X^d,(X-x_1)_+^d,\cdots,(X-x_k)_+^d )$

given some knots $\{x_1,\cdots,x_k\}$. To make things easier to understand, let us work with our previous dataset,

plot(db)

If we consider one knot, and an expansion of order 1,

attach(db)
library(splines)
B=bs(xr,knots=c(3),Boundary.knots=c(0,10),degre=1)
reg=lm(yr~B)
lines(xr[xr<=3],predict(reg)[xr<=3],col="red")
lines(xr[xr>=3],predict(reg)[xr>=3],col="blue")

The prediction obtained with this spline can be compared with regressions on subsets (the doted lines)

reg=lm(yr~xr,subset=xr<=3)
lines(xr[xr<=3],predict(reg)[xr<=3],col="red",lty=2)
reg=lm(yr~xr,subset=xr>=3)
lines(xr[xr>=3],predict(reg),col="blue",lty=2)

It is different, since we have here three parameters (and not four, as for the regressions on the two subsets). One degree of freedom is lost, when asking for a continuous model. Observe that it is possible to write, equivalently

reg=lm(yr~bs(xr,knots=c(3),Boundary.knots=c(0,10),degre=1),data=db)

So, what happened here?

B=bs(xr,knots=c(2,5),Boundary.knots=c(0,10),degre=1)
matplot(xr,B,type="l")
abline(v=c(0,2,5,10),lty=2)

Here, the functions that appear in the regression are the following

Now, if we run the regression on those two components, we get

B=bs(xr,knots=c(2,5),Boundary.knots=c(0,10),degre=1)
matplot(xr,B,type="l")
abline(v=c(0,2,5,10),lty=2)

If we add one knot, we get

the prediction is

reg=lm(yr~B)
lines(xr,predict(reg),col="red")

Of course, we can choose much more knots,

B=bs(xr,knots=1:9,Boundary.knots=c(0,10),degre=1)
reg=lm(yr~B)
lines(xr,predict(reg),col="red")

We can even get a confidence interval

reg=lm(yr~B)
P=predict(reg,interval="confidence")
plot(db,col="white")
polygon(c(xr,rev(xr)),c(P[,2],rev(P[,3])),col="light blue",border=NA)
points(db)
reg=lm(yr~B)
lines(xr,P[,1],col="red")
abline(v=c(0,2,5,10),lty=2)

And if we keep the  two knots we chose previously, but consider Taylor’s expansion of order 2, we get

B=bs(xr,knots=c(2,5),Boundary.knots=c(0,10),degre=2)
matplot(xr,B,type="l")
abline(v=c(0,2,5,10),lty=2)

So, what’s going on? If we consider the constant, and the first component of the spline based matrix, we get

k=2
plot(db)
B=cbind(1,B)
lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1)

If we add the constant term, the first term and the second term, we get the part on the left, before the first knot,

k=3
lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1)

and with three terms from the spline based matrix, we can get the part between the two knots,

k=4
lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1)

and finallty, when we sum all the terms, we get this time the part on the right, after the last knot,

k=5
lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1)

This is what we get using a spline regression, quadratic, with two (fixed) knots. And can can even get confidence intervals, as before

reg=lm(yr~B)
P=predict(reg,interval="confidence")
plot(db,col="white")
polygon(c(xr,rev(xr)),c(P[,2],rev(P[,3])),col="light blue",border=NA)
points(db)
reg=lm(yr~B)
lines(xr,P[,1],col="red")
abline(v=c(0,2,5,10),lty=2)

The great idea here is to use functions $(x-x_i)_+$, that will insure continuity at point $x_i$.

Of course, we can use those splines on our Dexter application,

Here again, using linear spline function, it is possible to impose a continuity constraint,

plot(data$no,data$mu,ylim=c(6,10))
abline(v=12*(0:8)+.5,lty=2)
reg=lm(mu~bs(no,knots=c(12*(1:7)+.5),Boundary.knots=c(0,97),
degre=1),data=db)
lines(c(1:94,96),predict(reg),col="red")

But we can also consider some quadratic splines,

plot(data$no,data$mu,ylim=c(6,10))
abline(v=12*(0:8)+.5,lty=2)
reg=lm(mu~bs(no,knots=c(12*(1:7)+.5),Boundary.knots=c(0,97),
degre=2),data=db)
lines(c(1:94,96),predict(reg),col="red")

# Some heuristics about local regression and kernel smoothing

In a standard linear model, we assume that $\mathbb{E}(Y\vert X=x)=\beta_0+\beta_1 x$. Alternatives can be considered, when the linear assumption is too strong.

• Polynomial regression

A natural extension might be to assume some polynomial function,

$\mathbb{E}(Y\vert X=x)=\beta_0+\beta_1 x+\beta_2 x^2 +\cdots +\beta_k x^k$

Again, in the standard linear model approach (with a conditional normal distribution using the GLM terminology), parameters $\boldsymbol{\beta}=(\beta_0,\beta_1,\cdots,\beta_k)$ can be obtained using least squares, where a regression of $Y$ on $\boldsymbol{X}=(1,X,X^2,\cdots,X^k)$ is considered.

Even if this polynomial model is not the real one, it might still be a good approximation for $\mathbb{E}(Y\vert X=x)=h(x)$. Actually, from Stone-Weierstrass theorem, if $h(\cdot)$ is continuous on some interval, then there is a uniform approximation of $h(\cdot)$ by polynomial functions.

Just to illustrate, consider the following (simulated) dataset

set.seed(1)
n=10
xr = seq(0,n,by=.1)
yr = sin(xr/2)+rnorm(length(xr))/2
db = data.frame(x=xr,y=yr)
plot(db)

with the standard regression line

reg = lm(y ~ x,data=db)
abline(reg,col="red")

Consider some polynomial regression. If the degree of the polynomial function is large enough, any kind of pattern can be obtained,

reg=lm(y~poly(x,5),data=db)

But if the degree is too large, then too many ‘oscillations’ are obtained,

reg=lm(y~poly(x,25),data=db)

and the estimation might be be seen as no longer robust: if we change one point, there might be important (local) changes

plot(db)
attach(db)
lines(xr,predict(reg),col="red",lty=2)
yrm=yr;yrm[31]=yr[31]-2
regm=lm(yrm~poly(xr,25))
lines(xr,predict(regm),col="red")

• Local regression

Actually, if our interest is to have locally a good approximation of  $h(\cdot)$, why not use a local regression?

This can be done easily using a weighted regression, where, in the least square formulation, we consider

$\min\left\{ \sum_{i=1}^n \omega_i [Y_i-(\beta_0+\beta_1 X_i)]^2 \right\}$

(it is possible to consider weights in the GLM framework, but let’s keep that for another post). Two comments here:

• here I consider a linear model, but any polynomial model can be considered. Even a constant one. In that case, the optimization problem is

$\min\left\{ \sum_{i=1}^n \omega_i [Y_i-\beta_0]^2 \right\}$which can be solve explicitly, since

$\widehat{\beta}_0=\frac{\sum \omega_i Y_i}{\sum \omega_i}$

• so far, nothing was mentioned about the weights. The idea is simple, here: if you can a good prediction at point $x_0$, then $\omega_i$ should be proportional to some distance between $X_i$ and $x_0$: if $X_i$ is too far from $x_0$, then it should not have to much influence on the prediction.

For instance, if we want to have a prediction at some point $x_0$, consider $\omega_i\propto \boldsymbol{1}(\vert X_i-x_0 \vert<1)$. With this model, we remove observations too far away,

Actually, here, it is the same as

reg=lm(yr~xr,subset=which(abs(xr-x0)<1)

A more general idea is to consider some kernel function $K(\cdot)$ that gives the shape of the weight function, and some bandwidth (usually denoted h) that gives the length of the neighborhood, so that

$\omega_i = K\left(\frac{x_0-X_i}{b}\right)$

This is actually the so-called Nadaraya-Watson estimator of function $h(\cdot)$.
In the previous case, we did consider a uniform kernel $K(x)=\boldsymbol{1}(x\in[-1/2,+1/2])$, with bandwith $2$,

But using this weight function, with a strong discontinuity may not be the best idea… Why not a Gaussian kernel,

$K(x)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right)$

This can be done using

fitloc0 = function(x0){
w=dnorm((xr-x0))
reg=lm(y~1,data=db,weights=w)
return(predict(reg,newdata=data.frame(x=x0)))}

On our dataset, we can plot

ul=seq(0,10,by=.01)
vl0=Vectorize(fitloc0)(ul)
u0=seq(-2,7,by=.01)
linearlocalconst=function(x0){
w=dnorm((xr-x0))
plot(db,cex=abs(w)*4)
lines(ul,vl0,col="red")
axis(3)
axis(2)
reg=lm(y~1,data=db,weights=w)
u=seq(0,10,by=.02)
v=predict(reg,newdata=data.frame(x=u))
lines(u,v,col="red",lwd=2)
abline(v=c(0,x0,10),lty=2)
}
linearlocalconst(2)

Here, we want a local regression at point 2. The horizonal line below is the regression (the size of the point is proportional to the wieght). The curve, in red, is the evolution of the local regression

Let us use an animation to visualize the construction of the curve. One can use

library(animate)

but for some reasons, I cannot install the package easily on Linux. And it is not a big deal. We can still use a loop to generate some graphs

vx0=seq(1,9,by=.1)
vx0=c(vx0,rev(vx0))
graphloc=function(i){
name=paste("local-reg-",100+i,".png",sep="")
png(name,600,400)
linearlocalconst(vx0[i])
dev.off()}

for(i in 1:length(vx0)) graphloc(i)

and then, in a terminal, I simply use

    convert -delay 25 /home/freak/local-reg-1*.png /home/freak/local-reg.gif

Of course, it is possible to consider a linear model, locally,

fitloc1 = function(x0){
w=dnorm((xr-x0))
reg=lm(y~poly(x,degree=1),data=db,weights=w)
return(predict(reg,newdata=data.frame(x=x0)))}

or even a quadratic (local) regression,

fitloc2 = function(x0){
w=dnorm((xr-x0))
reg=lm(y~poly(x,degree=2),data=db,weights=w)
return(predict(reg,newdata=data.frame(x=x0)))}

Of course, we can change the bandwidth

To conclude the technical part this post, observe that, in practise, we have to choose the shape of the weight function (the so-called kernel). But there are (simple) technique to select the “optimal” bandwidth h. The idea of cross validation is to consider

$\min\left\{ \sum_{i=1}^n [Y_i-\widehat{Y}_i(b)]^2 \right\}$

where $\widehat{Y}_i(b)$ is the prediction obtained using a local regression technique, with bandwidth $b$. And to get a more accurate (and optimal) bandwith $\widehat{Y}_i(b)$ is obtained using a model estimated on a sample where the ith observation was removed. But again, that is not the main point in this post, so let’s keep that for another one…

Perhaps we can try on some real data? Inspired from a great post on http://f.briatte.org/teaching/ida/092_smoothing.html, by François Briatte, consider the Global Episode Opinion Survey, from some TV show, http://geos.tv/index.php/index?sid=189 , like Dexter.

library(XML)
file = "geos-tww.csv"
html = htmlParse("http://www.geos.tv/index.php/list?sid=189&collection=all")
html = xpathApply(html, "//table[@id='collectionTable']")[[1]]
data = data[,-3]
names(data)=c("no",names(data)[-1])
data=data[-(61:64),]

Let us reshape the dataset,

data$no = 1:96 data$mu = as.numeric(substr(as.character(data$Mean), 0, 4)) data$se =  sd(data$mu,na.rm=TRUE)/sqrt(as.numeric(as.character(data$Count)))
data$season = 1 + (data$no - 1)%/%12
data$season = factor(data$season)
plot(data$no,data$mu,ylim=c(6,10))
segments(data$no,data$mu-1.96*data$se, data$no,data$mu+1.96*data$se,col="light blue")

As done by François, we compute some kind of standard error, just to reflect uncertainty. But we won’t really use it.

plot(data$no,data$mu,ylim=c(6,10))
abline(v=12*(0:8)+.5,lty=2)
for(s in 1:8){reg=lm(mu~no,data=db,subset=season==s)
lines((s-1)*12+1:12,predict(reg)[1:12],col="red") }

Henre, we assume that all seasons should be considered as completely independent… which might not be a great assumption.

db = data
NW = ksmooth(db$no,db$mu,kernel = "normal",bandwidth=5)
plot(data$no,data$mu)
lines(NW,col="red")

We can try to look the curve with a larger bandwidth. The problem is that there is a missing value, at the end. If we (arbitrarily) fill it, we can run a kernel regression,

db$mu[95]=7 NW = ksmooth(db$no,db$mu,kernel = "normal",bandwidth=12) plot(data$no,data\$mu,ylim=c(6,10))
lines(NW,col="red")