+ - 0:00:00
Notes for current slide
Notes for next slide

AlwaysR, Módulo III: Estadística en R

Clase 5: Introducción a modelos lineales y aditivos generalizados (GLM y GAM)

Dr. Giancarlo M. Correa

1 / 59

Modelos lineales generalizados (GLM)

2 / 59

Librerías a utilizar

A lo largo de esta clase, utilizaremos:

library(sjPlot)
library(sjmisc)
library(ggplot2)
library(faraway)
library(Sleuth3)
library(nlme)
library(lme4)
library(statmod)
library(MASS)
library(GLMsData)
library(ggeffects)
library(psych)
require(visibly)
3 / 59

GLM

Un modelo lineal simple puede ser expresado como:

yN(μ,σ2)

μ=β0+β1x1+β2x2+...+βpxp

4 / 59

GLM

Un modelo lineal simple puede ser expresado como:

yN(μ,σ2)

μ=β0+β1x1+β2x2+...+βpxp

Un GLM es una opción cuando la varianza de la variable respuesta no es constante o no está normalmente distribuida.

4 / 59

GLM

Un modelo lineal simple puede ser expresado como:

yN(μ,σ2)

μ=β0+β1x1+β2x2+...+βpxp

Un GLM es una opción cuando la varianza de la variable respuesta no es constante o no está normalmente distribuida.

GLM transforma la variable respuesta para que el método de estimación de parámetros siga funcionando bien (mediante una función link).

4 / 59

GLM

Un modelo lineal simple puede ser expresado como:

yN(μ,σ2)

μ=β0+β1x1+β2x2+...+βpxp

Un GLM es una opción cuando la varianza de la variable respuesta no es constante o no está normalmente distribuida.

GLM transforma la variable respuesta para que el método de estimación de parámetros siga funcionando bien (mediante una función link).

Al momento de implementar un GLM, debemos decidir dos cosas:

  1. La distribución a utilizar
  2. La función link (esto normalmente está por defecto en R).

yEDM(μ,etc)

g(μ)=β0+j=1pβjxj

4 / 59

GLM

Por ejemplo, podemos tener:

data(quilpie)
head(quilpie)
## Year Rain SOI Phase Exceed y
## 1 1921 38.4 2.7 2 Yes 1
## 2 1922 0.0 2.0 5 No 0
## 3 1923 0.0 -10.7 3 No 0
## 4 1924 24.4 6.9 2 Yes 1
## 5 1925 0.0 -12.5 3 No 0
## 6 1926 9.1 -1.0 4 No 0

Para este caso, si y es la variable respuesta, entonces podemos usar:

yiBinom(μi,mi)

logμi1μi=β0+β1SOIi

5 / 59

GLM: respuesta binaria

Para una variable de respuesta binaria (e.g. 0/1 - presencia/asusencia - hembra/macho), tipicamente la modelamos con una función link tipo logit y función de varianza μ(1μ). La respuesta modelada son las probabilidades de registro predichas (predicted log odds) del evento de interés.

6 / 59

GLM: respuesta binaria

Para una variable de respuesta binaria (e.g. 0/1 - presencia/asusencia - hembra/macho), tipicamente la modelamos con una función link tipo logit y función de varianza μ(1μ). La respuesta modelada son las probabilidades de registro predichas (predicted log odds) del evento de interés.

head(hsb)
## id gender race ses schtyp prog read write math science socst
## 1 70 male white low public general 57 52 41 47 57
## 2 121 female white middle public vocation 68 59 53 63 61
## 3 86 male white high public general 44 33 54 58 31
## 4 141 male white high public vocation 63 44 47 53 56
## 5 172 male white middle public academic 47 52 57 53 61
## 6 113 male white middle public academic 44 52 51 63 61
6 / 59

GLM: respuesta binaria

Para una variable de respuesta binaria (e.g. 0/1 - presencia/asusencia - hembra/macho), tipicamente la modelamos con una función link tipo logit y función de varianza μ(1μ). La respuesta modelada son las probabilidades de registro predichas (predicted log odds) del evento de interés.

