Annexe 3 : Scripts R des modèles d’association utilisés dans l’article
Plan
Haut de pageTexte intégral
1Voici les scripts permettant d’estimer les modèles depuis le logiciel R.
Installation logiciel et environnement
2Pour installer R : https://cran.r-project.org/
3Pour installer RStudio (optionel mais conseillé) : https://www.rstudio.com/products/rstudio/
Code R pour reproduire les modèles 1 à 13 (excepté les modèles 10 et 11)
4# Il faut préalablement installer les packages gnm et logmult
5library(gnm)
6library(logmult)
7# Fonction qui calcule le BIC basé sur la statistique du L2 --------------
8model.summary <- function(obj){
9aic <- obj$deviance - obj$df.residual * 2 # AIC(L2)
10bic <- obj$deviance - obj$df.residual * log(sum(obj$y)) # BIC(L2)
11# Index of Dissimilarity (%)
12delta <- 100 * sum(abs(obj$y - obj$fitted.values)) / (2 * sum(obj$y))
13p <- pchisq(obj$deviance, obj$df.residual, lower.tail = FALSE)
14# p <- ifelse(p < 0.001,"< 0.001", p)
15result <- matrix(0,1,6)
16rownames(result) <- ""
17colnames(result) <- c("Df","L2","p","AIC(L2)","BIC(L2)","Delta")
18# exp.freq=obj$fitted.values
19result[1,] <- c(obj$df.residual, obj$deviance, p, aic, bic, delta)
20print("Model Summary:")
21return(result)
22}
23# Données -----------------------------------------------------------------
24Freq<-c(1, 4, 5, 2,
259, 28, 15, 25,
2615, 20, 13, 6,
2720, 19, 26, 5,
2814, 13, 9, 1,
2918, 21, 7, 2,
3016, 12, 4, 7,
3153, 23, 18, 11,
3243, 24, 47, 7,
33136, 84, 160, 152,
34187, 65, 109, 133,
35228, 110, 237, 144,
36404, 111, 164, 51,
37307, 106, 162, 161,
3887, 37, 30, 52,
3980, 52, 49, 54,
4022, 18, 31, 5,
4156, 44, 42, 82,
4252, 42, 42, 45,
43109, 92, 76, 58,
44239, 101, 58, 62,
45143, 86, 87, 100,
4634, 24, 18, 42,
4710, 7, 3, 10)
48tab <- array(Freq, dim=c(4,8,3),
49dimnames = list("parti" = c("gauche", "centre.droit", "liberaux", "populiste"),
50"class" = c("dirigeants", "indep", "expert.tech", "manager", "sociocult", "employe", "ouvrier", "trav.nq"),
51"year" = c("1975", "1995", "2011")))
52tab_df <- as.data.frame(as.table(tab))
53tab_df
54# Modèle 1 - Indépendance complète {R, C, L} ------------------------------------------------------
55m1 <- glm(Freq ~ class + parti + year, data = tab_df, family = poisson)
56model.summary(m1)
57# Modèle 2 - Toutes les interactions de second-ordre {RC, RL, CL} -----------------------------------
58m2 <- glm(Freq ~ class * year + parti * year + parti * class, data = tab_df, family = poisson)
59model.summary(m2)
60# Modèle 3 - RC(0)L homogène: indépendance conditionnelle {RL, CL} -----------------------------
61m3 <- glm(Freq ~ class * year + parti * year, data = tab_df, family = poisson)
62model.summary(m3)
63# Modèle 4 - RC(1)L homogène ------------------------------------------------------
64m4 <- rcL(tab, nd = 1, layer.effect = "none", weighting = "none", se = "none")
65summary(m4)
66pchisq(m4$deviance, df = m4$df.residual, lower.tail = FALSE)
67# Modèle 5 - RC(1)L hétérogène ----------------------------------------------------
68m5 <- rcL(tab, nd = 1, layer.effect = "heterogeneous", weighting = "none", se = "none")
69summary(m5)
70pchisq(m5$deviance, df = m5$df.residual, lower.tail = FALSE)
71# Modèle 6 - RC(2)L homogène ------------------------------------------------------
72m6 <- rcL(tab, nd = 2, layer.effect = "none", weighting = "none", se = "none")
73summary(m6)
74pchisq(m6$deviance, df = m6$df.residual, lower.tail = FALSE)
75# Modèle 7 - RC(2)L hétérogène ----------------------------------------------------
76m7 <- rcL(tab, nd = 2, layer.effect = "heterogeneous", weighting = "none", se = "none")
77summary(m7)
78pchisq(m7$deviance, df = m7$df.residual, lower.tail = FALSE)
79# Modèle 8 - RC(3)L homogène ------------------------------------------------------
80m8 <- rcL(tab, nd = 3, layer.effect = "none", weighting ="none", se = "none")
81summary(m8)
82pchisq(m8$deviance, df = m8$df.residual, lower.tail = FALSE)
83# Modèle 9 - RC(3)L hétérogène ----------------------------------------------------
84m9 <- rcL(tab, nd = 3, layer.effect = "heterogeneous", weighting = "none", se = "none")
85summary(m9)
86pchisq(m9$deviance, df = m9$df.residual, lower.tail = FALSE)
87# Modèle 10 - RC(2)-L homogène dans µi1, νj1, µi2 et νj2 -------------------------------------------------
88#
89# Ne converge pas avec gnm et logmult
90#
91# Modèle 11 - RC(2)-L homogène dans µi1, νj1, µi2, νj2 et ϕ2 ---------------------------------------------
92#
93# Ne converge pas avec gnm et logmult
94#
95# Modèle 12 - RC(2)-L homogène dans µi1, νj1, µi2, νj2 et ϕ2, avec tendance linéaire dans ϕ1 --------
96Parti<-gl(4,1,8*4*3)
97Clas<-gl(8,4,8*4*3)
98L<-gl(3,32,8*4*3)
99levels(L)<-factor(1:3)
100Rscore<-as.numeric(Clas)
101Cscore<-as.numeric(Parti)
102Year<-L
103levels(Year)<-c(0,20,36)
104nYear<-as.numeric(as.vector(Year))
105options(scipen=999)
106Clas.c1<-Clas
107levels(Clas.c1)<-factor(c(1,2,3,4,5,6,7,8))
108Clas.c1
109Clas.c2<-Clas
110levels(Clas.c2)<-factor(c(1,2,3,4,5,6,7,8))
111Clas.c2
112Parti.c1<-Parti
113levels(Parti.c1)<-factor(c(1,2,3,4))
114Parti.c1
115Parti.c2<-Parti
116levels(Parti.c2)<-factor(c(1,2,3,4))
117Parti.c2
118## IL FAUT PARFOIS RELANCER PLUSIEURS FOIS LA LIGNE DE COMMANDE POUR QUE LE MODELE m12 CONVERGE ##
119m12 <- gnm(Freq~Clas+Parti+L+Clas*L+Parti*L+Mult(1+nYear,Clas.c1,Parti.c1)+Mult(1,Clas.c2,Parti.c2),family=poisson,tolerance = 1e-12,iterStart=10,iterMax=50000)
120summary(m12)
121model.summary(m12)
122mu1a<-getContrasts(m12,pickCoef(m12,"[.]Clas.c1"),
123ref = "mean", scaleRef = "mean",
124scaleWeights = "unit")
125nu1a<-getContrasts(m12,pickCoef(m12,"[.]Parti.c1"),
126ref = "mean", scaleRef = "mean",
127scaleWeights = "unit")
128mu2a<-getContrasts(m12,pickCoef(m12,"[.]Clas.c2"),
129ref = "mean", scaleRef = "mean",
130scaleWeights = "unit")
131nu2a<-getContrasts(m12,pickCoef(m12,"[.]Parti.c2"),
132ref = "mean", scaleRef = "mean",
133scaleWeights = "unit")
134cona<-c(mu1a$qvframe[,1][c(1,8)],nu1a$qvframe[,1][c(1,4)],mu2a$qvframe[,1][c(1,8)],nu2a$qvframe[,1][c(1,4)])
135## IL FAUT PARFOIS RELANCER PLUSIEURS FOIS LA LIGNE DE COMMANDE POUR QUE LE MODELE m12a CONVERGE ##
136m12a <- gnm(Freq~Clas+Parti+L+Clas*L+Parti*L+Mult(1+nYear,Clas.c1,Parti.c1)+Mult(1,Clas.c2,Parti.c2), constrain=c(36,43,44,47,49,56,57,60), constrainTo=cona, family=poisson,tolerance = 1e-12,iterStart=10,iterMax=50000)
137mu1a<-getContrasts(m12a,pickCoef(m12a,"[.]Clas.c1"),
138ref = "mean", scaleRef = "mean",
139scaleWeights = "unit")
140nu1a<-getContrasts(m12a,pickCoef(m12a,"[.]Parti.c1"),
141ref = "mean", scaleRef = "mean",
142scaleWeights = "unit")
143mu2a<-getContrasts(m12a,pickCoef(m12a,"[.]Clas.c2"),
144ref = "mean", scaleRef = "mean",
145scaleWeights = "unit")
146nu2a<-getContrasts(m12a,pickCoef(m12a,"[.]Parti.c2"),
147ref = "mean", scaleRef = "mean",
148scaleWeights = "unit")
149summary(m12a);mu1a;nu1a;mu2a;nu2a;model.summary(m12a)
150# Modèle 13 - RC(2)-L homogène dans µi1, νj1, µi2, νj2 et ϕ2, avec tendance linéaire dans ϕ1 et contraintes d’égalité -------------------------------------------------------------
151Clas.c1<-Clas
152levels(Clas.c1)<-factor(c(1,2,3,2,4,5,4,4))
153Clas.c1
154Clas.c2<-Clas
155levels(Clas.c2)<-factor(c(1,2,3,4,3,3,5,5))
156Clas.c2
157Parti.c1<-Parti
158levels(Parti.c1)<-factor(c(1,2,2,3))
159Parti.c1
160Parti.c2<-Parti
161levels(Parti.c2)<-factor(c(1,2,1,3))
162Parti.c2
163m13 <- gnm(Freq~Clas+Parti+L+Clas*L+Parti*L+Mult(1+nYear,Clas.c1,Parti.c1)
164+Mult(1,Clas.c2,Parti.c2)
165,family=poisson,tolerance = 1e-12,iterStart=10,iterMax=50000)
166summary(m13)
167model.summary(m13)
168mu1<-getContrasts(m13,pickCoef(m13,"[.]Clas.c1"),
169ref = c(1,2,1,3,1)/8, scaleRef = c(1,2,1,3,1)/8,
170scaleWeights = c(1,2,1,3,1))
171nu1<-getContrasts(m13,pickCoef(m13,"[.]Parti.c1"),
172ref = c(1,2,1)/4, scaleRef = c(1,2,1)/4,
173scaleWeights=c(1,2,1))
174mu2<-getContrasts(m13,pickCoef(m13,"[.]Clas.c2"),
175ref = c(1,1,3,1,2)/8, scaleRef = c(1,1,3,1,2)/8,
176scaleWeights = c(1,1,3,1,2))
177nu2<-getContrasts(m13,pickCoef(m13,"[.]Parti.c2"),
178ref = c(2,1,1)/4, scaleRef = c(2,1,1)/4,
179scaleWeights = c(2,1,1))
180con<-c(mu1$qvframe[,1][c(1,5)],nu1$qvframe[,1][c(1,3)],
181mu2$qvframe[,1][c(1,5)],nu2$qvframe[,1][c(1,3)])
182m13a <- gnm(Freq~Clas+Parti+L+Clas*L+Parti*L+Mult(1+nYear,Clas.c1,Parti.c1)
183+Mult(1,Clas.c2,Parti.c2),
184constrain=c(36,40,41,43,45,49,50,52),constrainTo=con,family=poisson,tolerance = 1e-12,iterStart=10,iterMax=100000)
185summary(m13a);mu1;nu1;mu2;nu2;model.summary(m13a)
Documents annexes
-
Texte intégral disponible sur le portail Cairn (Lien externe)
-
Annexe 1 : Les déterminants socioprofessionnels du vote en Suisse (Annexes électroniques)
-
Annexe 2 : Scripts LEM des modèles d’association utilisés dans l’article (Annexes électroniques)
-
Annexe 3 : Scripts R des modèles d’association utilisés dans l’article (Annexes électroniques)
Pour citer cet article
Référence électronique
« Annexe 3 : Scripts R des modèles d’association utilisés dans l’article », Sociologie [En ligne], N° 3, vol. 10 | 2019, mis en ligne le 26 septembre 2019, consulté le 21 mai 2025. URL : http://journals.openedition.org/sociologie/6301
Haut de pageDroits d’auteur
Le texte et les autres éléments (illustrations, fichiers annexes importés), sont « Tous droits réservés », sauf mention contraire.
Haut de page