##################################
# Importation
##################################

library(foreign)
he<-read.xport("d:/olivier/heran.xpt")


##################################
#Vérification de la table
##################################

str(he)


##################################
# Création de nouvelles variables 
##################################

he$txchat<-ifelse(he$CHAT=='00',0,1)

he$txchien<-ifelse(he$CHIEN=='00',0,1)

he$chatn<-as.numeric(as.character(he$CHAT))
he$chienn<-as.numeric(as.character(he$CHIEN))
he$cagen<-as.numeric(as.character(he$CAGE))

he$txchat2<-((he$A381E=='1')+(he$A382E=='1')+(he$A383E=='1')+(he$A384E=='1')+(he$A385E=='1')>0)*1
he$txchien2<-((he$A381E=='2')+(he$A382E=='2')+(he$A383E=='2')+(he$A384E=='2')+(he$A385E=='2')>0)*1
he$txoiseaux<-((he$A381E=='3')+(he$A382E=='3')+(he$A383E=='3')+(he$A384E=='3')+(he$A385E=='3')>0)*1
he$txrongeurs<-((he$A381E=='4')+(he$A382E=='4')+(he$A383E=='4')+(he$A384E=='4')+(he$A385E=='4')>0)*1
he$txinconnu<-((he$A381E=='5')+(he$A382E=='5')+(he$A383E=='5')+(he$A385E=='5')+(he$A385E=='5')>0)*1
he$txpoissons<-((he$A381E=='6')+(he$A382E=='6')+(he$A383E=='6')+(he$A385E=='6')+(he$A385E=='6')>0)*1
he$txtortues<-((he$A381E=='7')+(he$A382E=='7')+(he$A383E=='7')+(he$A385E=='7')+(he$A385E=='7')>0)*1
he$txinconnu2<-((he$A381E=='8')+(he$A382E=='8')+(he$A383E=='8')+(he$A385E=='8')+(he$A385E=='8')>0)*1

he$tx_anim<-(he$A38=='1')*1


# Méthode 1 : assigner des valeurs dans certains cas de figures

he$cdip_b[he$CDIP=="1"]<-"1.Aucun Dip"
he$cdip_b[he$CDIP=="2"]<-"2.CE"
he$cdip_b[he$CDIP=="3"]<-"3.CAP"
he$cdip_b[he$CDIP=="4"]<-"4.BEPC"
he$cdip_b[he$CDIP=="5"]<-"5.Bac"
he$cdip_b[he$CDIP=="6"]<-"6.BacPro"
he$cdip_b[he$CDIP=="7"]<-"7.Supérieur"
he$cdip_b<-factor(he$cdip_b)

# Méthode 2 : changer les étiquettes
he$st_b<-he$ST
# on copie la variable et ensuite on attribue les étiquettes
levels(he$st_b)<-c("0.CantonRural","1.CantonRurbain","2.Ville<20 000","3.Ville20-100 000","4.Ville>100 000","5.Banlieue","6.Paris")
# attention de respecter scrupuleusement l'ordre des modalités


# Méthode 3 : ifelse (pratique pour les variables dichotomiques, moins pratiques pour les variables complexes)
# Attention aux valeurs manquantes !
he$ti_b<-ifelse(he$TI=="1","a.Maison",ifelse(he$TI=="2","b.Immeuble","c.nsp"))
he$ti_b<-factor(he$ti_b)

#Retour à la méthode 1. Pratique ici pour un codage plus complèxe
he$tnlr_b[he$TNLR=="1"]<-"a.Maison Indiv"
he$tnlr_b[he$TNLR=="2"]<-"b.2-4 Logements"
he$tnlr_b[he$TNLR=="3"]<-"c.5-9 logements"
he$tnlr_b[he$TNLR=="4"]<-"d.10 à 19 logements"
he$tnlr_b[he$TNLR %in% c("5","6")]<-"e.20 logements ou plus"
he$tnlr_b[he$TNLR=="9"]<-"f.Nb log ND"
he$tnlr_b<-factor(he$tnlr_b)


#Variante de la méthode 2 : 
he$th_b<-factor(he$TH,labels=c("a.Habitat rural","b.Habitat Indiv Agglo","c.Immeubles sur rue","d.Cités, gds ensembles","e.Habitat mixte","f.Habitat ND"))