head(hsb)
## id gender race ses schtyp prog read write math science socst
## 1 70 male white low public general 57 52 41 47 57
## 2 121 female white middle public vocation 68 59 53 63 61
## 3 86 male white high public general 44 33 54 58 31
## 4 141 male white high public vocation 63 44 47 53 56
## 5 172 male white middle public academic 47 52 57 53 61
## 6 113 male white middle public academic 44 52 51 63 61

Para este ejemplo, modelaremos el registro de hsb$prog == academic.

6 / 59

GLM: respuesta binaria

En R vamos a usar:

hsb$obj_var = ifelse(test = hsb$prog == 'academic', yes = 1, no = 0)
mod1 = glm(obj_var ~ read, family = binomial, data = hsb)
7 / 59

GLM: respuesta binaria

En R vamos a usar:

hsb$obj_var = ifelse(test = hsb$prog == 'academic', yes = 1, no = 0)
mod1 = glm(obj_var ~ read, family = binomial, data = hsb)

Aquí, la variable respuesta es una variable indicadora (binaria):

hsb$obj_var[1:20]
## [1] 0 0 0 0 1 1 0 1 0 1 0 1 1 1 1 0 1 0 1 0
7 / 59
summary(mod1)
##
## Call:
## glm(formula = obj_var ~ read, family = binomial, data = hsb)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8817 -1.0329 0.4935 1.0182 1.8654
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.67739 0.90426 -5.173 2.31e-07 ***
## read 0.09207 0.01730 5.321 1.03e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 276.76 on 199 degrees of freedom
## Residual deviance: 241.58 on 198 degrees of freedom
## AIC: 245.58
##
## Number of Fisher Scoring iterations: 4
8 / 59

GLM: respuesta binaria

Vemos un par de términos nuevos:

  • Null deviance: nos indica que tan bien la variable respuesta puede ser predicha por un modelo que solo contenga el intercepto. Un valor pequeño nos indica que el modelo explica muy bien la variable respuesta.

  • Residual deviance: nos indica que tan bien la variable respuesta es predicha por el modelo implementado. Un valor pequeño nos indica que el modelo explica muy bien la variable respuesta.

9 / 59

GLM: respuesta binaria

Vemos un par de términos nuevos:

  • Null deviance: nos indica que tan bien la variable respuesta puede ser predicha por un modelo que solo contenga el intercepto. Un valor pequeño nos indica que el modelo explica muy bien la variable respuesta.

  • Residual deviance: nos indica que tan bien la variable respuesta es predicha por el modelo implementado. Un valor pequeño nos indica que el modelo explica muy bien la variable respuesta.

Se recomienda explorar ?family para conocer las funciones link de cada distribución.

9 / 59

GLM: respuesta binaria

La revisión de supuestos en un GLM no es tan sencilla. Podemos revisar la independencia de observaciones (especialmente importante cuando tenemos datos que vienen con posible dependencia temporal o espacial):

plot(statmod::qresid(mod1))

10 / 59

GLM: residuos

11 / 59

GLM: respuesta binaria

Revisar los valores observados y ajustados:

plot_model(model = mod1, type = 'pred', show.data = TRUE)
## $read

12 / 59

GLM: respuesta binaria

Revisar distribucion asumida:

qqnorm(qresid(mod1))
qqline(qresid(mod1))

13 / 59

GLM: respuesta binaria

Revisar observaciones influyentes:

plot(cooks.distance(mod1), type='h')

14 / 59

GLM: datos de conteo

Cuando tenemos una variable discreta (e.g. conteo de individuos, de especies, etc.), es recomendable implementar el modelo usando la familia de distribuciones Poisson.

15 / 59

GLM: datos de conteo

Cuando tenemos una variable discreta (e.g. conteo de individuos, de especies, etc.), es recomendable implementar el modelo usando la familia de distribuciones Poisson.

head(cases)
## Days Students
## 1 1 6
## 2 2 8
## 3 3 12
## 4 3 9
## 5 4 3
## 6 4 3
15 / 59

GLM: datos de conteo

Exploremos la que será la variable respuesta en nuestro modelo:

