Skip to navigation – Site map

HomeNumé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

Full text

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)

Top of page

References

Electronic reference

« Annexe 3 : Scripts R des modèles d’association utilisés dans l’article », Sociologie [Online], N° 3, vol. 10 |  2019, Online since 26 September 2019, connection on 25 June 2025. URL : http://journals.openedition.org/sociologie/6301

Top of page

Copyright

The text and other elements (illustrations, imported files) are “All rights reserved”, unless otherwise stated.

Top of page

Search OpenEdition Search

You will be redirected to OpenEdition Search