source("http://www.poleto.com/Catdata.r")
ex1cca.TF<-rbind(c(167,17,19,10,1,3,52,10,11),
                 c(120,22,19, 8,5,1,39,12,12))
ex1cca.catdata<-readCatdata(TF=ex1cca.TF)
ex1cca.catdata
summary(ex1cca.catdata)
ex2cca.TF<-c(7,11,2,3,9,5,0,10,4)
ex2cca.catdata<-readCatdata(TF=ex2cca.TF)
ex2cca.catdata
summary(ex2cca.catdata)

ex1.TF<-rbind(c(167,17,19,10,1,3,52,10,11, 176,24,121, 28,10,12),
              c(120,22,19, 8,5,1,39,12,12, 103, 3, 80, 31, 8,14))
ex1.Zp<-kronecker(t(rep(1,2)),cbind(kronecker(diag(3),rep(1,3)),
  kronecker(rep(1,3),diag(3))))
ex1.Rp<-rbind(c(3,3),c(3,3))
ex1.Zp
ex1.catdata<-readCatdata(TF=ex1.TF,Zp=ex1.Zp,Rp=ex1.Rp)
ex1.catdata
summary(ex1.catdata)
ex2.TF<-c(7,11,2,3,9,5,0,10,4, 8,7,3,0, 0,7,14,7)
ex2.Zp<-cbind(rbind( cbind(kronecker(rep(1,2),diag(3)),rep(0,6)),
                     cbind(matrix(0,3,3),rep(1,3)) ),
              rbind( cbind(rep(1,3),matrix(0,3,3)),
                     cbind(rep(0,6),kronecker(rep(1,2),diag(3))) ) )
ex2.Rp<-c(4,4)
ex2.Zp
ex2.catdata<-readCatdata(TF=ex2.TF,Zp=ex2.Zp,Rp=ex2.Rp)
ex2.catdata
summary(ex2.catdata)

ex1m.TF<-rbind(c(167,17,19,10,1,3,52,10,11, 28,10,12, 0,0, 0),
               c(120,22,19, 8,5,1,39,12,12,103, 3,80,31,8,14))
ex1m.Zp<-cbind(kronecker(rep(1,3),diag(3)), kronecker(diag(3),rep(1,3)),
  kronecker(rep(1,3),diag(3)))
ex1m.Rp<-rbind(c(3,0),c(3,3))
ex1m.Zp
ex1m.catdata<-readCatdata(TF=ex1m.TF,Zp=ex1m.Zp,Rp=ex1m.Rp)
summary(ex1m.catdata)

ex1.satmarml<-satMarML(ex1.catdata,method="NR/FS-MAR")
ex1.satmarml2<-satMarML(ex1.catdata,method="FS-MCAR")
ex1.satmcarml<-satMarML(ex1.catdata,missing="MCAR")
ex1.satmarml
summary(ex1.satmarml) #compare the number of iterations
summary(ex1.satmarml2) #compare the std.errors & augmented freq.
summary(ex1.satmcarml)
ex2.satmarml<-satMarML(ex2.catdata,method="NR/FS-MAR")
ex2.satmarml<-satMarML(ex2.catdata)
ex2.satmcarml<-satMarML(ex2.catdata,method="FS-MCAR",missing="MCAR")
ex2.TF2<-c(7,11,2,3,9,5,1e-5,10,4, 8,7,3,0, 0,7,14,7) #subst.by small value
ex2.catdata2<-readCatdata(TF=ex2.TF2,Zp=ex2.Zp,Rp=ex2.Rp)
ex2.satmarml2<-satMarML(ex2.catdata2,method="NR/FS-MAR")
summary(ex2.satmcarml)
summary(ex2.satmarml) #compare all the results
summary(ex2.satmarml2)

ex1.satmcarwls<-satMcarWLS(ex1.catdata)
ex1.satmcarwls
summary(ex1.satmcarwls) #compare with the ML results under MCAR
ex2.satmcarwls<-satMcarWLS(ex2.catdata)
ex2.satmcarwls2<-satMcarWLS(ex2.catdata2)
summary(ex2.satmcarwls) #negative estimates

ex2.A<-rbind(c(1,1,1,0,0,0,0,0,0),
             c(0,0,0,1,1,1,0,0,0),
             c(1,0,0,1,0,0,1,0,0),
             c(0,1,0,0,1,0,0,1,0) )