hist(cases$Students)

16 / 59
mod2 = glm(Students ~ Days, family = poisson, data = cases)
summary(mod2)
##
## Call:
## glm(formula = Students ~ Days, family = poisson, data = cases)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.00482 -0.85719 -0.09331 0.63969 1.73696
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.990235 0.083935 23.71 <2e-16 ***
## Days -0.017463 0.001727 -10.11 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 215.36 on 108 degrees of freedom
## Residual deviance: 101.17 on 107 degrees of freedom
## AIC: 393.11
##
## Number of Fisher Scoring iterations: 5
17 / 59

GLM: datos de conteo

Verificar independencia de residuos:

plot(qresid(mod2))

18 / 59

GLM: datos de conteo

Revisar los valores observados y ajustados:

plot_model(model = mod2, type = 'pred', show.data = TRUE)
## $Days

19 / 59

GLM: datos de conteo

Revisar distribucion asumida:

qqnorm(qresid(mod2))
qqline(qresid(mod2))

20 / 59

GLM: datos de conteo

Revisar observaciones influyentes:

plot(cooks.distance(mod2), type='h')

21 / 59

GLM: datos de conteo con sobredispersión

Exploremos una nueva base de datos:

head(disc)
## count year yearSqr
## 1 5 0 0
## 2 3 1 1
## 3 0 2 4
## 4 2 3 9
## 5 0 4 16
## 6 3 5 25
22 / 59

GLM: datos de conteo con sobredispersión

Exploremos la que ahora será nuestra variable explicativa:

hist(disc$count)

23 / 59
mod3 = glm(count ~ year + yearSqr, family = "poisson", data = disc)
summary(mod3)
##
## Call:
## glm(formula = count ~ year + yearSqr, family = "poisson", data = disc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9066 -0.8397 -0.2544 0.4776 3.3303
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.592e-01 1.814e-01 4.186 2.84e-05 ***
## year 3.356e-02 8.499e-03 3.948 7.87e-05 ***
## yearSqr -4.106e-04 8.699e-05 -4.720 2.35e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 164.68 on 99 degrees of freedom
## Residual deviance: 132.84 on 97 degrees of freedom
## AIC: 407.85
##
## Number of Fisher Scoring iterations: 5
24 / 59
mod4 = glm(count ~ year + yearSqr, family = "quasipoisson", data = disc)
summary(mod4)
##
## Call:
## glm(formula = count ~ year + yearSqr, family = "quasipoisson",
## data = disc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9066 -0.8397 -0.2544 0.4776 3.3303
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7592473 0.2072715 3.663 0.000406 ***
## year 0.0335569 0.0097112 3.455 0.000816 ***
## yearSqr -0.0004106 0.0000994 -4.131 7.66e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.305649)
##
## Null deviance: 164.68 on 99 degrees of freedom
## Residual deviance: 132.84 on 97 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
25 / 59

GLM: datos de conteo con sobredispersión

El parámetro de overdispersion ( ϕ ) nos indica si una distribución quasipoisson es adecuada. Si ϕ>1, entonces nos sugiere usa una quasipoisson, si ϕ<1, una quasipoisson no es necesaria.

26 / 59

GLM: datos de conteo con sobredispersión

Revisamos los residuos:

qqnorm(qresid(mod4))
qqline(qresid(mod4))

27 / 59

GLM: binomial negativa

mod5 = MASS::glm.nb(count ~ year + yearSqr, data = disc)
summary(mod5)
##
## Call:
## MASS::glm.nb(formula = count ~ year + yearSqr, data = disc, init.theta = 11.53157519,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6849 -0.7707 -0.2344 0.4232 2.7134
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.526e-01 2.019e-01 3.727 0.000194 ***
## year 3.386e-02 9.471e-03 3.575 0.000350 ***
## yearSqr -4.132e-04 9.627e-05 -4.292 1.77e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(11.5316) family taken to be 1)
##
## Null deviance: 131.91 on 99 degrees of freedom
## Residual deviance: 106.20 on 97 degrees of freedom
## AIC: 406.21
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 11.53
## Std. Err.: 7.42
##
## 2 x log-likelihood: -398.214
28 / 59

