Consider the following dataset, from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy, with 563 observations, and nine variables : eight covariates, and the variable of interest, the number of extramarital affairs, over a year,
base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE) str(base) 'data.frame': 563 obs. of 9 variables: $ SEX : int 1 0 0 1 1 0 0 1 0 1 ... $ AGE : num 37 27 32 57 22 32 22 57 32 22 ... $ YEARMARRIAGE: num 10 4 15 15 0.75 1.5 0.75 15 15 1.5 ... $ CHILDREN : int 0 0 1 1 0 0 0 1 1 0 ... $ RELIGIOUS : int 3 4 1 5 2 2 2 2 4 4 ... $ EDUCATION : int 18 14 12 18 17 17 12 14 16 14 ... $ OCCUPATION : int 7 6 1 6 6 5 1 4 1 4 ... $ SATISFACTION: int 4 4 4 5 3 5 3 4 2 5 ... $ Y : int 0 0 0 0 0 0 0 0 0 0 ... |
Let us focus on two categorical covariates, related to the importance of religion, and the occupation
df=data.frame(y=base$Y, religion=as.factor(base$RELIGIOUS), occupation=as.factor(base$OCCUPATION), expo = 1) (E=xtabs(expo~religion+occupation,data=df)) occupation religion 1 2 3 4 5 6 7 1 4 1 8 4 16 9 0 2 23 3 11 17 56 36 6 3 29 1 10 12 39 25 2 4 38 7 12 21 59 44 2 5 13 1 3 10 19 19 3 (N=xtabs(y~religion+occupation,data=df)) occupation religion 1 2 3 4 5 6 7 1 4 1 13 3 13 7 0 2 1 1 13 10 25 43 10 3 15 0 12 11 34 35 1 4 24 1 3 15 11 9 10 5 6 0 0 6 11 7 0 |
The two tables above are the exposure (number of observations) and the number of extramarital affairs, here as contingency tables. Without any covariate, one can assume that N\sim\mathcal{P}(\lambda\cdot E), where \lambda would be
sum(N)/sum(E) [1] 0.6305506 |
The idea with the margin method is to assume that N_{i,j}=E_{i,j}\cdot\lambda_{i,j} where \lambda_{i,j}=A_i\cdot B_j. Bailey (1963) added two series of constraints : per row, \sum_j N_{i,j}=\sum_j E_{i,j}\cdot A_i\cdot B_j for any i and similarly, for any j \sum_i N_{i,j}=\sum_i E_{i,j}\cdot A_i\cdot B_jFrom the first series of constraints, write A_i=\frac{\sum_j N_{i,j}}{\sum_j E_{i,j}\cdot B_j} and use the second series to write B_j=\frac{\sum_i N_{i,j}}{\sum_i E_{i,j}\cdot A_i}Because we need A_i‘s to compute B_j‘s, and conversely, it is natural to consider some iterative procedure to solve it. Observe that we do not have unicity…
Consider here some starting values for A_i‘s and B_j‘s
A=rep(1,length(levels(df$religion))) B=rep(1,length(levels(df$occupation)))*sum(N)/sum(E) A [1] 1 1 1 1 1 B [1] 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 |
The predicted number of extramarital affairs would be \hat N_{i,j}=E_{i,j}\cdot\hat A_i\cdot \hat B_j
E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.5222025 0.6305506 5.0444050 2.5222025 10.0888099 5.6749556 0.0000000 2 14.5026643 1.8916519 6.9360568 10.7193606 35.3108348 22.6998224 3.7833037 3 18.2859680 0.6305506 6.3055062 7.5666075 24.5914742 15.7637655 1.2611012 4 23.9609236 4.4138544 7.5666075 13.2415631 37.2024867 27.7442274 1.2611012 5 8.1971581 0.6305506 1.8916519 6.3055062 11.9804618 11.9804618 1.8916519 sum(B*E[1,]) [1] 26.48313 sum(B*E[2,]) [1] 95.84369 apply(t(B*t(E)),1,sum) 1 2 3 4 5 26.48313 95.84369 74.40497 115.39076 42.87744 sum(A*E[,1]) [1] 107 sum(A*E[,2]) [1] 13 apply(A*E,2,sum) 1 2 3 4 5 6 7 107 13 44 64 189 133 13 |
From expressions above, observe that one can very easily write expressions of A_i‘s and B_j‘s as functions of B_j‘s and A_i‘s respectively
A=apply(N,1,sum)/apply(t(B*t(E)),1,sum) B=apply(N,2,sum)/apply(A*E,2,sum) |
Let it iterate one thousand times
for(i in 1:1000){ A=apply(N,1,sum)/apply(t(B*t(E)),1,sum) B=apply(N,2,sum)/apply(A*E,2,sum) } |
We obtain here
A
1 2 3 4 5
1.5404346 1.0447195 1.4825650 0.6553159 0.6634763
B
1 2 3 4 5 6 7
0.4685515 0.2629769 0.8454435 0.7245310 0.4889697 0.7770553 1.6753750
E * A%*%t(B)
occupation
religion 1 2 3 4 5 6 7
1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000
2 11.2586111 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
3 20.1450811 0.3898804 12.5342484 12.8899708 28.2722423 28.8008726 4.9677044
4 11.6678702 1.2063307 6.6483904 9.9707299 18.9053460 22.4055332 2.1957997
5 4.0413463 0.1744790 1.6827951 4.8070914 6.1639760 9.7955975 3.3347148 |
That is our prediction, per category, of the number of affairs. Observe that here, sums per row are equal to observed numbers,
apply(N,1,sum) 1 2 3 4 5 41 103 108 73 30 apply(E * A%*%t(B),1,sum) 1 2 3 4 5 41 103 108 73 30 |
as well as sums per colums
apply(N,2,sum) 1 2 3 4 5 6 7 50 3 41 45 94 101 21 apply(E * A%*%t(B),2,sum) 1 2 3 4 5 6 7 50 3 41 45 94 101 21 |
Now, why should I mention that here, in the section on the Poisson regression in our course ? Because actually, this is exactly what we get if we run a Poisson regression on those two covariates
reg=glm(y~religion+occupation,data=df,family=poisson) summary(reg) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.32604 0.21325 -1.529 0.126285 religion2 -0.38832 0.18791 -2.066 0.038783 * religion3 -0.03829 0.18585 -0.206 0.836771 religion4 -0.85470 0.19757 -4.326 1.52e-05 *** religion5 -0.84233 0.24416 -3.450 0.000561 *** occupation2 -0.57758 0.59549 -0.970 0.332083 occupation3 0.59022 0.21349 2.765 0.005699 ** occupation4 0.43588 0.20603 2.116 0.034381 * occupation5 0.04265 0.17590 0.242 0.808399 occupation6 0.50587 0.17360 2.914 0.003569 ** occupation7 1.27415 0.26298 4.845 1.27e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
First of all, observe that the total sum of predictions equals the total sum of observations
yp = predict(reg,type="response") sum(yp) [1] 355 sum(df$y) [1] 355 |
But actually, the predicted number of affairs, for our 35 classes, is exactly what we got using Bailey’s technique
xtabs(yp~df$religion+df$occupation) df$occupation df$religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586112 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450813 0.3898804 12.5342484 12.8899708 28.2722424 28.8008726 4.9677044 4 11.6678703 1.2063307 6.6483904 9.9707300 18.9053460 22.4055332 2.1957997 5 4.0413464 0.1744790 1.6827951 4.8070914 6.1639761 9.7955975 3.3347148 E * A%*%t(B) occupation religion 1 2 3 4 5 6 7 1 2.8870914 0.4050987 10.4188024 4.4643702 12.0516123 10.7730250 0.0000000 2 11.2586111 0.8242113 9.7157637 12.8678376 28.6068235 29.2249717 10.5017811 3 20.1450811 0.3898804 12.5342484 12.8899708 28.2722423 28.8008726 4.9677044 4 11.6678702 1.2063307 6.6483904 9.9707299 18.9053460 22.4055332 2.1957997 5 4.0413463 0.1744790 1.6827951 4.8070914 6.1639760 9.7955975 3.3347148 |
To be more specific, up to a multiplicate constant, series of coefficients are equal here, e.g. for A_i‘s
a=exp(coefficients(reg)[1]+c(0,coefficients(reg)[2:5])) a/a[1] religion2 religion3 religion4 religion5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 A/A[1] 1 2 3 4 5 1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 |
but also for B_j‘s
b=exp(coefficients(reg)[1]+c(0,coefficients(reg)[6:11])) b/b[1] occupation2 occupation3 occupation4 occupation5 occupation6 occupation7 1.0000000 0.5612551 1.8043769 1.5463210 1.0435773 1.6584203 3.5756477 B/B[1] 1 2 3 4 5 6 7 1.0000000 0.5612551 1.8043770 1.5463210 1.0435773 1.6584203 3.5756478 |
This will have major implications in non-life insurance models (for claims reserving).