Chapter 7 Ejemplo 9

library(fda)

## Loading required package: splines

## Loading required package: Matrix

## Loading required package: fds

## Loading required package: rainbow

## Loading required package: MASS

## Loading required package: pcaPP

## Loading required package: RCurl

## 
## Attaching package: 'fda'

## The following object is masked from 'package:graphics':
## 
##     matplot

library(fda.usc)

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-36. For overview type 'help("mgcv-package")'.

## ----------------------------------------------------------------------------------
##  Functional Data Analysis and Utilities for Statistical Computing
##  fda.usc version 2.0.2 (built on 2020-02-17) is now loaded
##  fda.usc is running sequentially usign foreach package
##  Please, execute ops.fda.usc() once to run in local parallel mode
##  Deprecated functions: min.basis, min.np, anova.hetero, anova.onefactor, anova.RPm
##  New functions: optim.basis, optim.np, fanova.hetero, fanova.onefactor, fanova.RPm
## ----------------------------------------------------------------------------------

library(rainbow)
library(mgcv)
require(freqdom.fda)

## Loading required package: freqdom.fda

## Loading required package: mvtnorm

## Loading required package: freqdom

## Loading required package: matrixcalc

## 
## Attaching package: 'freqdom'

## The following object is masked from 'package:base':
## 
##     %*%

require(fdANOVA)

## Loading required package: fdANOVA

require(fds)
require(fdcov)

## Loading required package: fdcov

require(changepoint)

## Loading required package: changepoint

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Successfully loaded changepoint package version 2.2.2
##  NOTE: Predefined penalty values changed in version 2.2.  Previous penalty values with a postfix 1 i.e. SIC1 are now without i.e. SIC and previous penalties without a postfix i.e. SIC are now with a postfix 0 i.e. SIC0. See NEWS and help files for further details.

require(reshape2)

## Loading required package: reshape2

require(fChange)

## Loading required package: fChange

require(roahd)

## Loading required package: roahd

## 
## Attaching package: 'roahd'

## The following object is masked from 'package:fda':
## 
##     fbplot

require(refund)

## Loading required package: refund

library(rriskDistributions)

data("CanadianWeather")
f=function(x){
  n=length(x)
  i=0
  for (ii in 1:n) {
    if(x[ii]<=0){
      i=i+1
    }
  }
  return(i)
}

La variable respuesta son los conteos de días que tuvieron temperaturas menores a 0 para cada ciudad.

x=fdata(t(CanadianWeather$dailyAv[,,2]))
CanadianWeather$y=as.data.frame(apply(CanadianWeather$dailyAv[,,"Temperature.C"], 2, f))
tt=x[["argvals"]]
response=CanadianWeather$y
response=as.data.frame(response)
colnames(response)=c("response")
nbasis.x=5
nbasis.b=5
basis1=create.bspline.basis(rangeval=range(tt),nbasis=nbasis.x)
basis2=create.bspline.basis(rangeval=range(tt),nbasis=nbasis.b)
basis.x=list("x"=basis1)
basis.b=list("x"=basis2)
ldata=list("df"=response,"x"=x)
res2=fregre.gsam(formula=response~s(x,k=5),family=poisson(),data=ldata,basis.b=basis.b,basis.x=basis.x)
summary(res2)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## response ~ +s(x.bspl4.1, k = 5) + s(x.bspl4.2, k = 5) + s(x.bspl4.3, 
##     k = 5) + s(x.bspl4.4, k = 5) + s(x.bspl4.5, k = 5)
## <environment: 0x000001fa51a1aab8>
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.69178    0.03139   149.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                edf Ref.df Chi.sq  p-value    
## s(x.bspl4.1) 4.000  4.000  38.21  < 2e-16 ***
## s(x.bspl4.2) 3.742  3.936  65.42  < 2e-16 ***
## s(x.bspl4.3) 3.999  4.000  25.55 3.92e-05 ***
## s(x.bspl4.4) 3.707  3.932  62.81  < 2e-16 ***
## s(x.bspl4.5) 2.773  3.130  17.83 0.000299 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.857   Deviance explained = 94.8%
## UBRE = 1.9094  Scale est. = 1         n = 35

plot.fdata(x,type="l",col=rainbow(10))

plot(res2$residuals,ylab="", main="Residuales del modelo")
abline(h=0)

mean((res2$residuals)^2)

## [1] 0.06785139

anova(res2)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## response ~ +s(x.bspl4.1, k = 5) + s(x.bspl4.2, k = 5) + s(x.bspl4.3, 
##     k = 5) + s(x.bspl4.4, k = 5) + s(x.bspl4.5, k = 5)
## <environment: 0x000001fa51a1aab8>
## 
## Approximate significance of smooth terms:
##                edf Ref.df Chi.sq  p-value
## s(x.bspl4.1) 4.000  4.000  38.21  < 2e-16
## s(x.bspl4.2) 3.742  3.936  65.42  < 2e-16
## s(x.bspl4.3) 3.999  4.000  25.55 3.92e-05
## s(x.bspl4.4) 3.707  3.932  62.81  < 2e-16
## s(x.bspl4.5) 2.773  3.130  17.83 0.000299

Xie, Yihui. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. http://yihui.org/knitr/.

———. 2021. Bookdown: Authoring Books and Technical Documents with r Markdown. https://CRAN.R-project.org/package=bookdown.