GLM: datos de conteo con sobredispersión

Podemos decidir usar una distribución binomial negativa cuando, en un GLM Poisson, el ratio de Residual deviance/degrees of freedom es mayor a 1.

29 / 59

GLM: datos de conteo con sobredispersión

Podemos decidir usar una distribución binomial negativa cuando, en un GLM Poisson, el ratio de Residual deviance/degrees of freedom es mayor a 1.

Revisamos los residuos:

qqnorm(qresid(mod5))
qqline(qresid(mod5))

29 / 59

GLM: respuesta gamma

Dada nuestra base de datos:

head(case1202)
## Bsal Sal77 Sex Senior Age Educ Exper
## 1 5040 12420 Male 96 329 15 14.0
## 2 6300 12060 Male 82 357 15 72.0
## 3 6000 15120 Male 67 315 15 35.5
## 4 6000 16320 Male 97 354 12 24.0
## 5 6000 12300 Male 66 351 12 56.0
## 6 6840 10380 Male 92 374 15 41.5
30 / 59

GLM: respuesta gamma

Modelo lineal simple:

mod6 = lm(Exper ~ Age, data=case1202)
summary(mod6)
##
## Call:
## lm(formula = Exper ~ Age, data = case1202)
##
## Residuals:
## Min 1Q Median 3Q Max
## -150.974 -21.633 -3.205 24.866 146.490
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -144.59165 20.26830 -7.134 2.28e-10 ***
## Age 0.51754 0.04099 12.626 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 55.13 on 91 degrees of freedom
## Multiple R-squared: 0.6366, Adjusted R-squared: 0.6326
## F-statistic: 159.4 on 1 and 91 DF, p-value: < 2.2e-16
31 / 59

GLM: respuesta gamma

Veamos los valores ajustados:

plot_model(mod6, type = 'pred', show.data = TRUE)
## $Age

32 / 59

GLM: respuesta gamma

Exploremos la distribución de la variable respuesta:

hist(case1202$Exper)

33 / 59

GLM: respuesta gamma

¿Qué alternativas tenemos?

34 / 59

GLM: respuesta gamma

¿Qué alternativas tenemos?

Podemos implementar un GLM con respuesta gamma (solo valores positivos):

mod7 = glm(Exper + 1 ~ Age, data = case1202, family = Gamma(link = "log"))
34 / 59
summary(mod7)
##
## Call:
## glm(formula = Exper + 1 ~ Age, family = Gamma(link = "log"),
## data = case1202)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.22397 -0.60604 -0.02406 0.34153 1.07844
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4770108 0.2118855 6.971 4.87e-10 ***
## Age 0.0060060 0.0004285 14.016 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.3321014)
##
## Null deviance: 101.471 on 92 degrees of freedom
## Residual deviance: 46.044 on 91 degrees of freedom
## AIC: 970.8
##
## Number of Fisher Scoring iterations: 6
35 / 59

GLM: respuesta gamma

plot_model(mod7, type = 'pred', show.data = TRUE)
## $Age

36 / 59

Modelos aditivos generalizados (GAM)

37 / 59

GAM

Hay casos en donde la asociación lineal entre las variables no se cumple. Para casos donde no se tenga una relación lineal entre la variable respuesta y variables independientes, podemos utilizar un modelo aditivo generalizado:

yEDM(μ,etc)

g(μ)=β0+f(x1)+f(x2)+...+f(xp)

E(y)μ

38 / 59

GAM

Hay casos en donde la asociación lineal entre las variables no se cumple. Para casos donde no se tenga una relación lineal entre la variable respuesta y variables independientes, podemos utilizar un modelo aditivo generalizado:

yEDM(μ,etc)

g(μ)=β0+f(x1)+f(x2)+...+f(xp)

E(y)μ Variables independientes ahora incorporan una función de suavizamiento (basis function, f ), lo que permite relaciones no lineales entre las variables.

38 / 59

Por ejemplo:

39 / 59

