#sink('xie.output')
library (gnm)

## Xie 1992 paper, Table 3
table3<-structure(c(311, 161, 128, 88, 36, 43, 356, 150, 12, 130, 128,
109, 83, 45, 23, 375, 180, 14, 79, 66, 89, 43, 38, 25, 325, 187,
18, 24, 22, 26, 72, 27, 16, 108, 48, 5, 22, 11, 25, 41, 47, 14,
140, 74, 18, 7, 6, 3, 5, 3, 99, 5, 9, 10, 70, 112, 197, 112,
110, 86, 1506, 802, 96, 44, 47, 113, 64, 80, 81, 839, 685, 114,
1, 1, 4, 4, 4, 40, 22, 15, 56, 105, 59, 40, 38, 40, 27, 36, 22,
4, 72, 113, 86, 37, 68, 74, 138, 88, 18, 19, 37, 64, 17, 55,
77, 93, 79, 26, 9, 9, 10, 38, 38, 27, 22, 18, 9, 8, 14, 20, 23,
95, 52, 38, 24, 14, 3, 0, 4, 2, 10, 461, 5, 8, 19, 26, 54, 103,
36, 92, 156, 339, 235, 68, 11, 34, 61, 22, 74, 286, 189, 209,
107, 1, 2, 4, 1, 7, 73, 9, 11, 47, 52, 30, 10, 26, 8, 24, 33,
32, 5, 15, 27, 19, 24, 13, 47, 89, 49, 10, 13, 14, 10, 5, 6,
44, 40, 28, 3, 3, 3, 2, 20, 3, 17, 13, 14, 0, 2, 4, 4, 8, 9,
22, 18, 17, 6, 0, 0, 0, 1, 4, 92, 5, 5, 3, 11, 27, 16, 33, 31,
132, 188, 159, 33, 7, 12, 11, 22, 20, 144, 104, 109, 42, 0, 2,
1, 0, 1, 21, 5, 4, 8), dim = c(9, 9, 3), dimnames = structure(list(
    origin = c("I", "II", "III", "IVa", "IVb", "IVc", "V/VI",
    "VIIa", "VIIb"), destination = c("I", "II", "III", "IVa",
    "IVb", "IVc", "V/VI", "VIIa", "VIIb"), country = c("England",
    "France", "Sweden")), .Names = c("origin", "destination", "country"
)), class = "table")

table3 <- as.data.frame(table3)
lvl <- levels(table3$origin)
levels(table3$origin) <- levels(table3$destination) <- c(rep(paste(lvl[1:2], collapse = " + "), 2), lvl[3], rep(paste(lvl[4:5], collapse = " + "), 2), lvl[6:9])
table3 <- xtabs(Freq ~ origin + destination + country, data = table3)
levelMatrix <- matrix(c(2, 3, 4, 6, 5, 6, 6,
	 3, 3, 4, 6, 4, 5, 6,
	 4, 4, 2, 5, 5, 5, 5,
	 6, 6, 5, 1, 6, 5, 2,
	 4, 4, 5, 6, 3, 4, 5,
	 5, 4, 5, 5, 3, 3, 5,
	 6, 6, 5, 3, 5, 4, 1), 7, 7, byrow = TRUE)

### Fit the levels models given in Table 3 of Xie (1992)
## Null association between origin and destination
nullModel <- gnm(Freq ~ country:origin + country:destination, family = poisson, data = table3, verbose = FALSE)
nullModel

## H_0 Interaction specified by levelMatrix, common to all countries
commonTopo <- update(nullModel, ~ . + Topo(origin, destination, spec = levelMatrix), verbose = FALSE)
commonTopo

## H_x Interaction specified by levelMatrix, different multiplier for each country
multTopo <- update(nullModel, ~ . + Mult(Exp(country), Topo(origin, destination, spec = levelMatrix)), verbose = FALSE)
multTopo

## H_l Interaction specified by levelMatrix, different effects for each country
separateTopo <- update(nullModel, ~ . + country:Topo(origin, destination, spec = levelMatrix), verbose = FALSE)
separateTopo

# FI_0
FI0 <- gnm(Freq ~ country:origin + country:destination + origin:destination, family = poisson, data = table3)
FI0

# FI_x Cross-nationally log-multiplicative full two-way R and C interaction
FIx <- gnm(Freq ~ country:origin + country:destination + Mult(Exp(country), origin:destination), ofInterest="[.]country", family = poisson, data = table3)
FIx
# contrast with sum of squared = 1
# set ref="mean" is the same as recaling the results by normal score
getContrasts(FIx, ofInterest(FIx), ref="first", scaleRef="first", scaleWeights="unit")

anova(nullModel, commonTopo, multTopo, separateTopo)



# Xie's unidiff.do in r
library(foreign)
mob<-as.data.frame(read.dta('e:\\Huashan\\UM\\542\\Data\\xie\\mobility.dta'))

# R(u) Conditional independence model, df=48
Ru <- gnm(obs~ country:father + country:son, family = poisson, data = mob)
Ru
# df =33
Ru2 <- gnm(obs~ country:father + country:son + Diag(father, son):country, family = poisson, data = mob)
Ru2

# FI(x)
unidiff <- gnm(obs~ country*father + country*son + Diag(father, son):country + Mult(country, father:son), ofInterest="[.]country", family = poisson, data = mob, trace=TRUE)
unidiff
getContrasts(unidiff, pickCoef(unidiff, "[.]country"), ref="first", scaleRef="first", scaleWeights="unit")


#sink()