ex2.X<-rep(1,2)%x%diag(2)
ex2.linmlmar<-linML(ex2.satmarml2,A=ex2.A,X=ex2.X)
ex2.linmlmar
summary(ex2.linmlmar)
ex2.U<-t(c(1,-1))%x%diag(2)
ex2.linmlmar2<-linML(ex2.satmarml2,A=ex2.A,U=ex2.U)
ex2.linmlmar2
summary(ex2.linmlmar2)

ex1.E<-rbind(c(1,-1,0),c(0,1,-1))
ex1.A<-kronecker(diag(2),kronecker(ex1.E,ex1.E))
ex1.XL<-rep(1,8)
ex1.loglinmlcca<-loglinML(ex1cca.catdata,A=ex1.A,XL=ex1.XL)
ex1.loglinmlmar<-loglinML(ex1.satmarml,A=ex1.A,XL=ex1.XL)
ex1.loglinmlmcar<-loglinML(ex1.satmcarml,A=ex1.A,XL=ex1.XL)
summary(ex1.loglinmlcca)
summary(ex1.loglinmlmar)
summary(ex1.loglinmlmcar)
ex1.XL2<-rep(1,2)%x%diag(4)
ex1.loglinmlmar2<-loglinML(ex1.satmarml,A=ex1.A,XL=ex1.XL2)
ex1.UL<-t(c(1,-1))%x%diag(4)
ex1.loglinmlmar3<-loglinML(ex1.satmarml,A=ex1.A,UL=ex1.UL)
summary(ex1.loglinmlmar2)
summary(ex1.loglinmlmar3)

ex2.linwlsmar<-funlinWLS(model="lin",obj=ex2.satmarml2,A1=ex2.A,X=ex2.X)
ex2.linwlsmar
summary(ex2.linwlsmar)
ex2.linwlsmar2<-funlinWLS(model="lin",obj=ex2.satmarml2,A1=ex2.A,U=ex2.U)
summary(ex2.linwlsmar2)
ex1.loglinwlscca<-funlinWLS(model=c("lin","log"),obj=ex1cca.catdata,A1=ex1.A,XL=ex1.XL)
ex1.loglinwlsmar<-funlinWLS(model=c("lin","log"),obj=ex1.satmarml,A1=ex1.A,XL=ex1.XL)
ex1.loglinwlsmcar<-funlinWLS(model=c("lin","log"),obj=ex1.satmcarml,A1=ex1.A,XL=ex1.XL)
ex1.loglinwlsmar2<-funlinWLS(model=c("lin","log"),obj=ex1.satmarml,A1=ex1.A,XL=ex1.XL2)
ex1.loglinwlsmar3<-funlinWLS(model=c("lin","log"),obj=ex1.satmarml,A1=ex1.A,UL=ex1.UL)
summary(ex1.loglinwlscca)
summary(ex1.loglinwlsmar)
summary(ex1.loglinwlsmcar)
summary(ex1.loglinwlsmar2)
summary(ex1.loglinwlsmar3)
W1<-c(1,0.75,0,0.75,1,0.75,0,0.75,1) #squared weights, Fleiss e Cohen (1973)
W2<-c(1,0.5,0,0.5,1,0.5,0,0.5,1) #absolute weights, Agresti (2002)
ex2.kw1A1<-rbind(
 t(W1),
 rep(1,9),
 kronecker(diag(3),t(rep(1,3))),
 kronecker(t(rep(1,3)),diag(3)) )
ex2.kw2A1<-rbind(
 t(W2),
 rep(1,9),
 kronecker(diag(3),t(rep(1,3))),
 kronecker(t(rep(1,3)),diag(3)) )
ex2.kwA2<-rbind(
 cbind(diag(2),matrix(0,2,6)),
 cbind(matrix(0,9,2),  cbind(kronecker(diag(3),rep(1,3)) , kronecker(rep(1,3),diag(3)))  ) )
ex2.kw1A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W1)) )
ex2.kw2A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W2)) )
ex2.kA4<-t(c(1,-1))
ex2.kappaw1<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"),
  obj=ex2.satmarml,A1=ex2.kw1A1,A2=ex2.kwA2,A3=ex2.kw1A3,A4=ex2.kA4,PI1=-1,X=1)
ex2.kappaw2<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"),
  obj=ex2.satmarml,A1=ex2.kw2A1,A2=ex2.kwA2,A3=ex2.kw2A3,A4=ex2.kA4,PI1=-1,X=1)
ex2.kappaw1
ex2.kappaw2

waldTest(ex1.loglinmlmar2,C=diag(4))
waldTest(ex1.loglinwlsmar2,C=diag(4))