¿Podríamos implementar un ANOVA o algún modelo lineal?

40 / 59

Una función como esta nos puede ayudar mucho:

41 / 59

42 / 59

¿Qué es un nodo (knots)?

43 / 59

¿Qué es un nodo (knots)?

A mayor cantidad de nodos, mejor ajuste (cuidado con oversmoothing!).

43 / 59

De una forma general:

g(μ)=β0+f(x1)+f(x2)+...+f(xp)

44 / 59

De una forma general:

g(μ)=β0+f(x1)+f(x2)+...+f(xp)

f(x1)=β0+β1x1+...+βdxd

44 / 59

De una forma general:

g(μ)=β0+f(x1)+f(x2)+...+f(xp)

f(x1)=β0+β1x1+...+βdxd Podemos producir cualquier función polinomial!

44 / 59

GAM en R

library(mgcv)
## This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.
isit = read.csv('ISIT.csv')
isit$Season <- as.factor(isit$Season)
head(isit)
## SampleDepth Sources Station Time Latitude Longitude Xkm Ykm Month Year
## 1 517 28.73 1 3 50.1508 -14.4792 -34.106 16.779 4 2001
## 2 582 27.90 1 3 50.1508 -14.4792 -34.106 16.779 4 2001
## 3 547 23.44 1 3 50.1508 -14.4792 -34.106 16.779 4 2001
## 4 614 18.33 1 3 50.1508 -14.4792 -34.106 16.779 4 2001
## 5 1068 12.38 1 3 50.1508 -14.4792 -34.106 16.779 4 2001
## 6 1005 11.23 1 3 50.1508 -14.4792 -34.106 16.779 4 2001
## BottomDepth Season Discovery RelativeDepth
## 1 3939 1 252 3422
## 2 3939 1 252 3357
## 3 3939 1 252 3392
## 4 3939 1 252 3325
## 5 3939 1 252 2871
## 6 3939 1 252 2934
45 / 59
fit_gam = gam(Sources ~ Season + s(SampleDepth), data = isit)
summary(fit_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Sources ~ Season + s(SampleDepth)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.2493 0.3613 20.07 <2e-16 ***
## Season2 6.1631 0.4826 12.77 <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(SampleDepth) 8.829 8.991 184.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.692 Deviance explained = 69.5%
## GCV = 42.991 Scale est. = 42.401 n = 789
46 / 59

GAM en R

Aquí hay algunos términos nuevos:

  • edf: effective degrees of freedom. Es un resultados de la penalidad aplicada al grado de suavizamiento (detalles estadísticos en bibliografía!). Evita un over-suavizamiento o under-suavizamiento. Puede interpretarse como el nivel de suavizamiento aplicado a una variable.

  • GCV: generalized cross validation. Es un proxy del error de predicción cuadrático medio (mean square prediction error). Esto se basa en un método llamado leave-one-out cross validation. Utilizado para comparar modelos (similar a AIC), donde el modelo con menor GCV es el mejor.

  • Deviance explained: es una generalización del R2 (sin ajustar).

  • Scale est.: equivalente a la suma de cuadrado de los residuales.

47 / 59
plot(fit_gam, pages=1)

48 / 59

También podemos observarlo en escala de las variable respuesta:

myEffect = ggeffects::ggpredict(fit_gam, terms = "Season")
plot(myEffect)

49 / 59

También podemos observarlo en escala de las variable respuesta:

myEffect = ggeffects::ggpredict(fit_gam, terms = "SampleDepth")
plot(myEffect)

50 / 59

GAM en R

Existen diferentes formas de suavizar una variable (smooth term), los cuales pueden ser explorados aquí. El método por defecto es: thin plate regression splines (funciona bien en la mayoría de los casos) bs = "tp".

51 / 59

También debemos revisar los supuestos de un GAM:

par(mfrow = c(2,2))
gam.check(fit_gam)

##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradient at convergence was 0.0002508607 .
## The Hessian was positive definite.
## Model rank = 11 / 11
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(SampleDepth) 9.00 8.83 1.1 1
52 / 59

También revisamos el nivel de suavizamiento que hemos usado usando gam.check(). Es decir, el número de nodos (k) en la función s(). Recuerden que el valor por defecto es s(..., k = 10).

## k' edf k-index p-value
## s(SampleDepth) 9 8.82882 1.10425 1
53 / 59

También revisamos el nivel de suavizamiento que hemos usado usando gam.check(). Es decir, el número de nodos (k) en la función s(). Recuerden que el valor por defecto es s(..., k = 10).

## k' edf k-index p-value
## s(SampleDepth) 9 8.82882 1.10425 1

Lo mas importante es ver el pvalue. Un valor pvalue muy pequeño nos dice que no hemos usado un k suficientemente alto.

53 / 59

GAM en R: respuesta binomial

Hasta el momento hemos visto una respuesta Gaussiana (normal) con una función link igual a la identidad (identity). En GAMs también podemos incorporar diferentes familias:

54 / 59

GAM en R: respuesta binomial

Hasta el momento hemos visto una respuesta Gaussiana (normal) con una función link igual a la identidad (identity). En GAMs también podemos incorporar diferentes familias:

cricketer = read.csv('cricketr.csv')
cricketer$hand = NA
cricketer$hand = ifelse(test = cricketer$left == 'left', yes = 1, no = 0)
head(cricketer)
## left year life dead acd kia inbed cause hand
## 1 right 1890 102 0 0 0 0 alive 0
## 2 left 1892 100 0 0 0 0 alive 1
## 3 right 1893 99 0 0 0 0 alive 0
## 4 right 1894 98 0 0 0 0 alive 0
## 5 right 1896 96 0 0 0 0 alive 0
## 6 right 1896 96 0 0 0 0 alive 0
54 / 59

Para implementar el modelo:

hand.gam = gam(hand ~ s(year), data=cricketer, family=binomial)
summary(hand.gam)
##
## Family: binomial
## Link function: logit
##
## Formula:
## hand ~ s(year)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.49285 0.03361 -44.41 <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(year) 7.873 8.664 23.87 0.0066 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.00301 Deviance explained = 0.472%
## UBRE = -0.044557 Scale est. = 1 n = 5960
55 / 59

Revisamos los efectos:

plot(ggeffects::ggpredict(hand.gam), facets = TRUE)

56 / 59

GAM en R: respuesta Poisson

Podemos transformar los datos de la siguiente forma:

rtlef = data.frame(with(cricketer, as(table(year, left),"matrix")))
rtlef$year = as.numeric(rownames(rtlef))
head(rtlef)
## left right year
## 1840 1 6 1840
## 1841 4 16 1841
## 1842 5 16 1842
## 1843 3 25 1843
## 1844 4 25 1844
## 1845 2 24 1845
57 / 59

GAM en R: respuesta Poisson

Podemos transformar los datos de la siguiente forma:

rtlef = data.frame(with(cricketer, as(table(year, left),"matrix")))
rtlef$year = as.numeric(rownames(rtlef))
head(rtlef)
## left right year
## 1840 1 6 1840
## 1841 4 16 1841
## 1842 5 16 1842
## 1843 3 25 1843
## 1844 4 25 1844
## 1845 2 24 1845
hand.gam2 = gam(left ~ s(year), data=rtlef, family=poisson)
57 / 59
plot(ggeffects::ggpredict(hand.gam2), facets = TRUE)

58 / 59
plot(ggeffects::ggpredict(hand.gam2), facets = TRUE)

Hay mucha información sobre GAMs. Se recomienda esta fuente.

58 / 59

Gracias!

Contacto: cursos@cousteau-group.com

59 / 59

Modelos lineales generalizados (GLM)

2 / 59
Paused

Help

Keyboard shortcuts

, , Pg Up, k Go to previous slide
, , Pg Dn, Space, j Go to next slide
Home Go to first slide
End Go to last slide
Number + Return Go to specific slide
b / m / f Toggle blackout / mirrored / fullscreen mode
c Clone slideshow
p Toggle presenter mode
t Restart the presentation timer
?, h Toggle this help
Esc Back to slideshow