---
title: 'GAM'
output: html_document
date: "2023-11-15"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r simulated data}
n=200
x=sort(runif(n,0,1))
x=seq(0,1,length=n)
truef=function(x) sin(2*pi*x)
my=truef(x)
signaltonoise=2

simula=function(signalnoise=signaltonoise){
  y=rnorm(n,my,sd=sd(my)/signalnoise)
  return(data.frame(x=x,y=y,truef=my))
}
d=simula()
plot(d$x,d$y)
```

```{r}
library(mgcv)
fit=gam(y~s(x),data=d)
fit3=gam(y~s(x,k=3),data=d)
fit50=gam(y~s(x,k=50),data=d)
plot(fit)
```


```{r}
plot(fit,residuals=TRUE,pch=20)
```

- partial residuals are weighted working residuals from PIRLS added to the spline (systematic departure from fj indicates a problem)
- Rug plot shows values of predictors.
- EDF for term reported in y axis lable.
- 95% Bayesian CIs shown.

Some more options for {\tt plot.gam}

```{r}
plot(fit,shade=TRUE,seWithMean = TRUE,scale=0)
```


```{r}
gam.check(fit)
```
## choice of k

Increase k to see 
```{r}
fit=gam(y~s(x,k=10),data=d)
fit=gam(y~s(x,k=20),data=d)
```

Check for unmodeled patterns in the residuals

```{r}
fit=gam(y~s(x,k=3),data=d)
rsd <- residuals(fit,type="deviance")
fit.res=gam(rsd~s(x,k=20)-1,data=d,select=TRUE)
```

```{r}
fit=gam(y~s(x,k=20),data=d)
rsd <- residuals(fit,type="deviance")
fit.res=gam(rsd~s(x,k=20)-1,data=d,select=TRUE)
```


## Smoothness selection

```{r}
fit1=gam(y~s(x,k=20),data=d,method="GCV.Cp")
fit2=gam(y~s(x,k=20),data=d,method="REML")
fit3=gam(y~s(x,k=20),data=d,method="ML")
plot(d$x,d$y)
curve(predict(fit1,newdata=data.frame(x=x)),add=TRUE)
curve(predict(fit2,newdata=data.frame(x=x)),add=TRUE,col="red")
curve(predict(fit3,newdata=data.frame(x=x)),add=TRUE,col="green")
```

## Model inference

```{r}
summary(fit)
```

```{r}
anova(fit)
```

```{r}
beta.hat=coef(fit)
beta.vcov=vcov(fit)
beta.vcov.freq=vcov(fit,freq=TRUE)
```


## More (smooth) covariates

```{r}
dat <- gamSim(1,n=400,dist="normal",scale=2)
fit=gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
```
```{r}
par(mfrow=c(1,4))
plot(fit)
```

```{r}
par(mfrow=c(1,4))
plot(fit,scale=0)
```

```{r}
vis.gam(fit,view=c("x1","x2"),theta=30,se=2)
```