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.7334

Se 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.001

Un 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: 5
fitP = 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")