# %% %%%%%%%%%%%%%%%%%%%%%%%%%%% # Simulazione Monte Carlo: # distribuzione del Chi-quadrato # %% %%%%%%%%%%%%%%%%%%%%%%%%%%% # # # %% La somma di n variabili casuali indipendenti Normali standard (z) elevatea al quadrato # %% segue una distribuzione Chi-quadrato con gdl pari agli elementi indipendenti che costituiscono # %% la sommatoria meno eventuali "vincoli" (evidenti o meno) di dipendenza tra gli addendi. # # SIMULIAMO: # # Creiamo 50,000 campioni causali Chi-quadrato Chisq_1 <-c() R=0 k1=6 while(R<=50000){ R=R+1 Z1<-rnorm(k1) # k1 variabili random "z" Chisq_1 <-c(Chisq_1,sum(Z1^2)) # Chisq(df=k1) = somma k1 variabili random "z" al quadrato cat("i=",R," \n") } # # par(mfrow=c(1,2)) # due plot per riga (1 riga e 2 colonne) hist(Chisq_1,prob=TRUE,border="white") # # %% La somma di k1=6 valori z indipendenti elevati al quadrato segue una distribuzione # %% Chi-quadrato, vincolata tra (0,+Inf), con Media = k1 e Varianza = 2k1. # plot(density(Chisq_1),col="blue") x=seq(1,80,by=0.01) curve(dchisq(x,df=k1),add=TRUE,col="red") abline(v=k1,col="orange",lwd=2,lty="dashed") abline(v=mean(Chisq_1),col="orange",lwd=2,lty="dashed") mean(Chisq_1) # = k1 var(Chisq_1) # = 2*k1 # # # # %% %%%%%%%%%%%%%%%%%%%%%%%%%%% # Simulazione Monte Carlo: # Differenza tra Chi-quadrato # %% %%%%%%%%%%%%%%%%%%%%%%%%%%% # # # %% La differenza tra due variabili indipendenti con distribuzione Chi-quadrato # %% non segue una distribuzione Chi-quadrato. La differenza può assumere valori # %% nell'intervallo (-Inf, +Inf), a differenza della distribuzione chi-quadrato standard, # %% che ha come limite inferiore lo zero. Per gradi di libertà "k" elevati, # %% la distribuzione della differenza può essere approssimata da una distribuzione normale # %% con media = k2-k1 e varianza = 2(k2 + 2k1). # # SIMULIAMO: # # Creiamo 50,000 campioni causali Chisq1 e Chisq2 INDIPENDENTI Chisq_1 <-c() Chisq_2 <-c() R=0 k1=60 k2=80 while(R<=50000){ R=R+1 Z1<-rnorm(k1) # k1 variabili random "z" Z2<-rnorm(k2) # k2 variabili random "z" Chisq_1 <-c(Chisq_1,sum(Z1^2)) # Chisq(df=k1) = somma k1 variabili random "z" al quadrato Chisq_2 <-c(Chisq_2,sum(Z2^2)) # Chisq(df=k2) = somma k2 variabili random "z" al quadrato cat("i=",R," \n") } # # par(mfrow=c(1,2)) # due plot per riga (1 riga e 2 colonne) hist(Chisq_2-Chisq_1,prob=TRUE,border="white") # # %% Nel caso di Chisq_1 e Chisq_2 indipendenti, con gradi di libertà elevati, avremo che: # %% Approsimazione Normale con Valore atteso[Chisq_2 - Chisq1] = (k2-k1), e # %% Varianza[Chisq_2 - Chisq1]= Var(Chisq_2) + Var(Chisq_1) = 2k2 + 2k1 = 2(k2 + 2k1) # plot(density(Chisq_2-Chisq_1),col="blue") x=seq(1,80,by=0.01) curve(dnorm(x,mean=(k2-k1),sd=(2*(k2+k1))^0.5),add=TRUE,col="red") abline(v=k2-k1,col="orange",lwd=2,lty="dashed") mean(Chisq_2-Chisq_1) # = (k2-k1) var(Chisq_2-Chisq_1) # = 2(k2 + 2k1) # # # %% Idea Sbaglita: Nella pratica la differenze tra due Chi-quadrato indipendenti è raramente # %% una variabile Chi-quadrato, potendo assumere anche valori negativi ed essendo confinata tra -Inf e + Inf. # %% La differenza tra due variabili Chi-quadrato segue una distribuzione Chi-quadrato se e solo se # %% le variabili sono dipendenti in un modo specifico, che potremo definire "nested", # %% un requisito NECESSARIO ANCHE PER IL TEST "LRT". Vediamo di chiarire... # # SIMULIAMO: # # Creiamo 50,000 campioni causali Chisq1 e Chisq2 DIPENDENTI Chisq_1 <-c() Chisq_2 <-c() R=0 k1=6 # In questo caso non servono gdl elevati per avere una distribuzione nota k2=8 # In questo caso non servono gdl elevati per avere una distribuzione nota while(R<=50000){ R=R+1 Z2<-rnorm(k2) # k2 variabili random "z" Chisq_2 <-c(Chisq_2,sum(Z2^2)) # Chisq(df=k2) = somma k2 variabili random "z" al quadrato # %% Chisq_1 <-c(Chisq_1,sum(Z2[1:k1]^2)) # Chisq(df=k1) = somma k1 elementi vettore Z2 al quadrato # %% La somma Chisq_1 è "nidificata" nella somma Chisq_2, infatti, mediante un # %% vettore v di contrasti c(1,1,1,1,1,1,0,0) otteniamo lo stesso risultato! v=c(rep(1,k1),rep(0,k2-k1)) Chisq_1 <-c(Chisq_1,Z2^2%*%v) # otteniamo lo stesso risultato! cat("i=",R," \n") } # # # par(mfrow=c(1,2)) hist(Chisq_2-Chisq_1,prob=TRUE,border="white") # # %% Nel caso di Chisq_1 e Chisq_2 DIPENDENTI ("nested", so to speak), avremo che: # %% Distribuzione Chi-quadrato, con # %% Valore atteso[Chisq_2 - Chisq1] = gdl = (k2-k1), # %% Varianza[Chisq_2 - Chisq1] = 2gdl = 2(k2-k1) # plot(density(Chisq_2-Chisq_1),col="blue") x=seq(1,80,by=0.01) curve(dchisq(x,df=(k2-k1)),add=TRUE,col="red") abline(v=(k2-k1),col="blue",lwd=2) mean(Chisq_2-Chisq_1) # = gdl = (k2-k1) var(Chisq_2-Chisq_1) # = 2(k2-k1) # # # # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # Simulazione Monte Carlo: # Distrubuzione chi-quadrato LRT # modello ridotto (H0) vs completo (H1) # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # Modello teorico esatto per i dati: # # x1<-seq(-3,3,by=0.01) x2<-x1+runif(length(x1)) #variabile aggiuntiva (H1) x3<-x2+runif(length(x1)) #variabile aggiuntiva (H1) x4<-x3+runif(length(x1)) #variabile aggiuntiva (H1) x5<-x4+runif(length(x1)) #variabile aggiuntiva (H1) ssize=length(x1) # # Modello che verifica H0: B2 = B3 = B4 = B5 = 0 "Nested" rispetto a H1 # p.0<-exp(0.4 + 3*x1 + 0*x2 + 0*x3 + 0*x4 + 0*x5)/(1+exp(0.4 + 3*x1 + 0*x2 + 0*x3 + 0*x4 + 0*x5)) #(H0): x2-3-4-5 effetto nullo # # # Generazione di un vettore random y con distribuzione dicotomica Bernoulli: # "p.0" prevista dal modello teorico su base x1 (x2-3-4-5 non significative) y<-rbinom(n=ssize,size=1,p=p.0) plot(p.0) # # # Creiamo R=50,000 campioni causali y R=0 Y=c() while(R<=50000){ R=R+1 Y<-rbind(Y,rbinom(n=ssize,size=1,p=p.0)) cat("i=",R," \n") } # # # 50,000 analisi glm per i campioni causali y LRT=c() for(i in 1:dim(Y)[1]){ mod<-glm(Y[i,]~x1+x2+x3+x4+x5,family=binomial) #(H1): x2,x3,x4, e x5 "inserite" residual.deviance<- 2*(0-c(logLik(mod))) # 2*[delta log-likelihood "saturo vs. modello stimato"]. H0: # fit del modello simato non sig. diverso dal modello saturo mod.null<-glm(Y[i,]~x1,family=binomial) #(H0): SOLO x1 "inserita" null.deviance<- 2*(0-logLik(mod.null)) # 2*[delta log-likelihood "saturo vs. modello nullo"]. H0: # fit del modello simato non sig. diverso dal modello saturo LRT=c(LRT,null.deviance-residual.deviance) cat("sample = ",i,"\n") } summary(LRT) plot(density(LRT),xlim=c(0,max(LRT)),col="red") #hist(LRT,prob=TRUE) abline(v=4) curve(dchisq(x,df=4),from=0, to=25,add=TRUE,col="blue") # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ANALISI DATI SUICIDIO # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # Regressione logistica # VD = variabili suicidarie, PENSIERO SUICIDARIO-PIANIFICAZIONE SUICIDARIA-TENTATO SUICIDIO # Predittori di interesse: dipendenza da sostanze per ciascuna sostanza # 2008 # IRSEX = SEX (1-M; 2-F) # MHSUITHK = PENSIERO SUICIDARIO # MHSUIPLN = PIANIFICAZIONE SUICIDARIA # MHSUITRY = Tentamen # NDSSDNSP = Nicotine Dependence Syndrome Scale (NDSS) # DEPNDALC = ALCOHOL # DEPNDCOC = COCAINE # DEPNDHER = HEROINE # DEPNDANL = PAIN RELIEVERS # DEPNDSED = SEDATIVES # DEPNDSTM = STIMULANTS # DEPNDHAL = HALLUCINOGENS # DEPNDINH = INHALANTS # DEPNDMRJ = MARIJUANA # DEPNDTRN = TRANQUILIZERS ##VARS<-c("IRSEX","MHSUITHK","MHSUIPLN","MHSUITRY","NDSSDNSP","DEPNDALC","DEPNDCOC","DEPNDHER","DEPNDANL","DEPNDSED","DEPNDSTM","DEPNDHAL","DEPNDINH","DEPNDMRJ","DEPNDTRN") ##load("dataset08.RData") ##dataset08_reduced<-dataset08[VARS] ##save("dataset08_reduced",file="dataset08.RData") ## Ricaviamo la directory locale della sessione corrente R ## mettiamo qui dentro il file RData: getwd() ## [1] "/Users/michelegrassi" ## carichiamo i dati: load("dataset08.RData") ## Convertiamo SEX in fattore: 1-Maschio; 2-Femmina dataset08_reduced$IRSEX<-factor(dataset08_reduced$IRSEX,labels=names(attributes(dataset08_reduced$IRSEX)$labels)) summary(dataset08_reduced) ## Convertiamo le dipendenze in fattori "1=Si"/"0=no" (non necessario se già in formato Dummy numerica 0-1) ## VARS<-c("NDSSDNSP","DEPNDALC","DEPNDCOC","DEPNDHER","DEPNDANL","DEPNDSED","DEPNDSTM","DEPNDHAL","DEPNDINH","DEPNDMRJ","DEPNDTRN") ## dataset08_reduced[VARS] <- lapply(dataset08_reduced[VARS] , factor) summary(dataset08_reduced) ## ATTENZIONE: Convertiamo "MISSING = -9" in SPSS --> NA in R dataset08_reduced$MHSUITHK<-ifelse(dataset08_reduced$MHSUITHK=="-9",NA,dataset08_reduced$MHSUITHK) dataset08_reduced$MHSUIPLN<-ifelse(dataset08_reduced$MHSUIPLN=="-9",NA,dataset08_reduced$MHSUIPLN) dataset08_reduced$MHSUITRY<-ifelse(dataset08_reduced$MHSUITRY=="-9",NA,dataset08_reduced$MHSUITRY) summary(dataset08_reduced) ## ATTENZIONE: Rimuoviamo i missing data nelle variabili di "SUICIDAL OUTCOME" [MHSUITHK, MHSUIPLN, MHSUITRY] dataset08_reduced=na.omit(dataset08_reduced) summary(dataset08_reduced) # Regressione logistica: PENSIERO SUICIDARIO Mod.PS_08<-glm(MHSUITHK ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + DEPNDHAL + DEPNDINH + DEPNDMRJ + DEPNDTRN, data = dataset08_reduced, family=binomial) summary(Mod.PS_08) # Regressione logistica: PIANIFICAZIONE SUICIDARIA Mod.PianS_08<-glm(MHSUIPLN ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + DEPNDHAL + DEPNDINH + DEPNDMRJ + DEPNDTRN, data = dataset08_reduced, family=binomial) summary(Mod.PianS_08) # Regressione logistica: TENTATO SUICIDIO Mod.Tent_08<-glm(MHSUITRY ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + DEPNDHAL + DEPNDINH + DEPNDMRJ + DEPNDTRN, data = dataset08_reduced, family=binomial) summary(Mod.Tent_08) # %% Call: # %% glm(formula = MHSUITRY ~ IRSEX + NDSSDNSP + DEPNDALC + DEPNDCOC + # %% DEPNDHER + DEPNDANL + DEPNDSED + DEPNDSTM + DEPNDHAL + DEPNDINH + # %% DEPNDMRJ + DEPNDTRN, family = binomial, data = dataset08_reduced) # %% # %% Coefficients: # %% Estimate Std. Error z value Pr(>|z|) # %% (Intercept) -5.34664 0.10121 -52.829 < 2e-16 *** # %% IRSEX 0.49910 0.11260 4.432 9.32e-06 *** # %% NDSSDNSP 0.69949 0.12646 5.531 3.18e-08 *** # %% DEPNDALC 1.43131 0.13866 10.322 < 2e-16 *** # %% DEPNDCOC 0.92608 0.27006 3.429 0.000606 *** # %% DEPNDHER 0.66303 0.55403 1.197 0.231412 # %% DEPNDANL 0.97931 0.25954 3.773 0.000161 *** # %% DEPNDSED 0.58977 0.83398 0.707 0.479457 # %% DEPNDSTM 1.23972 0.42353 2.927 0.003422 ** # %% DEPNDHAL -0.46863 0.87252 -0.537 0.591195 # %% DEPNDINH -9.76738 187.19446 -0.052 0.958387 ## SPARSE DATA!! # %% DEPNDMRJ 0.67954 0.20777 3.271 0.001073 ** # %% DEPNDTRN -0.08888 0.57110 -0.156 0.876323 # %% --- # %% Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # %% # %% (Dispersion parameter for binomial family taken to be 1) # %% # %% Null deviance: 4012.5 on 37356 degrees of freedom # %% Residual deviance: 3735.1 on 37344 degrees of freedom # %% AIC: 3761.1 # %% # %% Number of Fisher Scoring iterations: 12 # %% # %% "SPARSE data" per la variabile indipendente "DEPNDINH"! (coefficiente/errore standard molto grandi) # %% Verifichiamo cosa succede con il nostro algoritmo Fisher scoring (...time consuming) # %% Costruiamo il vettore di osservazioni binarie "y" = Tentativo di suicidio y<-c(dataset08_reduced$MHSUITRY) # %% Convertiamo il fattore IRSEX a variabile Dummy numerica 0=maschio, 1=femmina Dummy_IRSEX<-ifelse(dataset08_reduced$IRSEX=="Male",0,1) # %% Costruiamo la matrice del modello X X<-dataset08_reduced[,-c(1:4)] X<-cbind(rep(1,dim(X)[1]),Dummy_IRSEX,X) is.matrix(X) # %% Rendiamola riconoscibile da R come "matrice" di dati: # %% in caso contrario, le operazioni matriciali non saranno possibili! X<-as.matrix(X) is.matrix(X) # %% PREAMBOLI j=dim(X)[2]-1 i=0 phit0<-matrix(rep(0,j+1)) # data set vuoto per i risultati: FisherScoring<-c() # Per macOS: Innalzare il limite di RAM utilizzabile... mem.maxVSize() mem.maxVSize(vsize = Inf) 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 dim(FisherScoring) # 31 cicli poi singolarità della matrice Hessiana --> Sparse Data! ## [1] 13 31 ## Al 12mo ciclo le stime coincidono con quelle glm(), che si ferma prima della singolarità della matrice Hessiana FisherScoring[,12] - Mod.Tent_08$coeff round(FisherScoring[,12] - Mod.Tent_08$coeff,6) ## arrotondamento sesto decimale ## (Intercept) IRSEXFemale NDSSDNSP DEPNDALC DEPNDCOC DEPNDHER DEPNDANL DEPNDSED DEPNDSTM DEPNDHAL DEPNDINH ## 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.354319 ## DEPNDMRJ DEPNDTRN ## 0.000000 0.000000 ## Nelle succesive iterazioni Fisher Scoring tutti i coefficienti rimangono stabili nella stima MLE tranne "DEPNDINH": ## per questa variabile il coefficiente B (riga 11) tende a diminuire verso -Inf... FisherScoring[,22:31] # %% Result_iter Result_iter Result_iter Result_iter Result_iter Result_iter Result_iter Result_iter Result_iter Result_iter # %% [1,] -5.34663638 -5.34663638 -5.34663638 -5.34663638 -5.34663638 -5.34663638 -5.34663638 -5.34663638 -5.34663638 -5.34663638 # %% [2,] 0.49909776 0.49909776 0.49909776 0.49909776 0.49909776 0.49909776 0.49909776 0.49909776 0.49909776 0.49909776 # %% [3,] 0.69948707 0.69948707 0.69948707 0.69948707 0.69948707 0.69948707 0.69948707 0.69948707 0.69948707 0.69948707 # %% [4,] 1.43131058 1.43131058 1.43131058 1.43131058 1.43131058 1.43131058 1.43131058 1.43131058 1.43131058 1.43131058 # %% [5,] 0.92608326 0.92608326 0.92608326 0.92608326 0.92608326 0.92608326 0.92608326 0.92608326 0.92608326 0.92608326 # %% [6,] 0.66302850 0.66302850 0.66302850 0.66302850 0.66302850 0.66302850 0.66302850 0.66302850 0.66302850 0.66302850 # %% [7,] 0.97930832 0.97930832 0.97930832 0.97930832 0.97930832 0.97930832 0.97930832 0.97930832 0.97930832 0.97930832 # %% [8,] 0.58977192 0.58977192 0.58977192 0.58977192 0.58977192 0.58977192 0.58977192 0.58977192 0.58977192 0.58977192 # %% [9,] 1.23971843 1.23971843 1.23971843 1.23971843 1.23971843 1.23971843 1.23971843 1.23971843 1.23971843 1.23971843 # %% [10,] -0.46863362 -0.46863362 -0.46863362 -0.46863362 -0.46863362 -0.46863362 -0.46863362 -0.46863362 -0.46863362 -0.46863362 # %% [11,] -19.41306910 -20.41306910 -21.41306910 -22.41306910 -23.41306910 -24.41306910 -25.41306910 -26.41306910 -27.41306910 -28.41306910 # %% [12,] 0.67953789 0.67953789 0.67953789 0.67953789 0.67953789 0.67953789 0.67953789 0.67953789 0.67953789 0.67953789 # %% [13,] -0.08888073 -0.08888073 -0.08888073 -0.08888073 -0.08888073 -0.08888073 -0.08888073 -0.08888073 -0.08888073 -0.08888073 ## Tabella: Tentativo di Suicidio (y) 0-1 (in colonna) per "DEPNDINH" 0-1 (in riga): table(dataset08_reduced$DEPNDINH,y) ## y ## 0 1 ## 0 36995 355 ## 1 7 0 ## ## Presente una cella con valore di frequenza 0 ==> SPARSE DATA! ## ODDS RATIO: [ p(y=1|DEPNDINH=1)/p(y=0|DEPNDINH=1) ]/[ p(y=1|DEPNDINH=0)/p(y=0|DEPNDINH=0) ] (0/7)/(355/36995) ## [1] 0 ## Coefficiente di regressione logistica B per variabile "DEPNDINH" = log(ODDS RATIO) log(0) ## [1] -Inf