{"word_count": 1052, "character_count": 8748}
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 occupationdf=data.frame(y=base$Y,              religion=as.factor(base$RELIGIOUS),              occupation=as.factor(base$OCCUPATION),              expo = 1)(E=xtabs(expo~religion+occupation,data=df))        occupationreligion  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))        occupationreligion  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  0The 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 [latex]N\sim\mathcal{P}(\lambda\cdot E)[/latex], where [latex]\lambda[/latex] would besum(N)/sum(E)[1] 0.6305506The idea with the margin method is to assume that [latex]N_{i,j}=E_{i,j}\cdot\lambda_{i,j}[/latex] where [latex]\lambda_{i,j}=A_i\cdot B_j[/latex]. Bailey (1963) added two series of constraints : per row, [latex display="true"]\sum_j N_{i,j}=\sum_j E_{i,j}\cdot A_i\cdot B_j[/latex] for any [latex]i[/latex] and similarly, for any [latex]j[/latex] [latex display="true"]\sum_i N_{i,j}=\sum_i E_{i,j}\cdot A_i\cdot B_j[/latex]From the first series of constraints, write [latex display="true"]A_i=\frac{\sum_j N_{i,j}}{\sum_j E_{i,j}\cdot B_j}[/latex] and use the second series to write [latex display="true"]B_j=\frac{\sum_i N_{i,j}}{\sum_i E_{i,j}\cdot A_i}[/latex]Because we need [latex]A_i[/latex]'s to compute [latex]B_j[/latex]'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 [latex]A_i[/latex]'s and [latex]B_j[/latex]'sA=rep(1,length(levels(df$religion)))B=rep(1,length(levels(df$occupation)))*sum(N)/sum(E)A[1] 1 1 1 1 1B[1] 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506The predicted number of extramarital affairs would be [latex]\hat N_{i,j}=E_{i,j}\cdot\hat A_i\cdot \hat B_j[/latex]E * A%*%t(B)        occupationreligion          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.8916519sum(B*E[1,])[1] 26.48313sum(B*E[2,])[1] 95.84369apply(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] 107sum(A*E[,2])[1] 13apply(A*E,2,sum)  1   2   3   4   5   6   7 107  13  44  64 189 133  13From expressions above, observe that one can very easily write expressions of [latex]A_i[/latex]'s and [latex]B_j[/latex]'s as functions of [latex]B_j[/latex]'s and [latex]A_i[/latex]'s respectivelyA=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 timesfor(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 hereA        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)        occupationreligion          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.3347148That 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 columsapply(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  21Now, 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 covariatesreg=glm(y~religion+occupation,data=df,family=poisson)summary(reg)Coefficients:            Estimate Std. Error z value Pr(&#62;|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 ‘ ’ 1First of all, observe that the total sum of predictions equals the total sum of observationsyp = predict(reg,type="response")sum(yp)[1] 355sum(df$y)[1] 355But actually, the predicted number of affairs, for our 35 classes, is exactly what we got using Bailey's techniquextabs(yp~df$religion+df$occupation)           df$occupationdf$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.3347148E * A%*%t(B)        occupationreligion          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.3347148To be more specific, up to a multiplicate constant, series of coefficients are equal here, e.g. for [latex]A_i[/latex]'sa=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.4307072but also for [latex]B_j[/latex]'sb=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.5756478This will have major implications in non-life insurance models (for claims reserving).