#Retour à la méthode 1
he$nphr_b[he$NPHR=="01"]<-"1 pièce"
he$nphr_b[he$NPHR=="02"]<-"2 pièces"
he$nphr_b[he$NPHR=="03"]<-"3 pièces"
he$nphr_b[he$NPHR=="04"]<-"4 pièces"
he$nphr_b[he$NPHR=="05"]<-"5 pièces"
he$nphr_b[he$NPHR=="06"]<-"6 pièces"
he$nphr_b[as.numeric(as.character(he$NPHR))>6 & as.numeric(as.character(he$NPHR))<=98]<-"7 pièces ou plus"
he$nphr_b[he$NPHR %in% c("00","99")]<-"Nb pièces ND"
he$nphr_b<-factor(he$nphr_b)


#Méthode 4 : On coupe la variable en morceaux
he$cs6_b<-cut(as.numeric(as.character(he$CCS82)),c(0,19,29,39,49,59,69,99))
levels(he$cs6_b)<-c("1.Agri","2.ArtC","3.Cadr","4.PInt","5.Empl","6.Ouvr","7.Inac")


#Méthode 4 : parfaite pour l'âge.
#Pour utiliser les cut, il est préférable de connaître les minimums et les maximums.
he$cage_b<-cut(as.numeric(as.character(he$CAGE)),c(17,34,54,64,98))
levels(he$cage_b)<-c("18-34 ans","35-54 ans","55-64 ans","65 ans et plus")

#Méthode 3 s'impose ici, comme il n'y a pas de troisième sexe
he$csex_b<-ifelse(he$CSEX=="1","homme","femme")
he$csex_b<-factor(he$csex_b)


#Pourquoi utiliser ici la méthode 1 ? Ce n'est pas la plus économe, certes... mais on est sûr de ne pas faire d'erreurs !
he$ccs82_b[he$CCS82=="00"]<-"CS_SO."
he$ccs82_b[he$CCS82=="11"]<-"AgriP"
he$ccs82_b[he$CCS82=="12"]<-"AgriM"
he$ccs82_b[he$CCS82=="13"]<-"AgriG"
he$ccs82_b[he$CCS82=="21"]<-"Arti"
he$ccs82_b[he$CCS82=="22"]<-"Commer"
he$ccs82_b[he$CCS82=="23"]<-"Indus"
he$ccs82_b[he$CCS82=="31"]<-"PrLib"
he$ccs82_b[he$CCS82=="33"]<-"CadrPub"
he$ccs82_b[he$CCS82=="34"]<-"Profs"
he$ccs82_b[he$CCS82=="35"]<-"Info"
he$ccs82_b[he$CCS82=="37"]<-"CadrPri"
he$ccs82_b[he$CCS82=="38"]<-"Ingénieurs"
he$ccs82_b[he$CCS82=="42"]<-"Instits"
he$ccs82_b[he$CCS82=="43"]<-"PISanté"
he$ccs82_b[he$CCS82=="44"]<-"Clergés"
he$ccs82_b[he$CCS82=="45"]<-"PIPub"
he$ccs82_b[he$CCS82=="46"]<-"PIPri"
he$ccs82_b[he$CCS82=="47"]<-"Techni"
he$ccs82_b[he$CCS82=="48"]<-"Contremaitres"
he$ccs82_b[he$CCS82=="52"]<-"EmPub"
he$ccs82_b[he$CCS82=="53"]<-"ArméePol"
he$ccs82_b[he$CCS82=="54"]<-"EmAdPri"
he$ccs82_b[he$CCS82=="55"]<-"EmCom"
he$ccs82_b[he$CCS82=="56"]<-"EmSePart"
he$ccs82_b[he$CCS82=="62"]<-"OQIndus"
he$ccs82_b[he$CCS82=="63"]<-"OQArtis"
he$ccs82_b[he$CCS82=="64"]<-"Chauffeurs"
he$ccs82_b[he$CCS82=="65"]<-"OQManut"
he$ccs82_b[he$CCS82=="67"]<-"ONQIndus"
he$ccs82_b[he$CCS82=="68"]<-"ONQArt"
he$ccs82_b[he$CCS82=="69"]<-"OAgri"
he$ccs82_b[he$CCS82=="71"]<-"AAgri"
he$ccs82_b[he$CCS82=="72"]<-"AArt"
he$ccs82_b[he$CCS82=="74"]<-"ACad"
he$ccs82_b[he$CCS82=="75"]<-"API"
he$ccs82_b[he$CCS82=="77"]<-"AEmp"
he$ccs82_b[he$CCS82=="78"]<-"AOuv"
he$ccs82_b[he$CCS82=="81"]<-"Chomeurs"
he$ccs82_b[he$CCS82=="83"]<-"Service"
he$ccs82_b[he$CCS82=="84"]<-"Etudes"
he$ccs82_b[he$CCS82=="85"]<-"Inact-60"
he$ccs82_b[he$CCS82=="86"]<-"Inact+60"
he$ccs82_b[he$CCS82=="99"]<-"N.D."
he$ccs82_b<-factor(he$ccs82_b)

