#################### Sparse Data problem ################### Sparse Data problem ################# Sparse Data problem ###############: Example 1 - x categoriale y<-c(rep(1,8),rep(0,2),rep(1,0),rep(0,10)) x<-factor(c(rep("A",10),rep("B",10))) table(x,y) n<-length(y) # Set A as base level X<-matrix(c(rep(1,n),ifelse(x=="A",0,1)),ncol=2,byrow=FALSE) X prop.table<-table(x,y)/10 Odds.Ratio<-(prop.table[2,2]/prop.table[2,1])/(prop.table[1,2]/prop.table[1,1]) Odds.Ratio log(Odds.Ratio) ###### Fisher Scoring j=1 i=0 phit0<-matrix(c(0,0)) # data set vuoto per i risultati: FisherScoring<-c() repeat{ i=i+1 # prima iterazione (la prima sarà nulla, quindi registreremo i-1) B<-phit0[1:(j+1)] p<-exp(X%*%B)/(1+exp(X%*%B)) u<-t(X)%*%(matrix(y)-p) W<-diag( # Matrice diagonale -1*p(1-p) c(-p*(1-p)) ) H<-t(X)%*%W%*%X phit1<- phit0 - solve(H)%*%u #registriamo i risultati nel data set FisherScoring: Result_iter<-c(phit1[1:(j+1)]) FisherScoring<-cbind(FisherScoring, Result_iter) #condizione logica per interrompere il ciclo: if(sum(abs(phit1 - phit0) < .0001)==(j+1)){ rownames(FisherScoring)<-c(paste("b",seq(1,j+1))) colnames(FisherScoring)<-seq(0,(i-1)) print(FisherScoring) break } phit0<-phit1 } FisherScoring summary(glm(y~x,family=binomial(link=logit))) ###############: Example 2 - x continua x = c(0,1,2,3,4,5,6,7) y = c(0,0,0,0,1,1,1,1) n<-length(y) X<-matrix(c(rep(1,n),x),ncol=2,byrow=FALSE) X ###### Fisher Scoring j=1 i=0 phit0<-matrix(c(0,0)) # data set vuoto per i risultati: FisherScoring<-c() repeat{ i=i+1 # prima iterazione (la prima sarà nulla, quindi registreremo i-1) B<-phit0[1:(j+1)] p<-exp(X%*%B)/(1+exp(X%*%B)) u<-t(X)%*%(matrix(y)-p) W<-diag( # Matrice diagonale -1*p(1-p) c(-p*(1-p)) ) H<-t(X)%*%W%*%X phit1<- phit0 - solve(H)%*%u #registriamo i risultati nel data set FisherScoring: Result_iter<-c(phit1[1:(j+1)]) FisherScoring<-cbind(FisherScoring, Result_iter) #condizione logica per interrompere il ciclo: if(sum(abs(phit1 - phit0) < .0001)==(j+1)){ rownames(FisherScoring)<-c(paste("b",seq(1,j+1))) colnames(FisherScoring)<-seq(0,(i-1)) print(FisherScoring) break } phit0<-phit1 } summary(glm(y~x,family=binomial(link=logit))) # Plot dei dati "sparsi": plot(x,y,pch=19,yaxt="n",bty="n") axis(2,at=c(0,0.5,1)) #probabilità previste dal modello "stimato" x.sim<-seq(0,7,by=0.001) x.sim[1:20] p1<-exp( -2.6666667 + x.sim* 0.7619048 )/( 1 + exp( -2.6666667 + x.sim* 0.7619048 ) ) p5<-exp( -17.730044 + x.sim*5.065727 )/( 1 + exp( -17.730044 + x.sim*5.065727 ) ) p25<-exp( -165.22410 + x.sim*47.20689 )/( 1 + exp( -165.22410 + x.sim*47.20689 ) ) # Notiamo che il rapporto tra alpha e beta nel modello stimato è sempre = -3.5: FisherScoring[1,c(1,5,25)]/FisherScoring[2,c(1,5,25)] # scegliamo quindi valori arbitrariamente grandi (approx Inf) pInf<-exp( -350 + x.sim*100 )/( 1 + exp( -350 + x.sim*100 ) ) # # # points(x.sim, p1,pch=19,type="l",col="green",lwd=2) points(x.sim, p5,pch=19,type="l",col="orange",lwd=2) points(x.sim, p25,pch=19,type="l",col="red",lwd=2) points(x.sim, pInf,pch=19,type="l",col="blue",lwd=2) ###############: Example 3 più variabili x continue (x casa) x1 = c(0.5, 2, 10, 6, 8, 1, 8, 4, 0, 7, 5, 1.1) x2 = c(7, 9.5, 6, 0, 8.1, 9, 10, 0, 8, 5, 3, 4) y = c(0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0) summary(glm(y~x1+x2,family=binomial(link=logit))) #install.packages("rgl") #install.packages("magick") mycolors <- ifelse(y==1,'red', 'blue') # Tahoe workaround to load rgl library: # Start a new session options(rgl.useNULL = TRUE) library(rgl) options(rgl.printRglwidget = TRUE) # Plot plot3d( x=x1, y=x2, z=y, col = mycolors, type = 's', radius = .25, xlab="X1", ylab="X2", zlab="Y") rglwidget() # To save to a html file htmlwidgets::saveWidget(rglwidget(width = 520, height = 520), file = "~/3dscatter.html", libdir = "libs", selfcontained = FALSE ) #################### Sparse Data problem ################### Sparse Data problem ################# Sparse Data problem x<-c(8,14,-7,6,5,6,-5,1,0,-17) y<-c(1,1,0,0,1,0,1,0,0,0) n<-length(y) X<-matrix(c(rep(1,n),x),ncol=2,byrow=FALSE) ############# GLM summary(glm(y~x,family=binomial(link=logit))) ###### Fisher Scoring j=1 i=0 phit0<-matrix(c(0,0)) n<-length(y) FisherScoring<-c() repeat{ i=i+1 # prima iterazione (la prima sarà nulla, quindi registreremo i-1) B<-phit0[1:(j+1)] p<-exp(X%*%B)/(1+exp(X%*%B)) u<-t(X)%*%(matrix(y)-p) W<-diag( c(-p*(1-p)) ) H<-t(X)%*%W%*%X phit1<- phit0 - solve(H)%*%u #registriamo i risultati nel data set FisherScoring: Result_iter<-c(phit1[1:(j+1)]) FisherScoring<-cbind(FisherScoring, Result_iter) #condizione logica per interrompere il ciclo: if(all(abs(phit1 - phit0) < .00001)){ rownames(FisherScoring)<-c(paste("b",seq(0,j))) colnames(FisherScoring)<-seq(0,(i-1)) print(FisherScoring) break } phit0<-phit1 } H # Test per B0 (j=1) var_B0 = -1*solve(H)[1,1] sqrt(var_B0) z_B0=(-0.7227534-0)/sqrt(var_B0) z_B0 pnorm(abs(z_B0),lower.tail=FALSE)*2 # Test per B1 (j=2) var_B1 = -1*solve(H)[2,2] sqrt(var_B1) z_B1=(0.1396 -0)/sqrt(var_B1) z_B1 pnorm(abs(z_B1),lower.tail=FALSE)*2 ## Intervalli di fiducia al 95% per i due coefficienti MLE -0.7228 + c(-1,+1)*1.96*sqrt(var_B0) ## [1] -2.3169131 0.8713131 0.1396 + c(-1,+1)*1.96*sqrt(var_B1) ## [1] -0.08672356 0.36592356 confint.default(glm(y~x,family=binomial(link=logit))) ## 2.5 % 97.5 % ## (Intercept) -2.31678189 0.871275 ## x -0.08667973 0.365936 ## Delta Method per il confronto tra coefficienti di regressione Logit x1<-c(8,14,-7,6,5,6,-5,1,0,-17) x2<-c(2,14,6,4,4,4,-7,2,10,-19) y<-c(1,1,0,0,1,0,1,0,0,0) ############# GLM mod<-glm(y~x1+x2,family=binomial(link=logit)) phi<-coefficients(mod) phi[2]-phi[3] v<-c(0,1,-1) # vettore di contrasti v%*%phi n_inv_H<-vcov(mod) Chisq<-t(v%*%phi)%*%solve(t(v)%*%n_inv_H%*%v)%*%(v%*%phi) Chisq # Valore della statistica di test, con 1 g.d.l. 1-pchisq(Chisq,df=1) # area della coda Chi-quadrato destra ## Delta Method per il rapporto tra coefficienti di regressione Logit mod$coefficients Ratio<-0.4801557/(-0.3801478) Ratio ## Calcoliamo il jacobiano della funzione h=phi1/phi2 h=expression(phi2/phi3) D(h,"phi2") D(h,"phi3") phi<-coefficients(mod) J=matrix(c(0,1/phi[3],-phi[2]/phi[3]^2)) ## Estrarre la matrice di varianza/covarianza delle stime MLE vcov(mod) ## Sandwich multiplication t(J)%*%vcov(mod)%*%J car::deltaMethod(mod,"x1/x2") # Intervallo di fiducia per i valori previsti (probabilità) x<-c(8,14,-7,6,5,6,-5,1,0,-17) y<-c(1,1,0,0,1,0,1,0,0,0) X<-matrix(c(rep(1,n),x),ncol=2,byrow=FALSE) mod<-glm(y~x,family=binomial(link=logit)) # Logit previsto per x=10 v=c(1,10) phi=coefficients(mod) t(v)%*%phi # Logit previsti per tutti i valori x Logit<-X%*%phi # Plot: valori previsti su scala lineare Logit: plot(x,Logit,pch=19,bty="l",xlim=c(-20,20),ylim=c(-10,10)) # modello lineare teorico points(x,Logit,pch=19,type="l",col="red",lty="dashed",lwd=2) # Limiti di fiducia al 95% per i valori lineari Logit = a + bx # Limiti di fiducia al 95% per i valori lineari Logit = a + bx # Limiti di fiducia al 95% per i valori lineari Logit = a + bx # inf 95% Inf95<-(v)%*%phi + 1.96*sqrt(t(v)%*%vcov(mod)%*%v) # sup 95% Sup95<-t(v)%*%phi - 1.96*sqrt(t(v)%*%vcov(mod)%*%v) # disegniamo i due punti estremi: points(x=c(10,10),y=c(Inf95,Sup95),pch=19col="green") ## per tutti i valori x ## formula matriciale: ## diag(X%*%vcov(mod)%*%t(X)) Var_Logit<-c() for(i in 1:length(x)){ v<-c(1,x[i]) Var_Logit <-c(Var_Logit,t(v)%*%vcov(mod)%*%v) } Inf95.logit<- Logit -1.96*sqrt(Var_Logit) Sup95.logit<- Logit +1.96*sqrt(Var_Logit) points(x, Inf95.logit,pch=19,col="green",lty="dotted") points(x, Sup95.logit,pch=19,col="green",lty="dotted") # soluzione algebrica, sfruttando la bilinearità della funzione di covarianza # Equazione 66 vcov(mod)[1,1] + x^2*vcov(mod)[2,2] + 2*x*vcov(mod)[1,2] # Var_Logit # calcolato in precedenza # diag(X%*%vcov(mod)%*%t(X)) # formula matriciale metodo Delta ## Plot: Valori previsti (probabilità) su scala non lineare logistica: plot(x,y,pch=19,bty="l",xlim=c(-20,20),ylim=c(0,1)) ## Antilogit: dal Logit alla probabilità p<-exp(Logit)/(1+exp(Logit)) points(x,p,pch=19,col="red") ## Costruiamo la logistica teorica logistic<-function(x){exp(phi[1]+phi[2]*x)/(1+exp(phi[1]+phi[2]*x))} curve(logistic,from=-20,to=20,col="blue",lty="dashed",add=TRUE) ## Limiti di fiducia al 95% per i valori non lineari "p" ## Limiti di fiducia al 95% per i valori non lineari "p" ## Limiti di fiducia al 95% per i valori non lineari "p" ## Applicazione formula inversa Logit (antilogit) ## inf 95% Inf95.antilogit<- exp(Inf95.logit)/(1+exp(Inf95.logit)) ## sup 95% Sup95.antilogit<- exp(Sup95.logit)/(1+exp(Sup95.logit)) points(x, Inf95.antilogit,pch=19,col="green",lty="dotted") points(x, Sup95.antilogit,pch=19,col="green",lty="dotted") ## Interpolazione continua funzione spline()/approx() points(spline(x, Inf95.antilogit),type="l",col="green",lty="dotted") points(spline(x, Sup95.antilogit),type="l",col="green",lty="dotted") ## Likelihood ratio test : LRT ## Likelihood ratio test : LRT ## Likelihood ratio test : LRT # dati: x1<-c(8,14,-7,6,5,6,-5,1,0,-17) x2<-c(2,14,6,4,4,4,-7,2,10,-19) y<-c(1,1,0,0,1,0,1,0,0,0) X<- matrix(c(rep(1,length(y)),x1,x2),byrow=FALSE,ncol=3) # modello con solo l’intercetta mod.H0<-glm(y~1,binomial(link = "logit")) summary(mod.H0) # Interpretazione: exp(Intercetta) = p/(1-p) per il campione intero exp(mod.H0$coeff[1]) mean(y)/(1-mean(y)) # modello con i parametri delle covariate diversi da 0 mod.H1<-glm(y~X[,2]+X[,3],binomial(link = "logit")) phi<-coefficients(mod.H1) summary(mod.H1) # Modello H1 # Deviance: -2*argmax(Log_Lik) Deviance.H1<- -2*(sum(y*X%*%phi-log(1+exp(X%*%phi)))) Deviance.H1 mod.H1$deviance # gradi di libertà n=10 # campione = length(y) k=2 # variabili x Df.H1<-n-k-1 Df.H1 # Modello H0 # Deviance: -2*argmax(Log_Lik) phi_H0<-c(coefficients(mod.H0),0,0) phi_H0 Deviance.H0<- -2*(sum(y*X%*%phi_H0-log(1+exp(X%*%phi_H0)))) Deviance.H0 mod.H0$deviance # gradi di libertà n=10 # campione = length(y) k=0 # solo intercetta Df.H0<-n-k-1 Df.H0 summary(mod.H1) # Likelihood ratio test (LRT): # distribuzione: Chi-quadrato, df = Df_H0 – Df.H1 # per il confronto tra modello nullo e alternativo LRT<- Deviance.H0 - Deviance.H1 LRT Df.H0-Df.H1 1-pchisq(LRT,df=2) # Confronto tra i due modelli stimati con glm() e anova() anova(mod.H0,mod.H1,test="Chisq") ## Modello nullo e saturo ## Approfondimento mod.H0<-glm(y~1,binomial(link = "logit")) # Nullo summary(mod.H0) mod.sat<-glm(y~factor(1:length(y)),binomial(link = "logit")) # saturo summary(mod.sat) mod.H1<-glm(y~X[,2]+X[,3],binomial(link = "logit")) # "buono"? summary(mod.H1) ## Deviance Chi-square tests of Goodness of Fit: logLik(mod.sat) ## df = 0 logLik(mod.H0) 1-pchisq(-2*c(logLik(mod.H0)),df=n-1) logLik(mod.H1) 1-pchisq(-2*c(logLik(mod.H1)),df=n-3) ##LRT LRT<- (-2*c(logLik(mod.H0))) - (-2*c(logLik(mod.H1))) LRT Df=9-7 1-pchisq(LRT,df=Df)