intro: GLM e overdispersion
workshop hands-on
enrico toffalini x psicostat
31 ottobre 2023
rm(list=ls())
library(flipscores)
library(effects)
library(lmtest)
library(ggplot2)
ts = 20
tc = "#555555"________________________________
POISSON
Esempio - I bambini leggono un brano molto lungo e occasionalmente commettono qualche errore. Contiamo gli errori commessi da ciascuno: questa è la distribuzione. Chiaramente, a causa del processo generativo, non ha valori continui e non appare distribuita come una normale.
set.seed(2)
N = 500; x = rpois(N,1.4)
ggplot()+
geom_histogram(aes(y=after_stat(density),x=x),fill="darkblue",alpha=.7,binwidth=0.5)+
scale_x_continuous(breaks=seq(0,100,1))+
theme(text=element_text(size=ts,color=tc))+xlab("Errors")
_______
Aggiungiamo una variabile età. Vediamo quanti errori vengono commessi da bambini in un range di età diverse. Chiaramente decrescono, ma non con una progressione lineare. Si nota anche che valore medio e varianza sono legati (nella Poisson: M = VAR (SD²) = λ).
set.seed(1)
N = 250; x = round(runif(N,6,10),1); y = rpois(N,exp(6.5-x*.8))
d = data.frame(y,x)
ggplot(d,aes(x=x,y=y))+
geom_point(size=4,alpha=.5,color="darkblue")+
theme(text=element_text(size=ts,color=tc))+
scale_y_continuous(breaks=seq(0,20,2))+scale_x_continuous(breaks=seq(0,20,.5))+
ylab("Errors")+xlab("Age (years)")
_______
Una semplice regressione lineare (sotto in rosso) non è adeguata a descrivere il fenomeno: non stima correttamente il valore atteso (medio) lungo il continuo dell’età, presenta eteroschedasticità, ed esce perfino dal range dei valori osservati.
fitL = glm(y~x,data=d)
effL = data.frame(allEffects(fitL,xlevels=list(x=seq(min(x),max(x),.05)))$"x")
ggplot(d,aes(x=x,y=y))+
geom_point(size=4,alpha=.5,color="darkblue")+
geom_line(data=effL,aes(y=fit),size=2,color="darkred")+
geom_ribbon(data=effL,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="darkred")+
theme(text=element_text(size=ts,color=tc))+
scale_y_continuous(breaks=seq(0,20,2))+scale_x_continuous(breaks=seq(0,20,.5))+
ylab("Errors")+xlab("Age (years)")
_______
Un GLM con family “poisson” descrive il fenomeno decisamente meglio.
fitP = glm(y~x,data=d,family="poisson")
effP = data.frame(allEffects(fitP,xlevels=list(x=seq(min(x),max(x),.05)))$"x")
ggplot(d,aes(x=x,y=y))+
geom_point(size=4,alpha=.5,color="darkblue")+
geom_line(data=effP,aes(y=fit),size=2,color="blue")+
geom_ribbon(data=effP,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="blue")+
geom_line(data=effL,aes(y=fit),size=2,color="darkred")+
geom_ribbon(data=effL,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="darkred")+
theme(text=element_text(size=ts,color=tc))+
scale_y_continuous(breaks=seq(0,20,2))+scale_x_continuous(breaks=seq(0,20,.5))+
ylab("Errors")+xlab("Age (years)")
_______
In verità i due modelli qui sopra, su dati come questi, presentano stessa potenza e stesso rischio di falsi positivi (type I e type II error).
_______
interazione o no?
Prendiamo questo altro caso in cui abbiamo due effetti principali (Age e un fattore categoriale Z), ma NON la loro interazione. Le due rette sembrano avvicinarsi (sulla scala lineare) ma non si avvicinano sulla scala della funzione legame del modello che ha generato i dati (log), quindi NON c’è interazione.
set.seed(3)
N = 250; x = round(runif(N,6,10),1); z = rbinom(N,1,.5); y = rpois(N,exp(5.5-x*.5+z*.8))
d = data.frame(y,x,z=as.factor(z))
fitP = glm(y~x*z,data=d,family="poisson")
effP = data.frame(allEffects(fitP,xlevels=list(x=seq(min(x),max(x),.05)))$"x:z")
effP$z = as.factor(effP$z)
ggplot(d,aes(x=x,y=y,shape=z,color=z))+
geom_point(size=4,alpha=.6)+
scale_color_manual(values=c("darkorange2","darkgreen"))+
scale_fill_manual(values=c("darkorange1","darkgreen"))+
geom_line(data=effP,aes(x=x,y=fit,group=z,linetype=z),size=2)+
geom_ribbon(data=effP,aes(x=x,y=fit,ymin=lower,ymax=upper,group=z,fill=z),alpha=.3,color=NA)+
theme(text=element_text(size=ts,color=tc))+scale_x_continuous(breaks=seq(0,20,.5))+
ylab("Errors")+xlab("Age (years)")
fitP1 = glm(y ~ x * z, data=d, family="poisson")
fitP0 = glm(y ~ x + z, data=d, family="poisson")
lrtest(fitP1, fitP0)## Likelihood ratio test
##
## Model 1: y ~ x * z
## Model 2: y ~ x + z
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -601.25
## 2 3 -601.31 -1 0.116 0.7334Se però sbagliando usiamo il modello lineare normale, questo stima un’interazione Age x Z che non c’è nel processo generativo dei dati. C’è quindi un elevato rischio di errore di I tipo (falsi positivi) in questa situazione.
fitL = glm(y~x*z,data=d)
effL = data.frame(allEffects(fitL,xlevels=list(x=seq(min(x),max(x),.05)))$"x:z")
ggplot(d,aes(x=x,y=y,shape=z,color=z))+
geom_point(size=4,alpha=.6)+
scale_color_manual(values=c("darkorange2","darkgreen"))+
scale_fill_manual(values=c("brown1","darkslategray"))+
geom_line(data=effL,aes(x=x,y=fit,group=z,linetype=z),size=2,color=rep(c("brown1","darkslategray"),each=nrow(effL)/2))+
geom_ribbon(data=effL,aes(x=x,y=fit,ymin=lower,ymax=upper,group=z,fill=z),alpha=.3,color=NA)+
theme(text=element_text(size=ts,color=tc))+scale_x_continuous(breaks=seq(0,20,.5))+
ylab("Errors")+xlab("Age (years)")
fitL1 = glm(y ~ x * z, data=d)
fitL0 = glm(y ~ x + z, data=d)
lrtest(fitL1, fitL0)## Likelihood ratio test
##
## Model 1: y ~ x * z
## Model 2: y ~ x + z
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -654.17
## 2 4 -678.62 -1 48.899 2.695e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1_______
Overdispersion
Riprendiamo l’esempio sopra ma trascuriamo il fattore Z (c’è, pesa sul fenomeno, ma semplicemente non lo abbiamo misurato durante la raccolta dati). Teniamo conto solo dell’Age. Apparentemente il modello va bene.
set.seed(3)
N = 250; x = round(runif(N,6,10),1); z = rbinom(N,1,.5); y = rpois(N,exp(5.5-x*.5+z*.8))
d = data.frame(y,x)
fitP = glm(y~x,data=d,family="poisson")
effP = data.frame(allEffects(fitP,xlevels=list(x=seq(min(x),max(x),.05)))$"x")
ggplot(d,aes(x=x,y=y))+
geom_point(size=4,alpha=.5,color="darkblue")+
geom_line(data=effP,aes(y=fit),size=2,color="chartreuse4")+
geom_ribbon(data=effP,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="chartreuse4")+
theme(text=element_text(size=ts,color=tc))+
ylab("Errors")+xlab("Age (years)")
Tuttavia la distribuzione (dispersione) dei valori osservati è più ampia di quella implicata dai parametri del modello –> overdispersion. Questa condizione aumenta drasticamente il rischio di falsi positivi (errore del I tipo).
nsim = 100
sim = simulate(fitP,nsim=nsim)
sim = reshape(sim,varying=colnames(sim),v.names="score",direction="long")
freq = data.frame(Errors=c(sim$score,d$y),type=c(rep("simulated",length(sim$score)),rep("observed",length(d$y))))
pd = position_dodge(width=0.4)
ggplot(data=freq,aes(x=Errors,group=type,fill=type))+
geom_histogram(aes(y=after_stat(density),alpha=type),position=pd,binwidth=1)+
scale_fill_manual(values=c("darkblue","chartreuse4"))+
scale_alpha_manual(values=c(0.7,0.6))+
theme(text=element_text(size=ts,color=tc),axis.text.y=element_blank())
Possiamo anche testare statisticamente l’overdispersion.
performance::check_overdispersion(fitP)## # Overdispersion test
##
## dispersion ratio = 2.243
## Pearson's Chi-Squared = 556.262
## p-value = < 0.001Un buon trick per correggere gli errori standard (e ridurre il rischio di falsi positivi), nel caso di overdispersion nella poisson, è usare la “quasipoisson”.
fitQP = glm(y ~ x, data=d, family="quasipoisson"); summary(fitQP)##
## Call:
## glm(formula = y ~ x, family = "quasipoisson", data = d)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.79407 0.23880 24.26 <2e-16 ***
## x -0.47083 0.03183 -14.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.242991)
##
## Null deviance: 1106.77 on 249 degrees of freedom
## Residual deviance: 574.25 on 248 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5fitP = glm(y ~ x, data=d, family="poisson"); summary(fitP)##
## Call:
## glm(formula = y ~ x, family = "poisson", data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.79407 0.15945 36.34 <2e-16 ***
## x -0.47083 0.02125 -22.16 <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: 1106.77 on 249 degrees of freedom
## Residual deviance: 574.25 on 248 degrees of freedom
## AIC: 1516.1
##
## Number of Fisher Scoring iterations: 5________________________________
GAMMA
Altra situazione frequente è quella della distribuzione Gamma, che potrei trovare ad esempio nelle distribuzioni dei tempi di risposta
set.seed(1)
N = 1e4; x = rgamma(N,1.7)
ggplot()+
geom_density(aes(y=after_stat(density),x=x),fill="darkblue",alpha=.7)+
scale_x_continuous(breaks=seq(0,100,1),limits=c(0,8))+
theme(text=element_text(size=ts,color=tc))+xlab("RT")
Un trick spesso usato per “normalizzare” i punteggi è la trasformazione in logaritmo, che però funziona fino a un certo punto.
xlog = log(x)
ggplot()+
geom_density(aes(y=after_stat(density),x=xlog),fill="darkblue",alpha=.7)+
scale_x_continuous(breaks=seq(-100,100,1))+
theme(text=element_text(size=ts,color=tc))+xlab("RT log")
_______
Anche qui possiamo inventare la situazione di un’interazione che non c’è, ma che un semplice modello lineare troverebbe.
set.seed(2)
N = 250; x = runif(N,-2,2); z = rbinom(N,1,.5); y = rgamma(N,shape=exp(-x*.5+z*.8))
d = data.frame(y,x,z); d$z = as.factor(d$z)
fitG = glm(y~x*z,data=d,family=Gamma(link="log"))
effG = data.frame(allEffects(fitG,xlevels=list(x=seq(min(x),max(x),.05)))$"x:z")
effG$z = as.factor(effG$z)
ggplot(d,aes(x=x,y=y,shape=z,color=z))+
geom_point(size=4,alpha=.6)+
scale_color_manual(values=c("darkorange2","darkgreen"))+
scale_fill_manual(values=c("darkorange1","darkgreen"))+
geom_line(data=effG,aes(x=x,y=fit,group=z,linetype=z),size=2)+
geom_ribbon(data=effG,aes(x=x,y=fit,ymin=lower,ymax=upper,group=z,fill=z),alpha=.3,color=NA)+
theme(text=element_text(size=ts,color=tc))+
ylab("RT")+xlab("x")
#### GLM con gamma ####
fitG1 = glm(y ~ x * z, data=d, family=Gamma(link="log"))
fitG0 = glm(y ~ x + z, data=d, family=Gamma(link="log"))
lrtest(fitG1, fitG0)## Likelihood ratio test
##
## Model 1: y ~ x * z
## Model 2: y ~ x + z
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -365.10
## 2 4 -365.46 -1 0.72 0.3961#### Linear model/GLM normale ####
fitL1 = glm(y ~ x * z, data=d)
fitL0 = glm(y ~ x + z, data=d)
lrtest(fitL1, fitL0)## Likelihood ratio test
##
## Model 1: y ~ x * z
## Model 2: y ~ x + z
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -467.71
## 2 4 -480.41 -1 25.411 4.633e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1_______
Anche qui, però, si pone il problema dell’overdispersion (non in questo semplice caso simulato, a meno di “dimenticare” uno dei due predittori veri, ma in qualsiasi caso reale dove almeno alcuni predittori rilevanti non sono stati misurati o comunque inclusi nel modello).
set.seed(2)
N = 300
x = runif(N,-2,2)
z = rnorm(N,0,1)
w = rnorm(N,0,1)
y = rgamma(N,shape=exp(1-x*.6+z*.7+w*.7))
d = data.frame(y,x)
fitG = glm(y~x,data=d,family=Gamma(link="log"))
effG = data.frame(allEffects(fitG,xlevels=list(x=seq(min(x),max(x),.05)))$"x")
ggplot(d,aes(x=x,y=y))+
geom_point(size=4,alpha=.5,color="darkblue")+
geom_line(data=effG,aes(y=fit),size=2,color="chartreuse4")+
geom_ribbon(data=effG,aes(y=fit,ymin=lower,ymax=upper),alpha=.3,fill="chartreuse4")+
theme(text=element_text(size=ts,color=tc))+
ylab("RT")+xlab("x")
nsim = 100
sim = simulate(fitG,nsim=nsim)
sim = reshape(sim,varying=colnames(sim),v.names="score",direction="long")
freq = data.frame(RT=c(sim$score,d$y),type=c(rep("simulated",length(sim$score)),rep("observed",length(d$y))))
pd = position_dodge(width=0.4)
ggplot(data=freq,aes(x=RT,group=type,fill=type))+
geom_histogram(aes(y=after_stat(density),alpha=type),position=pd,binwidth=1)+
scale_fill_manual(values=c("darkblue","chartreuse4"))+
scale_alpha_manual(values=c(0.7,0.6))+
theme(text=element_text(size=ts,color=tc),axis.text.y=element_blank())+
scale_x_continuous(limits=c(-1,40))
_______
Esempio di test del rischio di falsi positivi sotto H0 per la variabile x, a causa dell’overdispersion.
set.seed(1)
N = 50
niter = 2000
p = rep(NA,niter)
for(i in 1:niter){
tryCatch({
x = rnorm(N,0,1)
z = rnorm(N,0,1)
w = rnorm(N,0,1)
y = rgamma(N,shape=exp(1 + x*0 + z + w)) # l'effetto di x non c'è, ma ci sono forti effetti di z e w...
d = data.frame(y,x) # ...z e w, però, non sono stati misurati, quindi non possono essere inseriti nel modello
fitG = glm(y ~ x, data=d, family=Gamma(link="log"))
su = summary(fitG)
p[i] = su$coefficients["x","Pr(>|t|)"]
},error=function(e){})
}
round(mean(p<0.05,na.rm=T),3)## [1] 0.097