he$nbpd_b[he$NBPD=="01"]<-"Personne seule"
he$nbpd_b[02<=as.numeric(as.character(he$NBPD)) & as.numeric(as.character(he$NBPD))<=13]<-"Ménage plusieurs personnes"
he$nbpd_b<-factor(he$nbpd_b)

he$nbpd3_b[he$NBPD6=="1"]<-"1 personne /ménage "
he$nbpd3_b[he$NBPD6=="2"]<-"2 personnes /ménage"
he$nbpd3_b[as.numeric(as.character(he$NBPD6))>2]<-"3 personnes ou + / ménage"
he$nbpd3_b<-factor(he$nbpd3_b)

######################
# Premier paragraphe #
######################

#Nombre moyen d'animaux par ménage
mean(he$tx_anim,na.rm=TRUE)
mean(he$txchat,na.rm=TRUE)
mean(he$txchien,na.rm=TRUE)
mean(he$txoiseaux,na.rm=TRUE)
mean(he$txpoissons,na.rm=TRUE)
mean(he$txrongeurs,na.rm=TRUE)
mean(he$txtortues,na.rm=TRUE)

#Nombre de ménages en France (en utilisant la pondération)
sum(he$PONDER,na.rm=TRUE)

#Nombre de chats
sum(he$PONDER*he$chatn,na.rm=TRUE)

#Nombre de chiens en France
sum(he$PONDER*he$chienn,na.rm=TRUE)

######################
#    Figure 1        #
######################

#Construction de deux tableaux croisés
chat1<-prop.table(table(he$st_b,he$txchat),1)
chien1<-prop.table(table(he$st_b,he$txchien),1)

#Construction d'un tableau effectif
tot<-table(he$st_b)

#On fait les pourcentages en ligne
chat2<-prop.table(chat1,1)
chien2<-prop.table(chien1,1)

#On n'a pas besoin de la première colonne : on ne garde que la deuxième
chat3<-chat2[,2]
chien3<-chien2[,2]

#On assemble les deux colonnes conservées et l'effectif
c0<-cbind(chat3,chien3,tot)

#On transforme le tableau en base de données
c1<-data.frame(c0)

#On change les noms de variables dans la nouvelle base de données
colnames(c1)<-c("chat","chien","tot")

#On crée des nouvelles variables pour la mise en forme graphique
c1$color<-"red"
c1$bg<-"red"
c1$pch<-24

#On regarde la tête de l'ensemble
c1

#La même chose pour la deuxième variable en écrivant le tout de manière plus ramassée
chat<-prop.table(table(he$tnlr_b,he$txchat),1)
chien<-prop.table(table(he$tnlr_b,he$txchien),1)
c2<-data.frame(cbind(chat[,2],chien[,2],table(he$tnlr_b)))
colnames(c2)<-c("chat","chien","tot")
c2$color<-"blue"
c2$bg<-"blue"
c2$pch<-22
c2

