Navigation – Plan du site

AccueilNumérosN° 3, vol. 10EnquêtesLes transformations du vote de class(…)Annexe 3 : Scripts R des modèles d’a(…)

Les transformations du vote de classe en Suisse entre 1975 et 2011

Annexe 3 : Scripts R des modèles d’association utilisés dans l’article

Texte 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)

Haut de page

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 page

Droits 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

Rechercher dans OpenEdition Search

Vous allez être redirigé vers OpenEdition Search