#Ici on utilise tapply plutôt que table
chat<-tapply(he$txchat,he$th_b,mean,na.rm=TRUE)	
chien<-tapply(he$txchien,he$th_b,mean,na.rm=TRUE)	
c3<-data.frame(cbind(chat,chien,table(he$th_b)))
colnames(c3)<-c("chat","chien","tot")
c3$color<-"darkgreen"
c3$bg<-"darkgreen"
c3$pch<-21

#Ici aussi on utilise tapply et la moyenen plutôt que table
chat<-tapply(he$txchat,he$nphr_b,mean,na.rm=TRUE)	
chien<-tapply(he$txchien,he$nphr_b,mean,na.rm=TRUE)	
c4<-data.frame(cbind(chat,chien,table(he$nphr_b)))
colnames(c4)<-c("chat","chien","tot")
c4$bg<-"white"
c4$color<-"black"
c4$pch<-22
	
#On assemble le tout
fig1<-rbind(c1,c2,c3,c4)

#Le graphique progressivement enrichi
plot(fig1$chien,fig1$chat,xlab="Taux de possession de chiens",ylab="Taux de possession de chats",pch=fig1$pch,bg=fig1$bg,col=fig1$color,cex=0.001*fig1$tot)
text(fig1$chien,fig1$chat,row.names(fig1),pos=2,col=fig1$color)
lines(fig1$chien[8:12],fig1$chat[8:12],col="blue")
lines(fig1$chien[20:26],fig1$chat[20:26],col="black",lty=2)
abline(0,1,col="gray50",lty=2)
abline(0,0.5,col="gray50",lty=2)
text(0.4,0.4,"chien/chat=1",pos=3,col="gray50")
text(0.54,0.27,"chien/chat=2",col="gray50")


################
#   Figure 2   #
################

##Les taux de posssesion moyens par catégories. 
tapply(he$txchat,list(he$cage_b,he$csex_b,he$nbpd_b),mean)
tapply(he$txchien,list(he$cage_b,he$csex_b,he$nbpd_b),mean)


# On arrive pas vraiment à faire le même graphique


barplot(tapply(he$txchat[he$CSEX=="1" & he$NBPD =="01"],he$cage_b[he$CSEX=="1" & he$NBPD =="01"],mean),
xlab="Age",ylab="Taux de possession de chats",main="Homme seul")

barplot(tapply(he$txchien[he$CSEX=="1" & he$NBPD =="01"],he$cage_b[he$CSEX=="1" & he$NBPD =="01"],mean),
xlab="Age",ylab="Taux de possession de chiens",main="Homme seul")

barplot(tapply(he$txchat[he$CSEX=="2" & he$NBPD =="01"],he$cage_b[he$CSEX=="2" & he$NBPD =="01"],mean),
xlab="Age",ylab="Taux de possession de chats",main="Femme seule")

barplot(tapply(he$txchien[he$CSEX=="2" & he$NBPD =="01"],he$cage_b[he$CSEX=="2" & he$NBPD =="01"],mean),
xlab="Age",ylab="Taux de possession de chiens",main="Femme seule")

barplot(tapply(he$txchat[he$CSEX=="1" & he$NBPD !="01"],he$cage_b[he$CSEX=="1" & he$NBPD !="01"],mean),
xlab="Age",ylab="Taux de possession de chats",main="CM Homme. Ménage >1")

barplot(tapply(he$txchien[he$CSEX=="1" & he$NBPD !="01"],he$cage_b[he$CSEX=="1" & he$NBPD !="01"],mean),
xlab="Age",ylab="Taux de possession de chiens",main="CM Homme. Ménage >1")

barplot(tapply(he$txchat[he$CSEX=="2" & he$NBPD !="01"],he$cage_b[he$CSEX=="2" & he$NBPD !="01"],mean),
xlab="Age",ylab="Taux de possession de chats",main="CM Femme. Ménage >1")

barplot(tapply(he$txchien[he$CSEX=="2" & he$NBPD !="01"],he$cage_b[he$CSEX=="2" & he$NBPD !="01"],mean),
xlab="Age",ylab="Taux de possession de chiens",main="CM Femme. Ménage >1")



#################
#  Figure 3     #
#################

chienmen<-tapply(he$chienn,list(he$NBPD6,he$csex_b),mean)
chatmen<-tapply(he$chatn,list(he$NBPD6,he$csex_b),mean)
chienind<-tapply(he$chienn/as.numeric(as.character(he$NBPD6)),list(he$NBPD6,he$csex_b),mean)
chatind<-tapply(he$chatn/as.numeric(as.character(he$NBPD6)),list(he$NBPD6,he$csex_b),mean)
tablhaut<-cbind(chienmen[,2],chienind[,2],chatmen[,2],chatind[,2])
colnames(tablhaut)<-c("Nb chiens/ménage","Nb chiens/individus","Nb chats/ménage","Nb chats/individus")
rownames(tablhaut)<-paste("H / mén taille",row.names(tablhaut))

chienmen<-tapply(he$chienn,list(he$csex_b),mean)
chatmen<-tapply(he$chatn,list(he$csex_b),mean)
chienind<-tapply(he$chienn/as.numeric(as.character(he$NBPD6)),list(he$csex_b),mean)
chatind<-tapply(he$chatn/as.numeric(as.character(he$NBPD6)),list(he$csex_b),mean)
tablmil<-cbind(chienmen[2],chienind[2],chatmen[2],chatind[2])
rownames(tablmil)<-paste("H / Ensemble")

chienmen<-tapply(he$chienn,list(he$nbpd3_b,he$csex_b),mean)
chatmen<-tapply(he$chatn,list(he$nbpd3_b,he$csex_b),mean)
chienind<-tapply(he$chienn/as.numeric(as.character(he$NBPD6)),list(he$nbpd3_b,he$csex_b),mean)
chatind<-tapply(he$chatn/as.numeric(as.character(he$NBPD6)),list(he$nbpd3_b,he$csex_b),mean)
tablbas1<-cbind(chienmen[,1],chienind[,1],chatmen[,1],chatind[,1])
colnames(tablbas1)<-c("Nb chiens/ménage","Nb chiens/individus","Nb chats/ménage","Nb chats/individus")
rownames(tablbas1)<-paste("F / mén taille",row.names(tablbas1))

chienmen<-tapply(he$chienn,list(he$csex_b),mean)
chatmen<-tapply(he$chatn,list(he$csex_b),mean)
chienind<-tapply(he$chienn/as.numeric(as.character(he$NBPD6)),list(he$csex_b),mean)
chatind<-tapply(he$chatn/as.numeric(as.character(he$NBPD6)),list(he$csex_b),mean)
tablbas2<-cbind(chienmen[1],chienind[1],chatmen[1],chatind[1])
rownames(tablbas2)<-paste("F / Ensemble")

chienmen<-mean(he$chienn)
chatmen<-mean(he$chatn)
chienind<-mean(he$chienn/as.numeric(as.character(he$NBPD6)))
chatind<-mean(he$chatn/as.numeric(as.character(he$NBPD6)))
tablbas3<-cbind(chienmen[1],chienind[1],chatmen[1],chatind[1])
rownames(tablbas3)<-paste("Ensemble")

rbind(tablhaut,tablmil,tablbas1,tablbas2,tablbas3)

plot(c(as.numeric(as.factor(rownames(tablhaut))),as.numeric(as.factor(rownames(tablhaut)))),c(tablhaut[,1],tablhaut[,3])
,xlab="Nombre de personnes dans le ménage",ylab="Nombre moyen d'animaux par ménage")
lines(c(as.numeric(as.factor(rownames(tablhaut)))),c(tablhaut[,1]),col="blue")
lines(c(as.numeric(as.factor(rownames(tablhaut)))),c(tablhaut[,3]),col="red")
summary(lm(he$chienn~as.numeric(as.character(he$NBPD6))))
summary(lm(he$chatn~as.numeric(as.character(he$NBPD6))))
abline(0.081052,0.151760,col="gray50",lty=2)
abline(0.12904,0.08470,col="gray50",lty=2)


plot(c(as.numeric(as.factor(rownames(tablbas1))),as.numeric(as.factor(rownames(tablbas1)))),c(tablbas1[,1],tablbas1[,3])
,xlab="Nombre de personnes dans le ménage",ylab="Nombre moyen d'animaux par ménage")
lines(c(as.numeric(as.factor(rownames(tablbas1)))),c(tablbas1[,1]),col="blue")
lines(c(as.numeric(as.factor(rownames(tablbas1)))),c(tablbas1[,3]),col="red")
abline(0.081052,0.151760,col="gray50",lty=2)
abline(0.12904,0.08470,col="gray50",lty=2)

#################
#  Figure 4     #
#################

t1<-prop.table(table(he$cs6_b,he$tx_anim),1)
t2<-prop.table(table(he$cs6_b,he$txchat),1)
t3<-prop.table(table(he$cs6_b,he$txchien),1)
t4<-prop.table(addmargins(table(he$cdip_b,he$tx_anim),1),1)
t5<-prop.table(addmargins(table(he$cdip_b,he$txchat),1),1)
t6<-prop.table(addmargins(table(he$cdip_b,he$txchien),1),1)

tab<-rbind(cbind(t1[,2],t2[,2],t3[,2]),cbind(t4[,2],t5[,2],t6[,2]))

colnames(tab)<-c("Animal domestique","Chat","Chien")

tab

#################
#  Figure 5     #
#################

chien<-tapply(he$txchien,list(he$ccs82_b,he$csex_b),mean)
chat<-tapply(he$txchat,list(he$ccs82_b,he$csex_b),mean)


plot(c(chien[,1],chien[,2]),c(chat[,1],chat[,2]),xlab="Taux de possession de chiens",ylab="Taux de possession de chats",cex=0.01*table(list(he$ccs82_b,he$csex_b)))
text(chien[,1],chat[,1],paste("F-",row.names(chien)),col="magenta",pos=4,cex=0.5)
text(chien[,2],chat[,2],paste("H-",row.names(chien)),col="blue",pos=2,cex=0.5)
abline(0,1,col="gray70")
abline(0,0.5,col="gray50")
abline(0,0.75,col="gray70",lty=2)
text(0.4,0.4,"chien/chat=1",pos=3,col="gray50")
text(0.7,0.35,"chien/chat=2",col="gray50")
text(0.6,0.45,"chien/chat=1.5",col="gray50")


#################
#  Figure 6     #
#################

chat<-tapply(he$txchat[he$CSEX=="1"],he$ccs82_b[he$CSEX=="1"],mean)
chien<-tapply(he$txchien[he$CSEX=="1"],he$ccs82_b[he$CSEX=="1"],mean)
tab<-data.frame(cbind(chat,chien))
tab$ratio<-tab$chien/tab$chat
tab<-tab[order(tab$ratio),]
tab
tab2<-tab[order(1/tab$ratio),]
tab2$color[tab2$ratio<1]<-"pink"
tab2$color[tab2$ratio<1.5 & tab2$ratio>=1]<-"red"
tab2$color[tab2$ratio<2 & tab2$ratio>=1.5]<-"darkred"
tab2$color[tab2$ratio>=2 ]<-"purple"
tab2$z<-1
tab2$nrow<-cumsum(tab2$z)

barplot(tab2$ratio,horiz=TRUE,col=tab2$color)
text(0,tab2$nrow*1.18,row.names(tab2),pos=4)

#Etiquette des catégories pas géniales


#################
#  Figure 7     #
#################

chien<-tapply(he$txchien,list(he$ccs82_b,he$ti_b),mean)
chat<-tapply(he$txchat,list(he$ccs82_b,he$ti_b),mean)

plot(c(chien[,1],chien[,2]),c(chat[,1],chat[,2]),xlab="Taux de possession de chiens",ylab="Taux de possession de chats",cex=0.01*table(list(he$ccs82_b,he$ti_b)))
text(chien[,1],chat[,1],paste("M-",row.names(chien)),col="magenta",pos=3,cex=0.5)
text(chien[,2],chat[,2],paste("I-",row.names(chien)),col="blue",pos=1,cex=0.5)
abline(0,1,col="gray70")
abline(0,0.5,col="gray50")
abline(0,0.75,col="gray70",lty=2)
text(0.4,0.4,"chien/chat=1",pos=3,col="gray50")
text(0.7,0.35,"chien/chat=2",col="gray50")
text(0.6,0.45,"chien/chat=1.5",col="gray50")