#20/04/2007 #Seção 2 x<-2 #armazenando um escalar x x<-c(1,9,-15,8) #constrói um vetor x x<-c(-3,x,6) #amplia o vetor anterior x c(-3,x,6) #se não armazenar num objeto, o resultado já é apresentado 1:6 #seqüência de números 6:1 #ao contrário seq(0,1,0.1) #a seqüência começa em 0, vai até 1, com espaçamentos de 0.1 seq(0.1,1,0.2) #não necessariamente termina em 1 rep(2,5) #repete o 1º argumento o número de vezes indicado no 2º argumento rep(c(1,3),4) c(rep(1,3),rep(0,4)) rep(c(0,2,5),c(1,3,4)) #pode-se utilizar um vetor em cada argumento x[4] #extrai o 4º elemento de x x[2:4] #extrai do 2º ao 4º elemento de x x[c(1,3:4)] #extrai os elementos 1, 3 e 4 de x x[-3] #reproduz x sem o 3º elemento x[-c(1,4)] #reproduz x sem o 1º e o 4º elementos x1<-rbind(1:3,c(1,3,-1)) x1 x2<-cbind(c(1,3),c(0,4)) x2 cbind(x2,x1) matrix(1:6,nrow=3) #este comando utiliza o argumento byrow com F (False) como padrão... matrix(1:6,ncol=3) #que faz com que o vetor seja preenchido por coluna na matriz matrix(1:6,ncol=3,byrow=T) #utilizando a opção T (True), preenche-se por linha x1[2,] #extrai a 2ª linha de x1 x1[,2] #extrai a 2ª coluna de x1 x1[,2:3] #extrai uma submatriz de x1 x1[,-1] #outra maneira x1[,c(1,3)] #outra submatriz x1<-matrix(1:9,nrow=3) x2<-matrix(9:1,nrow=3) t(x1) x1 x2 x1*x2 x1%*%x2 solve(x1*x2) diag(x1) diag(c(1,3,4)) diag(2.5) #números quebrados são truncados diag(2) %x% x1 x1 %x% diag(2) #Seção 3 source("http://www.poleto.com/Catdata.r") #ou grave o arquivo deste sítio e leia-o de alguma outra localização, e.g., source("C:/temp/Catdata.r") #Seção 4 e15.TF<-rbind(c(19, 5, 4, 2), c( 5, 8, 0,17), c(11, 6, 7, 6), c( 2, 5, 1,22)) e15.catdata<-readCatdata(TF=e15.TF) e15.catdata #Apresenta as proporções e erros padrões print(e15.catdata,digits=3) #Pode-se escolher o nº de decimais do arredondamento e31.TF<-c(192,1,5,2,146,5,11,12,71) e31.catdata<-readCatdata(TF=e31.TF) e81.U<-rbind(c(0,-1, 0,1,0, 0,0,0), c(0, 0,-1,0,0, 0,1,0), c(0, 0, 0,0,0,-1,0,1)) e81.X<-rbind(c(1,0,0,0,0), c(0,1,0,0,0), c(0,0,1,0,0), c(0,1,0,0,0), c(0,0,0,1,0), c(0,0,0,0,1), c(0,0,1,0,0), c(0,0,0,0,1)) e81.linml1<-linML(obj=e31.catdata,U=e81.U) #formulação em termos de restrições: simetria e81.linml2<-linML(obj=e31.catdata,X=e81.X) #formulação em equações livres: simetria e81.linml1 #as estatísticas de ajuste das duas formulações coincidem, conforme o esperado... e81.linml2 #mas nesta formulação há informação sobre os parâmetros do modelo summary(e81.linml1) #a saída do summary apresenta mais informações summary(e81.linml2) e91.catdata<-readCatdata(TF=c(3,25,32,68)) e91.X<-rbind(c(0,0), c(0,1), c(1,0), c(1,1)) e91.loglinml1<-loglinML(obj=e91.catdata,U=c(1,-1,-1,1)) #formul.em termos de restrs.: indep. e91.loglinml2<-loglinML(obj=e91.catdata,X=e91.X) #formulação em equações livres: independ. e91.loglinml1 e91.loglinml2 summary(e91.loglinml2) e91.X2<-rbind(c(0,0,0), c(0,1,0), c(1,0,0), c(1,1,1)) e91.loglinml3<-loglinML(obj=e91.catdata,X=e91.X2) #modelo log-linear (LL) ordinário saturado e91.loglinml4<-loglinML(obj=e91.catdata,A=c(1,-1,-1,1),XL=1) #modelo LL generalizado saturado e91.loglinml3 e91.loglinml4 qnorm(0.975) #retorna o quantil 97.5% da distribuição normal padrão e91.loglinml4$beta+c(-1,1)*qnorm(0.975)*sqrt(e91.loglinml4$Vbeta) #IC(95%) para o parâmetro round(e91.loglinml4$beta+c(-1,1)*qnorm(0.975)*sqrt(e91.loglinml4$Vbeta),3) #arred.3 dígitos round(exp(e91.loglinml4$beta),3) #estimativa de MV da razão de chances (RR) round(exp(e91.loglinml4$beta+c(-1,1)*qnorm(0.975)*sqrt(e91.loglinml4$Vbeta)),3) #IC(95%) p/RR e112.linwls1<-funlinWLS(model="lin",obj=e31.catdata,U=e81.U) #form.em termos de restrs.: sim. e112.linwls2<-funlinWLS(model="lin",obj=e31.catdata,X=e81.X) #form.em eqs. livres: simetria e112.linwls1 e112.linwls2 summary(e112.linwls1) e91.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e91.catdata, U=c(1,-1,-1,1)) #formulação em termos de restrições: independência e91.loglinwls2<-funlinWLS(model=c("lin","log"),obj=e91.catdata, X=e91.X) #formulação em equações livres: independência e91.loglinwls1 e91.loglinwls2 summary(e91.loglinwls2) e91.loglinwls3<-funlinWLS(model=c("lin","log"),obj=e91.catdata, X=e91.X2) #modelo log-linear ordinário saturado e91.loglinwls4<-funlinWLS(model=c("lin","log"),obj=e91.catdata,A1=c(1,-1,-1,1), XL=1) #modelo log-linear generalizado saturado e91.loglinwls3 e91.loglinwls4 e1112.TF<-c(11,5,0,14,34,7,2,13,11) e1112.catdata<-readCatdata(TF=e1112.TF) e1112.A1<-rbind( c(rep(c(1,0,0,0),2),1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e1112.A2<-rbind( cbind(diag(2),matrix(0,2,6)), cbind(matrix(0,3,2),kronecker(t(rep(1,2)),diag(3))) ) e1112.A3<-cbind( c(1,0),c(1,1),-c(2,1)%*%t(rep(1,3)) ) e1112.A4<-t(c(1,-1)) e1112.A1 e1112.A2 e1112.A3 e1112.A4 e1112.kappawls<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e1112.catdata,A1=e1112.A1,A2=e1112.A2,A3=e1112.A3,A4=e1112.A4,PI1=-1,X=1) e1112.kappawls round(pnorm((e1112.kappawls$beta-0.35)/sqrt(e1112.kappawls$Vbeta)),3) round(e1112.kappawls$beta+c(-1,1)*qnorm(0.975)*sqrt(e1112.kappawls$Vbeta),3) e117.TF<-rbind(c(28,40,68), c( 5,21,49), c( 1, 4,15)) e117.catdata<-readCatdata(TF=e117.TF) e117.A<-kronecker(diag(3),cbind(diag(2),rep(-1,2))) e117.X1<-rbind(c(1,0,0,0), c(0,1,0,0), c(1,0,2,0), c(0,1,1,0), c(1,0,0,2), c(0,1,0,1)) e117.loglinwls<-funlinWLS(model=c("lin","log"),obj=e117.catdata,A1=e117.A, XL=e117.X1) #efeito de linha e117.loglinwls cbind(0*diag(2),diag(2)) waldTest(obj=e117.loglinwls,C=cbind(0*diag(2),diag(2))) #Seção 5 e134.TF<-c(12,4,5,2, 50,31, 27,12) e134.Zp<-cbind(kronecker(diag(2),rep(1,2)),kronecker(rep(1,2),diag(2))) e134.Rp<-c(2,2) e134.catdata<-readCatdata(TF=e134.TF,Zp=e134.Zp,Rp=e134.Rp) e134.catdata summary(e134.catdata) e132.TF<-c(7,11,2,3,9,5,0,10,4, 8,7,3,0, 0,7,14,7) e132.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))) ) ) e132.Rp<-c(4,4) e132.catdata<-readCatdata(TF=e132.TF,Zp=e132.Zp,Rp=e132.Rp) summary(e132.catdata) e134.satmcarml<-satMarML(e134.catdata,missing="MCAR") e134.satmarml<-satMarML(e134.catdata,method="FS-MCAR") e134.satmarml2<-satMarML(e134.catdata,method="NR/FS-MAR") e134.satmcarml summary(e134.satmcarml)#compare as estimativas das probabilidades e summary(e134.satmarml) #freqüencias ampliadas, erros padrões e summary(e134.satmarml2)#os números de iterações e132.satmarml<-satMarML(e132.catdata,method="NR/FS-MAR") e132.satmarml<-satMarML(e132.catdata) e132.satmcarml<-satMarML(e132.catdata,method="FS-MCAR",missing="MCAR") e132.TF2<-c(7,11,2,3,9,5,1e-5,10,4, 8,7,3,0, 0,7,14,7) #subst.zero por valor peq. e132.catdata2<-readCatdata(TF=e132.TF2,Zp=e132.Zp,Rp=e132.Rp) e132.satmarml2<-satMarML(e132.catdata2,method="NR/FS-MAR") summary(e132.satmcarml)#compare todos os resultados summary(e132.satmarml) summary(e132.satmarml2)#avalie o efeito da substituição pelo valor pequeno e132.satmarml$alphast #Tabela 13.3 - EMV das prob.condicionais de omissão e134.satmcarwls<-satMcarWLS(e134.catdata) e134.satmcarwls summary(e134.satmcarwls) #compare com a abordagem por MV e132.satmcarwls<-satMcarWLS(e132.catdata) e132.satmcarwls2<-satMcarWLS(e132.catdata2) #mesmo com a subst.do zero por valor peq. summary(e132.satmcarwls2) #estimativas negativas e132.U<-rbind(c(0, 1,1,-1,0,0,-1, 0), c(0,-1,0, 1,0,1, 0,-1) ) e132.linml<-linML(e132.satmarml2,U=e132.U) #homog.marg.: formul.restrições e132.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) ) e132.X<-rbind(c(1,0), c(0,1), c(1,0), c(0,1) ) e132.linml2<-linML(e132.satmarml2,A=e132.A,X=e132.X) #homog.marg.: formul.eqs.livres e132.linml summary(e132.linml) summary(e132.linml2) e133.TF<-c(77,87,94,70,67,36,143,78, 14,8,3,9, 25,18,43,16, 14,12) e133.Zp<-cbind(kronecker(diag(4),rep(1,2)), kronecker(diag(2),kronecker(rep(1,2),diag(2))), kronecker(diag(2),rep(1,4)) ) e133.Rp<-c(4,4,2) e133.catdata<-readCatdata(TF=e133.TF,Zp=e133.Zp,Rp=e133.Rp) e133.satmcarml<-satMarML(e133.catdata,missing="MCAR") e133.satmarml<-satMarML(e133.catdata) e133.X<-rbind(c( 1, 1, 1, 1, 1, 1), c( 1, 1,-1, 1,-1,-1), c( 1,-1, 1, -1, 1,-1), c( 1,-1,-1, -1,-1, 1), c(-1, 1, 1, -1,-1, 1), c(-1, 1,-1, -1, 1,-1), c(-1,-1, 1, 1,-1,-1), c(-1,-1,-1, 1, 1, 1)) e133.U<-c(1,-1,-1,1, -1,1,1,-1) e133.loglinml<-loglinML(obj=e133.satmcarml,X=e133.X) e133.loglinml2<-loglinML(obj=e133.satmcarml,U=e133.U) e133.loglinmlmar<-loglinML(obj=e133.satmarml,X=e133.X) e133.loglinml #Tabela 13.6 summary(e133.loglinml2) #Tabela 13.7 summary(e133.loglinml) #Tabela 13.6/13.7 summary(e133.loglinmlmar) e132.linwls<-funlinWLS(model="lin",obj=e132.satmarml2,U=e132.U) e132.linwls2<-funlinWLS(model="lin",obj=e132.satmarml2,A1=e132.A,X=e132.X) e132.linwls summary(e132.linwls) summary(e132.linwls2) e133.satmcarwls<-satMcarWLS(e133.catdata) e133.wlswlsmcar<-funlinWLS(model=c("lin","log"),obj=e133.satmcarwls,X=e133.X) e133.mlwlsmcar<-funlinWLS(model=c("lin","log"),obj=e133.satmcarml,X=e133.X) e133.mlwlsmar<-funlinWLS(model=c("lin","log"),obj=e133.satmarml,X=e133.X) summary(e133.wlswlsmcar) summary(e133.mlwlsmcar) summary(e133.mlwlsmar) e132.kA1<-rbind( c(rep(c(1,0,0,0),2),1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e132.kA2<-rbind( cbind(diag(2),matrix(0,2,6)), cbind(matrix(0,3,2),kronecker(t(rep(1,2)),diag(3))) ) e132.kA3<-cbind( c(1,0),c(1,1),-c(2,1)%*%t(rep(1,3)) ) e132.kA4<-t(c(1,-1)) e132.kappa<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e132.satmarml,A1=e132.kA1,A2=e132.kA2,A3=e132.kA3,A4=e132.kA4,PI1=-1,X=1) W1<-c(1,0.75,0,0.75,1,0.75,0,0.75,1) #pesos quadráticos Fleiss e Cohen (1973) W2<-c(1,0.5,0,0.5,1,0.5,0,0.5,1) #pesos absolutos Agresti (2002) e132.kw1A1<-rbind( t(W1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e132.kw2A1<-rbind( t(W2), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e132.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))))) e132.kw1A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W1)) ) e132.kw2A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W2)) ) e132.kappaw1<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e132.satmarml,A1=e132.kw1A1,A2=e132.kwA2,A3=e132.kw1A3,A4=e132.kA4,PI1=-1,X=1) e132.kappaw2<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e132.satmarml,A1=e132.kw2A1,A2=e132.kwA2,A3=e132.kw2A3,A4=e132.kA4,PI1=-1,X=1) e132.kappa e132.kappaw1 e132.kappaw2 e134.TF2<-c(e134.TF,24) #para mecanismos MNAR, cenários de omissão total trazem inf.na estim. mnarsat.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1,\alpha_2 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7];a2<-p[8] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))) )+ n121*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (exp(a20+a1+a2)/(1+exp(a20+a1+a2))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))) )+ n22*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (1/(1+exp(a20+a1+a2))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))* (exp(a30+a1+a2)/(1+exp(a30+a1+a2))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))) + t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) ) value } require(geoR) #.nlmP adapta nlm p/restringir o espaço paramétrico. Isso é importante, inipars<-c(0.25,0.25,0.25,0,0,0,0,0) #pois mecanismos MNAR resultam facilmente em estims. minpars<-c(0,0,0,-Inf,-Inf,-Inf,-Inf,-Inf) #p/probs. >1 ou <0 quando não se usa o EM ou maxpars<-c(1,1,1,Inf,Inf,Inf,Inf,Inf) #funções ligações próprias para probs.(e.g.,logito) mnarsat<-.nlmP(objfunc=mnarsat.mlv,params=inipars,lower=minpars,upper=maxpars,hessian=T) round(mnarsat$est,3) #estimativas de MV round( -mnarsat$min ,3) #valor máximo da log-verossimilhança mnarsat$code #veja o significado dos códigos de retorno em ?nlm mnarsat$it #número de iterações mnarsat.der<-deriv3(~-( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))) )+ n121*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (exp(a20+a1+a2)/(1+exp(a20+a1+a2))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))) )+ n22*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (1/(1+exp(a20+a1+a2))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))* (exp(a30+a1+a2)/(1+exp(a30+a1+a2))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))) + t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) ),c("t11","t12","t21","a10","a20","a30","a1","a2"), c("t11","t12","t21","a10","a20","a30","a1","a2", "n111","n112","n121","n122","n21","n22","n31","n32","N4") ) #obtém o gradiente e a hessiana analiticamente p<-mnarsat$est mnarsat.infobs<-attr(mnarsat.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],e134.TF2[1],e134.TF2[2], e134.TF2[3],e134.TF2[4],e134.TF2[5],e134.TF2[6],e134.TF2[7],e134.TF2[8],e134.TF2[9]), "hessian")[1,,] mnarsat.infobs2<-mnarsat$hess mnarsat.infobs #matriz informação observada estimada obtida analiticamente round(mnarsat.infobs2,7)#matriz informação observada estimada obtida numericamente mnarsat.cov<-solve(mnarsat.infobs) mnarsat.cov2<-solve(mnarsat.infobs2) round(mnarsat.cov,6) #compare os resultados round(mnarsat.cov2,6)#sugere-se usar a matriz analítica round(sqrt(diag(mnarsat.cov)),4) round(sqrt(diag(mnarsat.cov2)),4) mnarsat.esp<-function(p,N){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1,\alpha_2 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7];a2<-p[8] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))), t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))*(exp(a20+a1+a2)/(1+exp(a20+a1+a2))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))), t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))*(1/(1+exp(a20+a1+a2))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))), t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(exp(a30+a1+a2)/(1+exp(a30+a1+a2))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))), t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))), (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) value } mnarsat.esp(p=p,sum(e134.TF2)) #precisa-se organizar para ficar como na Tabela 13.11 matrix(mnarsat.esp(p=p,sum(e134.TF2))[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) b<-c(rep(0,3),1) B<-rbind(diag(3),rep(-1,3)) c(b+B%*%p[1:3]) #vetor completo das estimativas das probs, incluindo 1-p[1]-p[2]-p[3] B%*%mnarsat.cov[1:3,1:3]%*%t(B) #correspondente vetor de covariâncias sqrt(diag(B%*%mnarsat.cov[1:3,1:3]%*%t(B))) e134.A<-rbind(c(1,1,0,0),c(1,0,1,0)) mnarsat.wls<-funlinWLS(model="lin",theta=c(b+B%*%p[1:3]),Vtheta= B%*%mnarsat.cov[1:3,1:3]%*%t(B),A1=e134.A,X=rep(1,2)) mnarsat.wls2<-funlinWLS(model="lin",theta=c(b+B%*%p[1:3]),Vtheta= B%*%mnarsat.cov[1:3,1:3]%*%t(B),A1=t(c(0,1,-1,0)),X=1) mnarsat.wls summary(mnarsat.wls) mnarsat.wls2 #Seção 6 #Exemplo 8.1 (p.228) / 3.1 (p.47): Problema da intenção de voto e81.TF<-c(192,1,5,2,146,5,11,12,71) e81.catdata<-readCatdata(TF=e81.TF) e81.U<-rbind(c(0,-1, 0,1,0, 0,0,0), c(0, 0,-1,0,0, 0,1,0), c(0, 0, 0,0,0,-1,0,1)) e81.X<-rbind(c(1,0,0,0,0), c(0,1,0,0,0), c(0,0,1,0,0), c(0,1,0,0,0), c(0,0,0,1,0), c(0,0,0,0,1), c(0,0,1,0,0), c(0,0,0,0,1)) e81.linml1<-linML(e81.catdata,U=e81.U) #simetria e81.linml2<-linML(e81.catdata,X=e81.X) #simetria #Exemplo 8.2 (p.233) / 3.2 (p.49) / 1.2 (p.4): Problema do risco de cárie dentária e82.TF<-c(11,5,0,14,34,7,2,13,11) e82.catdata<-readCatdata(TF=e82.TF) e82.U<-rbind(c(0, 1,1,-1,0,0,-1, 0), c(0,-1,0, 1,0,1, 0,-1)) e82.X<-rbind(c(1, 0, 0,0,0,0), c(0, 1, 0,0,0,0), c(0,-1, 1,0,1,0), c(0, 0, 1,0,0,0), c(0, 0, 0,1,0,0), c(0, 1,-1,0,0,1), c(0, 0, 0,0,1,0), c(0, 0, 0,0,0,1)) e82.linml1<-linML(e82.catdata,U=e82.U) #homogeneidade marginal e82.linml2<-linML(e82.catdata,X=e82.X) #homogeneidade marginal e82.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) ) e82.U2<-rbind(c(1,0,-1, 0), c(0,1, 0,-1)) e82.X2<-rbind(c(1,0), c(0,1), c(1,0), c(0,1) ) e82.linml3<-linML(e82.catdata,A=e82.A,U=e82.U2) #homogeneidade marginal e82.linml4<-linML(e82.catdata,A=e82.A,X=e82.X2) #homogeneidade marginal #Exemplo 8.3 (p.236) / 3.3 (p.50) / 1.9 (p.12): Problema do tamanho da ninhada e83.TF<-rbind(c(10,21, 96,23), c( 4, 6, 28, 8), c( 9, 7, 58, 7), c( 8,19, 44, 1), c( 5,17, 56, 1), c( 1, 5, 20, 2), c(22,95,103, 4), c(18,49, 62, 0), c( 4,12, 16, 2)) e83.catdata<-readCatdata(TF=e83.TF) e83.A<-kronecker(diag(9),t(c(0,1,2,3))) e83.U<-rbind(c(-1,1,0,1,-1, 0,0, 0, 0), c(-1,1,0,0, 0, 0,1,-1, 0), c(-1,0,1,1, 0,-1,0, 0, 0), c(-1,0,1,0, 0, 0,1, 0,-1)) e83.X<-rbind(c(1,0,0,0,0), c(1,1,0,0,0), c(1,0,1,0,0), c(1,0,0,1,0), c(1,1,0,1,0), c(1,0,1,1,0), c(1,0,0,0,1), c(1,1,0,0,1), c(1,0,1,0,1)) e83.linml1<-linML(e83.catdata,A=e83.A,U=e83.U,epsilon2=1e-4) e83.linml2<-linML(e83.catdata,A=e83.A,X=e83.X,epsilon2=1e-5) waldTest(e83.linml2,rbind(c(0,0,0,1,0),c(0,0,0,0,1))) waldTest(e83.linml2,rbind(c(0,1,0,0,0),c(0,0,1,0,0))) #Exemplo 9.1 (p.263): Problema da anemia e91.TF<-c(3,25,32,68) e91.catdata<-readCatdata(TF=e91.TF) e91.U<-c(1,-1,-1,1) e91.X<-rbind(c(0,0), c(0,1), c(1,0), c(1,1)) e91.X2<-rbind(c(0,0,0), c(0,1,0), c(1,0,0), c(1,1,1)) e91.loglinml1<-loglinML(e91.catdata,U=e91.U) #independência e91.loglinml2<-loglinML(e91.catdata,X=e91.X) #independência e91.loglinml3<-loglinML(e91.catdata,X=e91.X2) #modelo saturado e91.loglinml4<-loglinML(e91.catdata,A=c(1,-1,-1,1),XL=1) #modelo saturado round(e91.loglinml4$beta+c(-1,1)*qnorm(0.975)*sqrt(e91.loglinml4$Vbeta),3) round(exp(e91.loglinml4$beta),3) #razão de chances round(exp(e91.loglinml4$beta+c(-1,1)*qnorm(0.975)*sqrt(e91.loglinml4$Vbeta)),3) #Exemplo 9.2 (p.267): Problema da acuidade visual e92.TF<-c(1520,266,124,66, 234,1512,432,78, 117,362,1772,205, 36,82,179,492) e92.catdata<-readCatdata(TF=e92.TF) e92.X1<-rbind(c(1,0,0,0,0,0,0,0,0), c(0,1,0,0,0,0,0,0,0), c(0,0,1,0,0,0,0,0,0), c(0,0,0,1,0,0,0,0,0), c(0,1,0,0,0,0,0,0,0), c(0,0,0,0,1,0,0,0,0), c(0,0,0,0,0,1,0,0,0), c(0,0,0,0,0,0,1,0,0), c(0,0,1,0,0,0,0,0,0), c(0,0,0,0,0,1,0,0,0), c(0,0,0,0,0,0,0,1,0), c(0,0,0,0,0,0,0,0,1), c(0,0,0,1,0,0,0,0,0), c(0,0,0,0,0,0,1,0,0), c(0,0,0,0,0,0,0,0,1)) e92.linml1<-linML(e92.catdata,X=e92.X1) #simetria em formulação linear e92.A1<-rbind(c(0,1,0,0, -1,0,0,0, 0, 0,0,0, 0, 0, 0,0), c(0,0,1,0, 0,0,0,0, -1, 0,0,0, 0, 0, 0,0), c(0,0,0,1, 0,0,0,0, 0, 0,0,0, -1, 0, 0,0), c(0,0,0,0, 0,0,1,0, 0,-1,0,0, 0, 0, 0,0), c(0,0,0,0, 0,0,0,1, 0, 0,0,0, 0,-1, 0,0), c(0,0,0,0, 0,0,0,0, 0, 0,0,1, 0, 0,-1,0)) e92.linml2<-linML(e92.catdata,U=e92.A1[,1:15]) #simetria em formulação linear #u_1,u_2,u_3, u_{11},u_{12},u_{13},u_{22},u_{23},u_{33} e92.X2<-rbind(c( 2, 0, 0, 1, 0, 0, 0, 0, 0), c( 1, 1, 0, 0, 1, 0, 0, 0, 0), c( 1, 0, 1, 0, 0, 1, 0, 0, 0), c( 0,-1,-1, -1,-1,-1, 0, 0, 0), c( 1, 1, 0, 0, 1, 0, 0, 0, 0), c( 0, 2, 0, 0, 0, 0, 1, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0), c(-1, 0,-1, 0,-1, 0,-1,-1, 0), c( 1, 0, 1, 0, 0, 1, 0, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0), c( 0, 0, 2, 0, 0, 0, 0, 0, 1), c(-1,-1, 0, 0, 0,-1, 0,-1,-1), c( 0,-1,-1, -1,-1,-1, 0, 0, 0), c(-1, 0,-1, 0,-1, 0,-1,-1, 0), c(-1,-1, 0, 0, 0,-1, 0,-1,-1), c(-2,-2,-2, 1, 2, 2, 1, 2, 1)) #análogo à matriz da pág.71 e92.loglinml1<-loglinML(e92.catdata,X=e92.X2) #simetria em formulação log-linear e92.A2<-rbind(cbind(kronecker(diag(3),t(rep(1,4))),matrix(0,3,4)), kronecker(t(rep(1,4)),cbind(diag(3),rep(0,3)))) e92.linml3<-linML(e92.catdata,A=e92.A2,X=kronecker(rep(1,2),diag(3))) #homogeneidade marginal #u_1^A,u_2^A,u_3^A, u_1^B,u_2^B,u_3^B, u_{11},u_{12},u_{13},u_{22},u_{23},u_{33} e92.X3<-rbind(c( 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0), c( 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0), c( 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0), c( 1, 0, 0, -1,-1,-1, -1,-1,-1, 0, 0, 0), c( 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0), c( 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0), c( 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0), c( 0, 1, 0, -1,-1,-1, 0,-1, 0,-1,-1, 0), c( 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0), c( 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0), c( 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1), c( 0, 0, 1, -1,-1,-1, 0, 0,-1, 0,-1,-1), c(-1,-1,-1, 1, 0, 0, -1,-1,-1, 0, 0, 0), c(-1,-1,-1, 0, 1, 0, 0,-1, 0,-1,-1, 0), c(-1,-1,-1, 0, 0, 1, 0, 0,-1, 0,-1,-1), c(-1,-1,-1, -1,-1,-1, 1, 2, 2, 1, 2, 1)) #e92.X2<-cbind(e92.X3[,1]+e92.X3[,4],e92.X3[,2]+e92.X3[,5],e92.X3[,3]+e92.X3[,6],e92.X3[,7:12]) e92.loglinml2<-loglinML(e92.catdata,X=e92.X3) #quasi-simetria #u_1,u_2,u_3, u_{11},u_{12},u_{13},u_{22},u_{23},u_{33},gama e92.X4<-rbind(c( 2, 0, 0, 1, 0, 0, 0, 0, 0, 0), c( 1, 1, 0, 0, 1, 0, 0, 0, 0, 1), c( 1, 0, 1, 0, 0, 1, 0, 0, 0, 1), c( 0,-1,-1, -1,-1,-1, 0, 0, 0, 1), c( 1, 1, 0, 0, 1, 0, 0, 0, 0, 0), c( 0, 2, 0, 0, 0, 0, 1, 0, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0, 1), c(-1, 0,-1, 0,-1, 0,-1,-1, 0, 1), c( 1, 0, 1, 0, 0, 1, 0, 0, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0, 0), c( 0, 0, 2, 0, 0, 0, 0, 0, 1, 0), c(-1,-1, 0, 0, 0,-1, 0,-1,-1, 1), c( 0,-1,-1, -1,-1,-1, 0, 0, 0, 0), c(-1, 0,-1, 0,-1, 0,-1,-1, 0, 0), c(-1,-1, 0, 0, 0,-1, 0,-1,-1, 0), c(-2,-2,-2, 1, 2, 2, 1, 2, 1, 0)) e92.loglinml3<-loglinML(e92.catdata,X=e92.X4) #simetria condicional e92.loglinml4<-loglinML(e92.catdata,A=e92.A1,XL=rep(1,6)) #simetria condicional round(e92.loglinml1$QvH-e92.loglinml3$QvH,3) round(1-pchisq(e92.loglinml1$QvH-e92.loglinml3$QvH,1),3) round(e92.loglinml1$QwH-e92.loglinml3$QwH,3) round(1-pchisq(e92.loglinml1$QwH-e92.loglinml3$QwH,1),3) round(exp(e92.loglinml4$beta),3) round(kronecker(diag(2),rbind(c(1,0,0),c(1,1,0),c(1,1,1)))%*%e92.A2%*%e92.loglinml4$thetaH,3) #Exemplo 9.3 (p.269) / 1.6 (p.11): Problema dos defeitos de fibras têxteis e93.TF<-c(28,40,68,5,21,49,1,4,15) e93.catdata<-readCatdata(TF=e93.TF) e93.U1<-rbind(c(1,-1, 0,-1, 1, 0, 0, 0,0), c(0, 1,-1, 0,-1, 1, 0, 0,0), c(0, 0, 0, 1,-1, 0,-1, 1,0), c(0, 0, 0, 0, 1,-1, 0,-1,1)) e93.X1<-rbind(c( 1, 0, 1, 0), c( 1, 0, 0, 1), c( 1, 0,-1,-1), c( 0, 1, 1, 0), c( 0, 1, 0, 1), c( 0, 1,-1,-1), c(-1,-1, 1, 0), c(-1,-1, 0, 1), c(-1,-1,-1,-1)) e93.loglinml1<-loglinML(e93.catdata,U=e93.U1) #indep. e93.loglinml2<-loglinML(e93.catdata,X=e93.X1) #indep. e93.X2<-rbind(c( 1, 0, 1, 0,-1, 0), c( 1, 0, 0, 1, 0, 0), c( 1, 0,-1,-1, 1, 0), c( 0, 1, 1, 0, 0,-1), c( 0, 1, 0, 1, 0, 0), c( 0, 1,-1,-1, 0, 1), c(-1,-1, 1, 0, 1, 1), c(-1,-1, 0, 1, 0, 0), c(-1,-1,-1,-1,-1,-1)) e93.loglinml3<-loglinML(e93.catdata,X=e93.X2) #efeito de linha round(-2*e93.loglinml3$beta[5]-e93.loglinml3$beta[6],3) #w_3=-w_1-w_2 round(exp(-2*e93.loglinml3$beta[5]-e93.loglinml3$beta[6]),2) round(e93.loglinml1$QvH-e93.loglinml3$QvH,2) round(1-pchisq(e93.loglinml1$QvH-e93.loglinml3$QvH,2),2) round(e93.loglinml1$QwH-e93.loglinml3$QwH,2) round(1-pchisq(e93.loglinml1$QwH-e93.loglinml3$QwH,2),2) e93.X3<-rbind(c( 1, 0, 1, 0, 1), c( 1, 0, 0, 1, 0), c( 1, 0,-1,-1,-1), c( 0, 1, 1, 0, 0), c( 0, 1, 0, 1, 0), c( 0, 1,-1,-1, 0), c(-1,-1, 1, 0,-1), c(-1,-1, 0, 1, 0), c(-1,-1,-1,-1, 1)) e93.loglinml4<-loglinML(e93.catdata,X=e93.X3) #associação uniforme round(e93.loglinml4$QvH-e93.loglinml3$QvH,2) round(1-pchisq(e93.loglinml4$QvH-e93.loglinml3$QvH,1),2) round(e93.loglinml4$QwH-e93.loglinml3$QwH,2) round(1-pchisq(e93.loglinml4$QwH-e93.loglinml3$QwH,1),2) round(exp(2*2*e93.loglinml4$beta[5]),2) round(exp(e93.loglinml4$beta[5]),2) round(exp(e93.loglinml4$beta[5])*exp(c(-1,1)*qnorm(0.975)*sqrt(e93.loglinml4$Vbeta[5,5])),2) round(e93.loglinml1$QvH-e93.loglinml4$QvH,2) round(1-pchisq(e93.loglinml1$QvH-e93.loglinml4$QvH,1),2) round(e93.loglinml1$QwH-e93.loglinml4$QwH,2) round(1-pchisq(e93.loglinml1$QwH-e93.loglinml4$QwH,1),2) #Exemplo 9.4 (p.274) / 1.11 (p.13): Problema da fobia em alcoólatras e94.TF<-c(10,24,6,12,13,17,4,7) e94.catdata<-readCatdata(TF=e94.TF) e94.X<-rbind(c(0,0,0, 0,0,0, 0), c(0,0,1, 0,0,0, 0), c(0,1,0, 0,0,0, 0), c(0,1,1, 0,0,1, 0), c(1,0,0, 0,0,0, 0), c(1,0,1, 0,1,0, 0), c(1,1,0, 1,0,0, 0), c(1,1,1, 1,1,1, 1)) e94.loglinml1<-loglinML(e94.catdata,X=e94.X) #ABC e94.loglinml2<-loglinML(e94.catdata,X=e94.X[,-7]) #AB,AC,BC e94.loglinml3<-loglinML(e94.catdata,X=e94.X[,-(6:7)]) #AB,AC e94.loglinml4<-loglinML(e94.catdata,X=e94.X[,-(5:7)]) #AB,C e94.loglinml5<-loglinML(e94.catdata,X=e94.X[,-(4:7)]) #A,B,C exp(e94.loglinml1$beta[4:6]) #Assoc.parcial ABC nível 1 exp(e94.loglinml1$beta[4:6]+e94.loglinml1$beta[7]) #Assoc.parcial ABC nível 2 exp(e94.loglinml2$beta[4:6]) #Assoc.parcial AB,AC,BC exp(e94.loglinml3$beta[4:5]) #Assoc.parcial AB,AC exp(e94.loglinml4$beta[4]) #Assoc.parcial AB,C exp(t(c(1,-1,-1,1))%*%log(e94.loglinml1$thetaH[c(1,3,5,7)]+e94.loglinml1$thetaH[c(2,4,6,8)]) ) #Assoc.marg.ABC - AB exp(t(c(1,-1,-1,1))%*%log(e94.loglinml1$thetaH[c(1,2,5,6)]+e94.loglinml1$thetaH[c(3,4,7,8)]) ) #Assoc.marg.ABC - AC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml1$thetaH[1:4]+e94.loglinml1$thetaH[5:8]) ) #Assoc.marg.ABC - BC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml2$thetaH[c(1,3,5,7)]+e94.loglinml2$thetaH[c(2,4,6,8)]) ) #Assoc.marg.AB,AC,BC - AB exp(t(c(1,-1,-1,1))%*%log(e94.loglinml2$thetaH[c(1,2,5,6)]+e94.loglinml2$thetaH[c(3,4,7,8)]) ) #Assoc.marg.AB,AC,BC - AC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml2$thetaH[1:4]+e94.loglinml2$thetaH[5:8]) ) #Assoc.marg.AB,AC,BC - BC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml3$thetaH[c(1,3,5,7)]+e94.loglinml3$thetaH[c(2,4,6,8)]) ) #Assoc.marg.AB,AC - AB exp(t(c(1,-1,-1,1))%*%log(e94.loglinml3$thetaH[c(1,2,5,6)]+e94.loglinml3$thetaH[c(3,4,7,8)]) ) #Assoc.marg.AB,AC - AC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml3$thetaH[1:4]+e94.loglinml3$thetaH[5:8]) ) #Assoc.marg.AB,AC - BC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml4$thetaH[c(1,3,5,7)]+e94.loglinml4$thetaH[c(2,4,6,8)]) ) #Assoc.marg.AB,C - AB exp(t(c(1,-1,-1,1))%*%log(e94.loglinml4$thetaH[c(1,2,5,6)]+e94.loglinml4$thetaH[c(3,4,7,8)]) ) #Assoc.marg.AB,C - AC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml4$thetaH[1:4]+e94.loglinml4$thetaH[5:8]) ) #Assoc.marg.AB,C - BC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml5$thetaH[c(1,3,5,7)]+e94.loglinml5$thetaH[c(2,4,6,8)]) ) #Assoc.marg.A,B,C - AB exp(t(c(1,-1,-1,1))%*%log(e94.loglinml5$thetaH[c(1,2,5,6)]+e94.loglinml5$thetaH[c(3,4,7,8)]) ) #Assoc.marg.A,B,C - AC exp(t(c(1,-1,-1,1))%*%log(e94.loglinml5$thetaH[1:4]+e94.loglinml5$thetaH[5:8]) ) #Assoc.marg.A,B,C - BC e94.X2<-rbind(c( 1, 1, 1, 1, 1, 1, 1), c( 1, 1,-1, 1,-1,-1,-1), c( 1,-1, 1,-1, 1,-1,-1), c( 1,-1,-1,-1,-1, 1, 1), c(-1, 1, 1,-1,-1, 1,-1), c(-1, 1,-1,-1, 1,-1, 1), c(-1,-1, 1, 1,-1,-1, 1), c(-1,-1,-1, 1, 1, 1,-1)) #para obter os resultados da Tabela 9.10 e94.loglinml6<-loglinML(e94.catdata,X=e94.X2) #ABC e94.loglinml7<-loglinML(e94.catdata,X=e94.X2[,-(6:7)]) #AB,AC e94.loglinml8<-loglinML(e94.catdata,X=e94.X2[,-(4:7)]) #A,B,C #Exemplo 9.5 (p.278) / Exercício 8.12 (p.244): Problema da obesidade juvenil e95.TF<-c(300,17,18,19,18,7,17,52) e95.catdata<-readCatdata(TF=e95.TF) #u_1^A,u_1^B,u_1^C,u_{11},u_{111} e95.X1<-rbind(c( 1, 1, 1, 3, 1), c( 1, 1,-1,-1,-1), c( 1,-1, 1,-1,-1), c( 1,-1,-1,-1, 1), c(-1, 1, 1,-1,-1), c(-1, 1,-1,-1, 1), c(-1,-1, 1,-1, 1), c(-1,-1,-1, 3,-1)) e95.loglinml1<-loglinML(e95.catdata,X=e95.X1) #quasi-simetria round(1-rbind(c(rep(1,4),rep(0,4)),rep(c(1,1,0,0),2),rep(c(1,0),4))%*%e95.loglinml1$thetaH,3) #u_1,u_{11},u_{111} e95.X2<-rbind(c( 3, 3, 1), c( 1,-1,-1), c( 1,-1,-1), c(-1,-1, 1), c( 1,-1,-1), c(-1,-1, 1), c(-1,-1, 1), c(-3, 3,-1)) e95.loglinml2<-loglinML(e95.catdata,X=e95.X2) #simetria completa (ver exerc.4.18-c) round(e95.loglinml2$QvH-e95.loglinml1$QvH,2) round(1-pchisq(e95.loglinml2$QvH-e95.loglinml1$QvH,2),2) round(e95.loglinml2$QwH-e95.loglinml1$QwH,2) round(1-pchisq(e95.loglinml2$QwH-e95.loglinml1$QwH,2),2) #Exemplo 9.6 (p.281) / 9.3 (p.269) / 1.6 (p.11): Problema dos defeitos de fibras têxteis e96.TF<-c(28,40,68,5,21,49,1,4,15, 31,70,69,5,12,10,0,1,2) e96.catdata<-readCatdata(TF=e96.TF) e96.X1<-rbind(c(0,0,0,0,0, 0,0,0,0,0,0,0,0, 1), c(0,1,0,0,0, 0,0,0,0,0,0,0,0, 0), c(0,0,1,0,0, 0,0,0,0,0,0,0,0,-1), c(0,0,0,1,0, 0,0,0,0,0,0,0,0, 0), c(0,1,0,1,0, 0,0,0,0,1,0,0,0, 0), c(0,0,1,1,0, 0,0,0,0,0,0,1,0, 0), c(0,0,0,0,1, 0,0,0,0,0,0,0,0,-1), c(0,1,0,0,1, 0,0,0,0,0,1,0,0, 0), c(0,0,1,0,1, 0,0,0,0,0,0,0,1, 1), c(1,0,0,0,0, 0,0,0,0,0,0,0,0,-1), c(1,1,0,0,0, 1,0,0,0,0,0,0,0, 0), c(1,0,1,0,0, 0,1,0,0,0,0,0,0, 1), c(1,0,0,1,0, 0,0,1,0,0,0,0,0, 0), c(1,1,0,1,0, 1,0,1,0,1,0,0,0, 0), c(1,0,1,1,0, 0,1,1,0,0,0,1,0, 0), c(1,0,0,0,1, 0,0,0,1,0,0,0,0, 1), c(1,1,0,0,1, 1,0,0,1,0,1,0,0, 0), c(1,0,1,0,1, 0,1,0,1,0,0,0,1,-1)) e96.loglinml1<-loglinML(e96.catdata,X=e96.X1) #associação uniforme comum e96.loglinml2<-loglinML(e96.catdata,X=e96.X1[,-14]) #sob ausência de interação de 2a.ordem e96.loglinml3<-loglinML(e96.catdata,X=e96.X1[,-c(8,9,14)]) #AC,BC e96.loglinml4<-loglinML(e96.catdata,X=e96.X1[,-(10:14)]) #AB,AC e96.loglinml5<-loglinML(e96.catdata,X=e96.X1[,-c(6,7,14)]) #AB,BC e96.X2<-rbind(c( 1, 1, 0, 1, 0, 1, 1, 1), c( 1, 1, 0, 0, 1, 1, 0, 0), c( 1, 1, 0,-1,-1, 1,-1,-1), c( 1, 0, 1, 1, 0, 0, 1, 0), c( 1, 0, 1, 0, 1, 0, 0, 0), c( 1, 0, 1,-1,-1, 0,-1, 0), c( 1,-1,-1, 1, 0, -1, 1,-1), c( 1,-1,-1, 0, 1, -1, 0, 0), c( 1,-1,-1,-1,-1, -1,-1, 1), c(-1, 1, 0, 1, 0, -1,-1, 1), c(-1, 1, 0, 0, 1, -1, 0, 0), c(-1, 1, 0,-1,-1, -1, 1,-1), c(-1, 0, 1, 1, 0, 0,-1, 0), c(-1, 0, 1, 0, 1, 0, 0, 0), c(-1, 0, 1,-1,-1, 0, 1, 0), c(-1,-1,-1, 1, 0, 1,-1,-1), c(-1,-1,-1, 0, 1, 1, 0, 0), c(-1,-1,-1,-1,-1, 1, 1, 1)) e96.loglinml6<-loglinML(e96.catdata,X=e96.X2) round(e96.loglinml6$QvH-e96.loglinml2$QvH,2) round(1-pchisq(e96.loglinml6$QvH-e96.loglinml2$QvH,5),2) round(e96.loglinml6$QwH-e96.loglinml2$QwH,2) round(1-pchisq(e96.loglinml6$QwH-e96.loglinml2$QwH,5),2) round(exp(e96.loglinml6$beta[8]),2) round(exp(e96.loglinml6$beta[8])*exp(c(-1,1)*qnorm(0.975)*sqrt(e96.loglinml6$Vbeta[8,8])),2) round(exp(-2*e96.loglinml6$beta[7]),2) round(exp(-2*e96.loglinml6$beta[6]),2) round(exp(-2*e96.loglinml6$beta[7])*exp(c(-1,1)*qnorm(0.975)*2*sqrt(e96.loglinml6$Vbeta[7,7])),2) round(exp(-2*e96.loglinml6$beta[6])*exp(c(-1,1)*qnorm(0.975)*2*sqrt(e96.loglinml6$Vbeta[6,6])),2) #Exemplo 9.7 (p.285): Problema da toxicodependência e97.TF<-c(267,577,55, 32, 97,11, 29, 82,6, 13,39,6, 321,622,35, 55,149,16, 49,113,4, 14,26,3, 118,289,45, 21, 53, 9, 14, 35,3, 6,20,0, 135,284,26, 27, 80,13, 20, 36,5, 4,20,5) e97.catdata<-readCatdata(TF=e97.TF) #A,B,C,C,C,D,D,AB,AC,AC,AC,AD,AD,BC,BC,BC,BD,BD,CD,CD,CD,CD,CD,CD, #ABC,ABC,ABC,ABD,ABD,ACD,ACD,ACD,ACD,ACD,ACD,BCD,BCD,BCD,BCD,BCD,BCD e97.X1<-rbind( c(0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,1,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,1,0,0,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,1,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,1,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,1,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,0,0,0,0,1, 0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,1,0,0,0,0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,1,0,0,1,0, 0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0), c(0,1,1,0,0,0,1, 0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0), c(0,1,0,1,0,0,0, 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,0,1,0,1,0, 0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0), c(0,1,0,1,0,0,1, 0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0), c(0,1,0,0,1,0,0, 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,1,0,0,1,1,0, 0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0), c(0,1,0,0,1,0,1, 0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1), c(1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,0,0,0,1,0, 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,0,0,0,0,1, 0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,1,0,0,0,0, 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,1,0,0,1,0, 0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0), c(1,0,1,0,0,0,1, 0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0), c(1,0,0,1,0,0,0, 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,0,1,0,1,0, 0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0), c(1,0,0,1,0,0,1, 0,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0), c(1,0,0,0,1,0,0, 0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,0,0,1,1,0, 0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0), c(1,0,0,0,1,0,1, 0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0), c(1,1,0,0,0,0,0, 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,1,0,0,0,1,0, 1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0, 0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,1,0,0,0,0,1, 1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0, 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0), c(1,1,1,0,0,0,0, 1,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0, 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,1,1,0,0,1,0, 1,1,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0, 1,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0), c(1,1,1,0,0,0,1, 1,1,0,0,0,1,1,0,0,0,1,0,1,0,0,0,0, 1,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0), c(1,1,0,1,0,0,0, 1,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0, 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,1,0,1,0,1,0, 1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,0, 0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0), c(1,1,0,1,0,0,1, 1,0,1,0,0,1,0,1,0,0,1,0,0,0,1,0,0, 0,1,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0), c(1,1,0,0,1,0,0, 1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0, 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(1,1,0,0,1,1,0, 1,0,0,1,1,0,0,0,1,1,0,0,0,0,0,1,0, 0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0), c(1,1,0,0,1,0,1, 1,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0,1, 0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,1)) e97.loglinml1<-loglinML(e97.catdata,X=e97.X1) #(ABC,ABD,ACD,BCD) e97.loglinml2<-loglinML(e97.catdata,X=e97.X1[,1:24]) #(AB,AC,AD,BC,BD,CD) e97.loglinml3<-loglinML(e97.catdata,X=e97.X1[,c(1:7,12:24)]) #(AD,BC,BD,CD) e97.loglinml4<-loglinML(e97.catdata,X=e97.X1[,c(1:7,14:24)]) #(A,BC,BD,CD) e97.loglinml5<-loglinML(e97.catdata,X=e97.X1[,c(1:7,12:13,17:24)]) #(AD,BD,CD) e97.loglinml6<-loglinML(e97.catdata,X=e97.X1[,c(1:7,12:16,19:24)]) #(AD,BC,CD) e97.loglinml7<-loglinML(e97.catdata,X=e97.X1[,c(1:7,12:18)]) #(AD,BC,BD) e97.X2<-rbind( c( 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), c( 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0), c( 1, 1, 1, 0, 0,-1,-1, 1, 1, 0, 0,-1,-1, 1, 0, 0,-1,-1,-1,-1, 0, 0, 0, 0, 1, 0, 0,-1,-1,-1,-1, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0), c( 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0), c( 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0), c( 1, 1, 0, 1, 0,-1,-1, 1, 0, 1, 0,-1,-1, 0, 1, 0,-1,-1, 0, 0,-1,-1, 0, 0, 0, 1, 0,-1,-1, 0, 0,-1,-1, 0, 0, 0, 0,-1,-1, 0, 0), c( 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0), c( 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1), c( 1, 1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,-1, 0, 0, 1,-1,-1, 0, 0, 0, 0,-1,-1, 0, 0, 1,-1,-1, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0,-1,-1), c( 1, 1,-1,-1,-1, 1, 0, 1,-1,-1,-1, 1, 0,-1,-1,-1, 1, 0,-1, 0,-1, 0,-1, 0,-1,-1,-1, 1, 0,-1, 0,-1, 0,-1, 0,-1, 0,-1, 0,-1, 0), c( 1, 1,-1,-1,-1, 0, 1, 1,-1,-1,-1, 0, 1,-1,-1,-1, 0, 1, 0,-1, 0,-1, 0,-1,-1,-1,-1, 0, 1, 0,-1, 0,-1, 0,-1, 0,-1, 0,-1, 0,-1), c( 1, 1,-1,-1,-1,-1,-1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), c( 1,-1, 1, 0, 0, 1, 0,-1, 1, 0, 0, 1, 0,-1, 0, 0,-1, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0,-1, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0), c( 1,-1, 1, 0, 0, 0, 1,-1, 1, 0, 0, 0, 1,-1, 0, 0, 0,-1, 0, 1, 0, 0, 0, 0,-1, 0, 0, 0,-1, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0), c( 1,-1, 1, 0, 0,-1,-1,-1, 1, 0, 0,-1,-1,-1, 0, 0, 1, 1,-1,-1, 0, 0, 0, 0,-1, 0, 0, 1, 1,-1,-1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0), c( 1,-1, 0, 1, 0, 1, 0,-1, 0, 1, 0, 1, 0, 0,-1, 0,-1, 0, 0, 0, 1, 0, 0, 0, 0,-1, 0,-1, 0, 0, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0, 0), c( 1,-1, 0, 1, 0, 0, 1,-1, 0, 1, 0, 0, 1, 0,-1, 0, 0,-1, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0,-1, 0, 0, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0), c( 1,-1, 0, 1, 0,-1,-1,-1, 0, 1, 0,-1,-1, 0,-1, 0, 1, 1, 0, 0,-1,-1, 0, 0, 0,-1, 0, 1, 1, 0, 0,-1,-1, 0, 0, 0, 0, 1, 1, 0, 0), c( 1,-1, 0, 0, 1, 1, 0,-1, 0, 0, 1, 1, 0, 0, 0,-1,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0,-1,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,-1, 0), c( 1,-1, 0, 0, 1, 0, 1,-1, 0, 0, 1, 0, 1, 0, 0,-1, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,-1), c( 1,-1, 0, 0, 1,-1,-1,-1, 0, 0, 1,-1,-1, 0, 0,-1, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0,-1, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0, 1, 1), c( 1,-1,-1,-1,-1, 1, 0,-1,-1,-1,-1, 1, 0, 1, 1, 1,-1, 0,-1, 0,-1, 0,-1, 0, 1, 1, 1,-1, 0,-1, 0,-1, 0,-1, 0, 1, 0, 1, 0, 1, 0), c( 1,-1,-1,-1,-1, 0, 1,-1,-1,-1,-1, 0, 1, 1, 1, 1, 0,-1, 0,-1, 0,-1, 0,-1, 1, 1, 1, 0,-1, 0,-1, 0,-1, 0,-1, 0, 1, 0, 1, 0, 1), c( 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1), c(-1, 1, 1, 0, 0, 1, 0,-1,-1, 0, 0,-1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0,-1, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), c(-1, 1, 1, 0, 0, 0, 1,-1,-1, 0, 0, 0,-1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0,-1, 0, 0, 0,-1, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0), c(-1, 1, 1, 0, 0,-1,-1,-1,-1, 0, 0, 1, 1, 1, 0, 0,-1,-1,-1,-1, 0, 0, 0, 0,-1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0), c(-1, 1, 0, 1, 0, 1, 0,-1, 0,-1, 0,-1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0,-1, 0,-1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0), c(-1, 1, 0, 1, 0, 0, 1,-1, 0,-1, 0, 0,-1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0,-1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0), c(-1, 1, 0, 1, 0,-1,-1,-1, 0,-1, 0, 1, 1, 0, 1, 0,-1,-1, 0, 0,-1,-1, 0, 0, 0,-1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0), c(-1, 1, 0, 0, 1, 1, 0,-1, 0, 0,-1,-1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,-1,-1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 1, 0), c(-1, 1, 0, 0, 1, 0, 1,-1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 1), c(-1, 1, 0, 0, 1,-1,-1,-1, 0, 0,-1, 1, 1, 0, 0, 1,-1,-1, 0, 0, 0, 0,-1,-1, 0, 0,-1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1), c(-1, 1,-1,-1,-1, 1, 0,-1, 1, 1, 1,-1, 0,-1,-1,-1, 1, 0,-1, 0,-1, 0,-1, 0, 1, 1, 1,-1, 0, 1, 0, 1, 0, 1, 0,-1, 0,-1, 0,-1, 0), c(-1, 1,-1,-1,-1, 0, 1,-1, 1, 1, 1, 0,-1,-1,-1,-1, 0, 1, 0,-1, 0,-1, 0,-1, 1, 1, 1, 0,-1, 0, 1, 0, 1, 0, 1, 0,-1, 0,-1, 0,-1), c(-1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1), c(-1,-1, 1, 0, 0, 1, 0, 1,-1, 0, 0,-1, 0,-1, 0, 0,-1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,-1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0), c(-1,-1, 1, 0, 0, 0, 1, 1,-1, 0, 0, 0,-1,-1, 0, 0, 0,-1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0,-1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0), c(-1,-1, 1, 0, 0,-1,-1, 1,-1, 0, 0, 1, 1,-1, 0, 0, 1, 1,-1,-1, 0, 0, 0, 0, 1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0), c(-1,-1, 0, 1, 0, 1, 0, 1, 0,-1, 0,-1, 0, 0,-1, 0,-1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0,-1, 0, 0, 0), c(-1,-1, 0, 1, 0, 0, 1, 1, 0,-1, 0, 0,-1, 0,-1, 0, 0,-1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0,-1, 0, 0), c(-1,-1, 0, 1, 0,-1,-1, 1, 0,-1, 0, 1, 1, 0,-1, 0, 1, 1, 0, 0,-1,-1, 0, 0, 0, 1, 0,-1,-1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0), c(-1,-1, 0, 0, 1, 1, 0, 1, 0, 0,-1,-1, 0, 0, 0,-1,-1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0,-1, 0), c(-1,-1, 0, 0, 1, 0, 1, 1, 0, 0,-1, 0,-1, 0, 0,-1, 0,-1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0,-1), c(-1,-1, 0, 0, 1,-1,-1, 1, 0, 0,-1, 1, 1, 0, 0,-1, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0, 1,-1,-1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1), c(-1,-1,-1,-1,-1, 1, 0, 1, 1, 1, 1,-1, 0, 1, 1, 1,-1, 0,-1, 0,-1, 0,-1, 0,-1,-1,-1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0), c(-1,-1,-1,-1,-1, 0, 1, 1, 1, 1, 1, 0,-1, 1, 1, 1, 0,-1, 0,-1, 0,-1, 0,-1,-1,-1,-1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), c(-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1)) e97.loglinml8<-loglinML(e97.catdata,X=e97.X2[,c(1:7,12:24)]) #(AD,BC,BD,CD) #Exemplo 9.10 (p.304) / 1.4 (p.5): Problema dos grafiteiros e910.TF<-rbind(c( 3,28,39,19, 1,75, 79), c( 0,14,24,25,30,94, 69), c(173,69,19,28, 1,45,165), c( 3,34, 3,47,90,36,136)) e910.catdata<-readCatdata(TF=e910.TF) e910.A<-kronecker(diag(4),cbind(diag(6),rep(-1,6))) e910.XL<-kronecker(rbind(c(1,1,1),c(1,1,-1),c(1,-1,1),c(1,-1,-1)),diag(6)) e910.loglinml<-loglinML(e910.catdata,A=e910.A,XL=e910.XL) round(exp(2*e910.loglinml$beta[7:18]),2) round(exp(2*e910.loglinml$beta[7:18]-qnorm(0.975)*2*sqrt(diag(e910.loglinml$Vbeta))[7:18]),2) round(exp(2*e910.loglinml$beta[7:18]+qnorm(0.975)*2*sqrt(diag(e910.loglinml$Vbeta))[7:18]),2) #Exemplo 9.11 (p.307) / 1.5 (p.5): Problema do uso do fio dental e911.TF<-rbind(c(19,5,4, 2), c( 5,8,0,17), c(11,6,7, 6), c( 2,5,1,22)) e911.catdata<-readCatdata(TF=e911.TF) e911.A<-kronecker(diag(4),t(c(1,-1,-1,1))) e911.XL<-cbind(rep(1,4),c(1,1,-1,-1),c(1,-1,1,-1)) e911.loglinml1<-loglinML(e911.catdata,A=e911.A,XL=e911.XL) #ABC,ABD,ACD,BCD e911.loglinml2<-loglinML(e911.catdata,A=e911.A,XL=e911.XL[,-2]) #ABC,ABD,BCD e911.loglinml3<-loglinML(e911.catdata,A=e911.A,XL=e911.XL[,-3]) #ABC,ABD,ACD e911.loglinml4<-loglinML(e911.catdata,A=e911.A,XL=e911.XL[,1]) #ABC,ABD,CD #A,B,C,D, AB,AC,AD, BC,BD,CD,BCD e911.X<-rbind(c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), c( 1, 1, 1,-1, 1, 1,-1, 1,-1,-1,-1), c( 1, 1,-1, 1, 1,-1, 1, -1, 1,-1,-1), c( 1, 1,-1,-1, 1,-1,-1, -1,-1, 1, 1), c( 1,-1, 1, 1, -1, 1, 1, -1,-1, 1,-1), c( 1,-1, 1,-1, -1, 1,-1, -1, 1,-1, 1), c( 1,-1,-1, 1, -1,-1, 1, 1,-1,-1, 1), c( 1,-1,-1,-1, -1,-1,-1, 1, 1, 1,-1), c(-1, 1, 1, 1, -1,-1,-1, 1, 1, 1, 1), c(-1, 1, 1,-1, -1,-1, 1, 1,-1,-1,-1), c(-1, 1,-1, 1, -1, 1,-1, -1, 1,-1,-1), c(-1, 1,-1,-1, -1, 1, 1, -1,-1, 1, 1), c(-1,-1, 1, 1, 1,-1,-1, -1,-1, 1,-1), c(-1,-1, 1,-1, 1,-1, 1, -1, 1,-1, 1), c(-1,-1,-1, 1, 1, 1,-1, 1,-1,-1, 1), c(-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1)) #sob multinom.ao invés de prod.de multinom. #C,D, AC,AD, BC,BD,CD,BCD e911.X<-rbind(c( 1, 1, 1, 1, 1, 1, 1, 1), c( 1,-1, 1,-1, 1,-1,-1,-1), c(-1, 1, -1, 1, -1, 1,-1,-1), c(-1,-1, -1,-1, -1,-1, 1, 1), c( 1, 1, 1, 1, -1,-1, 1,-1), c( 1,-1, 1,-1, -1, 1,-1, 1), c(-1, 1, -1, 1, 1,-1,-1, 1), c(-1,-1, -1,-1, 1, 1, 1,-1), c( 1, 1, -1,-1, 1, 1, 1, 1), c( 1,-1, -1, 1, 1,-1,-1,-1), c(-1, 1, 1,-1, -1, 1,-1,-1), c(-1,-1, 1, 1, -1,-1, 1, 1), c( 1, 1, -1,-1, -1,-1, 1,-1), c( 1,-1, -1, 1, -1, 1,-1, 1), c(-1, 1, 1,-1, 1,-1,-1, 1), c(-1,-1, 1, 1, 1, 1, 1,-1)) e911.loglinml5<-loglinML(e911.catdata,X=e911.X) #AB,AC,AD,BCD e911.loglinml6<-loglinML(e911.catdata,X=e911.X[,-4]) #AB,AC,BCD e911.loglinml7<-loglinML(e911.catdata,X=e911.X[,-(3:4)]) #AB,BCD e911.loglinml8<-loglinML(e911.catdata,X=e911.X[,-c(4,8)]) #AB,AC,BC,BD,CD e911.loglinml9<-loglinML(e911.catdata,X=e911.X[,-8]) #AB,AC,AD,BC,BD,CD round(e911.loglinml9$QvH-e911.loglinml5$QvH,2) #BCD=0|H_2 round(1-pchisq(e911.loglinml9$QvH-e911.loglinml5$QvH,1),2) round(e911.loglinml6$QvH-e911.loglinml5$QvH,2) #AD=0|H_2 round(1-pchisq(e911.loglinml6$QvH-e911.loglinml5$QvH,1),2) round(e911.loglinml8$QvH-e911.loglinml6$QvH,2) #BCD=0|H_3 round(1-pchisq(e911.loglinml8$QvH-e911.loglinml6$QvH,1),2) round(exp(4*(e911.loglinml6$beta[4:6]+e911.loglinml6$beta[7])),2) #RPC 5-8 anos: AD, AC, AB round(exp(4*(e911.loglinml6$beta[4:6]-e911.loglinml6$beta[7])),2) #RPC 9-12 anos: AD, AC, AB round(exp(4*e911.loglinml6$beta[3]),2) #RPC BD #Exemplo 9.12 (pp.316, 317, 319, 324, 326): Problema da aterosclerose coronariana e912.TF<-c(31,17,42,27,55,42,94,104,80,112,70,130,74,188,68,314) e912.catdata<-readCatdata(TF=e912.TF) #A,B,C,D, AB,AC,AD,BC,BD,CD, ABC,ABD,ACD,BCD e912.X<-rbind(c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), c( 1, 1, 1,-1, 1, 1,-1, 1,-1,-1, 1,-1,-1,-1), c( 1, 1,-1, 1, 1,-1, 1,-1, 1,-1, -1, 1,-1,-1), c( 1, 1,-1,-1, 1,-1,-1,-1,-1, 1, -1,-1, 1, 1), c( 1,-1, 1, 1, -1, 1, 1,-1,-1, 1, -1,-1, 1,-1), c( 1,-1, 1,-1, -1, 1,-1,-1, 1,-1, -1, 1,-1, 1), c( 1,-1,-1, 1, -1,-1, 1, 1,-1,-1, 1,-1,-1, 1), c( 1,-1,-1,-1, -1,-1,-1, 1, 1, 1, 1, 1, 1,-1), c(-1, 1, 1, 1, -1,-1,-1, 1, 1, 1, -1,-1,-1, 1), c(-1, 1, 1,-1, -1,-1, 1, 1,-1,-1, -1, 1, 1,-1), c(-1, 1,-1, 1, -1, 1,-1,-1, 1,-1, 1,-1, 1,-1), c(-1, 1,-1,-1, -1, 1, 1,-1,-1, 1, 1, 1,-1, 1), c(-1,-1, 1, 1, 1,-1,-1,-1,-1, 1, 1, 1,-1,-1), c(-1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1, 1), c(-1,-1,-1, 1, 1, 1,-1, 1,-1,-1, -1, 1, 1, 1), c(-1,-1,-1,-1, 1, 1, 1, 1, 1, 1, -1,-1,-1,-1)) e912.loglinml1<-loglinML(e912.catdata,X=e912.X) e912.loglinml2<-loglinML(e912.catdata,X=e912.X[,1:10]) e912.loglinml3<-loglinML(e912.catdata,X=e912.X[,1:4]) round(e912.loglinml2$QvH-e912.loglinml1$QvH,2) round(1-pchisq(e912.loglinml2$QvH-e912.loglinml1$QvH,e912.loglinml2$glH-e912.loglinml1$glH),2) round(e912.loglinml3$QvH-e912.loglinml2$QvH,2) round(1-pchisq(e912.loglinml3$QvH-e912.loglinml2$QvH,e912.loglinml3$glH-e912.loglinml2$glH),2) #continuação p.317 e912.loglinml4<-loglinML(e912.catdata,X=cbind(e912.X,e912.X[,1]*e912.X[,2]*e912.X[,3]*e912.X[,4])) e912.loglinml5<-loglinML(e912.catdata,X=e912.X[,-14]) e912.loglinml6<-loglinML(e912.catdata,X=e912.X[,-13]) e912.loglinml7<-loglinML(e912.catdata,X=e912.X[,-12]) e912.loglinml8<-loglinML(e912.catdata,X=e912.X[,-11]) e912.loglinml9<-loglinML(e912.catdata,X=e912.X[,-c(10:14)]) e912.loglinml10<-loglinML(e912.catdata,X=e912.X[,-c(9,11:14)]) e912.loglinml11<-loglinML(e912.catdata,X=e912.X[,-c(8,11:14)]) e912.loglinml12<-loglinML(e912.catdata,X=e912.X[,-c(7,11:14)]) e912.loglinml13<-loglinML(e912.catdata,X=e912.X[,-c(6,11:14)]) e912.loglinml14<-loglinML(e912.catdata,X=e912.X[,-c(5,11:14)]) round(c(e912.loglinml14$QvH-e912.loglinml2$QvH, e912.loglinml13$QvH-e912.loglinml2$QvH, e912.loglinml12$QvH-e912.loglinml2$QvH, e912.loglinml11$QvH-e912.loglinml2$QvH, e912.loglinml10$QvH-e912.loglinml2$QvH, e912.loglinml9$QvH-e912.loglinml2$QvH, e912.loglinml8$QvH-e912.loglinml1$QvH, e912.loglinml7$QvH-e912.loglinml1$QvH, e912.loglinml6$QvH-e912.loglinml1$QvH, e912.loglinml5$QvH-e912.loglinml1$QvH, e912.loglinml1$QvH),2) #parcial round(1-pchisq(c(e912.loglinml14$QvH-e912.loglinml2$QvH, e912.loglinml13$QvH-e912.loglinml2$QvH, e912.loglinml12$QvH-e912.loglinml2$QvH, e912.loglinml11$QvH-e912.loglinml2$QvH, e912.loglinml10$QvH-e912.loglinml2$QvH, e912.loglinml9$QvH-e912.loglinml2$QvH, e912.loglinml8$QvH-e912.loglinml1$QvH, e912.loglinml7$QvH-e912.loglinml1$QvH, e912.loglinml6$QvH-e912.loglinml1$QvH, e912.loglinml5$QvH-e912.loglinml1$QvH, e912.loglinml1$QvH),1),2) #parcial e912.loglinml15<-loglinML(e912.catdata,X=e912.X[,c(1:10,14)]) e912.loglinml16<-loglinML(e912.catdata,X=e912.X[,c(1:10,13)]) e912.loglinml17<-loglinML(e912.catdata,X=e912.X[,c(1:10,12)]) e912.loglinml18<-loglinML(e912.catdata,X=e912.X[,c(1:11)]) #ABC,AD,BD,CD: modelo final e912.loglinml19<-loglinML(e912.catdata,X=e912.X[,c(1:4,10)]) e912.loglinml20<-loglinML(e912.catdata,X=e912.X[,c(1:4,9)]) e912.loglinml21<-loglinML(e912.catdata,X=e912.X[,c(1:4,8)]) e912.loglinml22<-loglinML(e912.catdata,X=e912.X[,c(1:4,7)]) e912.loglinml23<-loglinML(e912.catdata,X=e912.X[,c(1:4,6)]) e912.loglinml24<-loglinML(e912.catdata,X=e912.X[,c(1:4,5)]) round(c(e912.loglinml3$QvH-e912.loglinml24$QvH, e912.loglinml3$QvH-e912.loglinml23$QvH, e912.loglinml3$QvH-e912.loglinml22$QvH, e912.loglinml3$QvH-e912.loglinml21$QvH, e912.loglinml3$QvH-e912.loglinml20$QvH, e912.loglinml3$QvH-e912.loglinml19$QvH, e912.loglinml2$QvH-e912.loglinml18$QvH, e912.loglinml2$QvH-e912.loglinml17$QvH, e912.loglinml2$QvH-e912.loglinml16$QvH, e912.loglinml2$QvH-e912.loglinml15$QvH, e912.loglinml1$QvH),2) #marginal round(1-pchisq(c(e912.loglinml3$QvH-e912.loglinml24$QvH, e912.loglinml3$QvH-e912.loglinml23$QvH, e912.loglinml3$QvH-e912.loglinml22$QvH, e912.loglinml3$QvH-e912.loglinml21$QvH, e912.loglinml3$QvH-e912.loglinml20$QvH, e912.loglinml3$QvH-e912.loglinml19$QvH, e912.loglinml2$QvH-e912.loglinml18$QvH, e912.loglinml2$QvH-e912.loglinml17$QvH, e912.loglinml2$QvH-e912.loglinml16$QvH, e912.loglinml2$QvH-e912.loglinml15$QvH, e912.loglinml1$QvH),1),2) #marginal #continuação p.326 e9122.TF<-rbind(c(31,17),c(42,27),c(55,42),c(94,104),c(80,112),c(70,130),c(74,188),c(68,314)) e9122.catdata<-readCatdata(TF=e9122.TF) #ABC,AD,BD,CD e9122.XL<-rbind(c(1,1,1,1), c(1,1,1,0), c(1,1,0,1), c(1,1,0,0), c(1,0,1,1), c(1,0,1,0), c(1,0,0,1), c(1,0,0,0)) e9122.loglinml<-loglinML(e9122.catdata,XL=e9122.XL) round(exp(e9122.loglinml$beta[2:4]),2) round(exp(e9122.loglinml$beta[2:4]-qnorm(0.975)*sqrt(diag(e9122.loglinml$Vbeta))[2:4]),2) round(exp(e9122.loglinml$beta[2:4]+qnorm(0.975)*sqrt(diag(e9122.loglinml$Vbeta))[2:4]),2) round(kronecker(diag(8),t(c(1,0)))%*%e912.loglinml18$thetaH/kronecker(diag(8),t(c(1,1)))%*% e912.loglinml18$thetaH,3) #o mesmo que se encontra nas prob.estim. de summary(e9122.loglinml) #Exemplo 10.1 (p.348) / 6.1 (p.149): Problema da intoxicação de besouros e101.TF<-cbind(c(6,13,18,28,52,53,61,60),c(59,60,62,56,63,59,62,60)) e101.TF[,2]<-e101.TF[,2]-e101.TF[,1] e101.catdata<-readCatdata(TF=e101.TF) e101.XL<-cbind(rep(1,8),c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)) e101.loglinml<-loglinML(e101.catdata,XL=e101.XL) round(exp(e101.loglinml$beta[2]*log10(2)),2) round(exp((e101.loglinml$beta[2]+c(-1,1)*qnorm(0.975)*sqrt(e101.loglinml$Vbeta[2,2]))*log10(2)),2) #Exemplo 10.2 (p.349) / 6.3 (p.156) / 1.2 (p.4): Problema do risco de cárie dentária e102.TF<-c(11,5,0,14,34,7,2,13,11) e102.catdata<-readCatdata(TF=e102.TF) e102.B<-rbind(c(1,-1,0),c(0,1,-1)) e102.A<-kronecker(e102.B,e102.B) e102.XL<-rep(1,4) e102.loglinml1<-loglinML(e102.catdata,A=e102.A,XL=e102.XL) e102.X<-rbind(c(0,0,0,0, 1), c(1,0,0,0, 0), c(0,1,0,0,-1), c(0,0,1,0, 0), c(1,0,1,0, 0), c(0,1,1,0, 0), c(0,0,0,1,-1), c(1,0,0,1, 0), c(0,1,0,1, 1)) e102.loglinml2<-loglinML(e102.catdata,X=e102.X) round(e102.loglinml1$beta,2) round(e102.loglinml1$beta+c(-1,1)*qnorm(0.975)*sqrt(e102.loglinml1$Vbeta),2) round(exp(e102.loglinml1$beta),2) round(exp(e102.loglinml1$beta+c(-1,1)*qnorm(0.975)*sqrt(e102.loglinml1$Vbeta)),2) #Exemplo 10.3 (p.349) / 6.4 (p.157) / 1.5 (p.5): Problema do uso do fio dental e103.TF<-rbind(c(19,5,4, 2), c( 5,8,0,17), c(11,6,7, 6), c( 2,5,1,22)) e103.catdata<-readCatdata(TF=e103.TF) e103.A<-kronecker(diag(4),t(c(1,-1,-1,1))) e103.XL<-cbind(rep(1,4),c(1,1,0,0),c(1,0,1,0)) e103.loglinml<-loglinML(e103.catdata,A=e103.A,XL=e103.XL) round(exp(e103.loglinml$beta),2) round(exp(e103.loglinml$beta-qnorm(0.975)*sqrt(diag(e103.loglinml$Vbeta))),2) round(exp(e103.loglinml$beta+qnorm(0.975)*sqrt(diag(e103.loglinml$Vbeta))),2) round(exp(t(c(0,1,1))%*%e103.loglinml$beta),2) round(exp(t(c(0,1,1))%*%e103.loglinml$beta+ c(-1,1)*qnorm(0.975)*sqrt(t(c(0,1,1))%*%e103.loglinml$Vbeta%*%c(0,1,1))),2) #Exemplo 10.4 (p.351) / 6.5 (p.160): Problema da complicação pulmonar e104.catdata<-readCatdata(TF=cbind(c(737,243,39),c(48,74,21))) e104.A<-rbind(c(0,-1,0,1,0,0),c(0,0,0,-1,0,1)) e104.loglinml<-loglinML(e104.catdata,A=e104.A,XL=c(1,1)) #A rotina loglinML não ajusta este tipo de modelo log-linear generalizado, #portanto utiliza-se a funlinWLS para o ajuste por MQG e104.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e104.catdata,A1=e104.A,X=c(1,1)) e104.loglinwls2<-funlinWLS(model=c("lin","log"),obj=e104.catdata,A1=e104.A,X=c(2,1)) exp(e104.loglinwls2$beta) exp(e104.loglinwls2$beta+c(-1,1)*qnorm(0.975)*sqrt(e104.loglinwls2$Vbeta)) exp(2*e104.loglinwls2$beta) exp(2*(e104.loglinwls2$beta+c(-1,1)*qnorm(0.975)*sqrt(e104.loglinwls2$Vbeta))) #Exemplo 10.5 (p.353) / 6.6 (p.163) / 1.3 (p.4): Problema do peso de recém-nascidos e105.TF<-rbind(c( 2, 11, 31),c( 5, 24, 95), c( 3, 32, 91),c( 11, 57, 238), c( 15, 58,134),c( 25,105, 445), c(130,362,695),c(231,694,2485), c( 94,225,340),c(105,339,1053)) e105.catdata<-readCatdata(TF=e105.TF) e105.A1<-kronecker(diag(10),rbind(c(1,0,0),c(0,1,1),c(0,1,0),c(0,0,1))) e105.A2<-kronecker(diag(10),kronecker(diag(2),t(c(1,-1)))) e105.X1<-kronecker(rbind(c(1, 1, 0, 0, 0, 1), c(1, 1, 0, 0, 0,-1), c(1, 0, 1, 0, 0, 1), c(1, 0, 1, 0, 0,-1), c(1, 0, 0, 1, 0, 1), c(1, 0, 0, 1, 0,-1), c(1, 0, 0, 0, 1, 1), c(1, 0, 0, 0, 1,-1), c(1,-1,-1,-1,-1, 1), c(1,-1,-1,-1,-1,-1)),diag(2)) e105.X2<-kronecker(rbind(c(1,-2, 1), c(1,-2,-1), c(1,-1, 1), c(1,-1,-1), c(1, 0, 1), c(1, 0,-1), c(1, 1, 1), c(1, 1,-1), c(1, 2, 1), c(1, 2,-1)),diag(2)) e105.funlinwls1<-funlinWLS(model=c("lin","log","lin"),obj=e105.catdata, A1=e105.A1,A2=e105.A2,X=e105.X1) e105.funlinwls2<-funlinWLS(model=c("lin","log","lin"),obj=e105.catdata, A1=e105.A1,A2=e105.A2,X=e105.X2) e105.TF2<-matrix(c(e105.A1%*%c(t(e105.TF))),20,2,byrow=T) e105.catdata2<-readCatdata(TF=e105.TF2) e105.A3<-kronecker(diag(20),t(c(1,-1))) e105.loglinml1<-loglinML(e105.catdata2,A=e105.A3,XL=e105.X1) e105.loglinml2<-loglinML(e105.catdata2,A=e105.A3,XL=e105.X2) round(exp(-e105.loglinml2$beta[3:4]),2) round(exp(-e105.loglinml2$beta[3:4]-qnorm(0.975)*sqrt(diag(e105.loglinml2$Vbeta))[3:4]),2) round(exp(-e105.loglinml2$beta[3:4]+qnorm(0.975)*sqrt(diag(e105.loglinml2$Vbeta))[3:4]),2) round(exp(e105.loglinml2$beta[5:6]),2) round(exp(e105.loglinml2$beta[5:6]-qnorm(0.975)*sqrt(diag(e105.loglinml2$Vbeta))[5:6]),2) round(exp(e105.loglinml2$beta[5:6]+qnorm(0.975)*sqrt(diag(e105.loglinml2$Vbeta))[5:6]),2) #Exemplo 10.6 (p.357) / 6.7 (p.167) / 1.3 (p.4): Problema do peso de recém-nascidos e106.TF<-rbind(c( 2, 11, 31),c( 5, 24, 95), c( 3, 32, 91),c( 11, 57, 238), c( 15, 58,134),c( 25,105, 445), c(130,362,695),c(231,694,2485), c( 94,225,340),c(105,339,1053)) e106.catdata<-readCatdata(TF=e106.TF) e106.A1<-kronecker(diag(10),rbind(c(1,0,0),c(0,1,1),c(1,1,0),c(0,0,1))) e106.A2<-kronecker(diag(10),kronecker(diag(2),t(c(1,-1)))) e106.X2<-rbind(c(1,0,-2, 0, 1), c(0,1,-2, 0, 1), c(1,0, 0,-2,-1), c(0,1, 0,-2,-1), c(1,0,-1, 0, 1), c(0,1,-1, 0, 1), c(1,0, 0,-1,-1), c(0,1, 0,-1,-1), c(1,0, 0, 0, 1), c(0,1, 0, 0, 1), c(1,0, 0, 0,-1), c(0,1, 0, 0,-1), c(1,0, 1, 0, 1), c(0,1, 1, 0, 1), c(1,0, 0, 1,-1), c(0,1, 0, 1,-1), c(1,0, 2, 0, 1), c(0,1, 2, 0, 1), c(1,0, 0, 2,-1), c(0,1, 0, 2,-1)) e106.X1<-cbind(e106.X2[,1:2],e106.X2[,3]+e106.X2[,4],e106.X2[,5]) #Estes resultados (MQG) não são exatamente iguais ao do livro (MV) e106.funlinwls1<-funlinWLS(model=c("lin","log","lin"),obj=e106.catdata, A1=e106.A1,A2=e106.A2,X=e106.X1) e106.funlinwls2<-funlinWLS(model=c("lin","log","lin"),obj=e106.catdata, A1=e106.A1,A2=e106.A2,X=e106.X2) round(e106.funlinwls1$QwH-e106.funlinwls2$QwH,2) round(1-pchisq(e106.funlinwls1$QwH-e106.funlinwls2$QwH,1),2) #Exemplo 10.8 (pp.360, 367): Problema da susceptibilidade a malária cerebral e108.TF<-rbind(c(35,10), c(25,23), c(27,21), c( 9,40)) e108.catdata<-readCatdata(TF=e108.TF) e108.X<-rbind(c(1,0,0), c(1,0,1), c(1,1,0), c(1,1,1)) e108.linml<-linML(e108.catdata,X=e108.X) #aditivo e108.A<-kronecker(diag(4),t(c(1,0))) #Este resultado (MQG) não é exatamente igual ao do livro (MV) e108.funlinwls<-funlinWLS(model=c("lin","log"),obj=e108.catdata, A1=e108.A,X=e108.X) #multiplicativo e108.loglinml<-loglinML(e108.catdata,XL=e108.X) #logístico #Exemplo 11.1 (p.376) / 7.1 (p.202) / 1.7 (p.11): Problema do grupo sanguíneo ABO e111.TF<-cbind(4219,890,313,4578) e111.catdata<-readCatdata(TF=e111.TF) e111.A1<-rbind(diag(4), c(0,0,0,0.5)) e111.A2<-rbind(c(1,0,0 ,1,0), c(0,1,0 ,1,0), c(0,0,0.5,0,1)) e111.A3<-cbind(1,1,-2) e111.funlinwls<-funlinWLS(model=c("lin","log","lin","exp","lin","log"), obj=e111.catdata,A1=e111.A1,A2=e111.A2,A3=e111.A3,X=1) e111.TF<-rbind(c(4219,890,313,4578), c( 96, 18, 5, 181), c( 214, 39, 13, 298)) e111.catdata<-readCatdata(TF=e111.TF) e111.A1<-kronecker(diag(3),rbind(diag(4), c(0,0,0,0.5))) e111.A2<-kronecker(diag(3),rbind(c(1,0,0 ,1,0), c(0,1,0 ,1,0), c(0,0,0.5,0,1))) e111.A3<-kronecker(diag(3),cbind(1,1,-2)) e111.funlinwls<-funlinWLS(model=c("lin","log","lin","exp","lin","log"), obj=e111.catdata,A1=e111.A1,A2=e111.A2,A3=e111.A3,X=diag(3)) waldTest(e111.funlinwls,diag(3)) #Exemplo 11.2 (p.382) / 8.1 (p.228) / 3.1 (p.47): Problema da intenção de voto e112.TF<-c(192,1,5,2,146,5,11,12,71) e112.catdata<-readCatdata(TF=e112.TF) e112.U<-rbind(c(0,-1, 0,1,0, 0,0,0), c(0, 0,-1,0,0, 0,1,0), c(0, 0, 0,0,0,-1,0,1)) e112.X<-rbind(c(1,0,0,0,0), c(0,1,0,0,0), c(0,0,1,0,0), c(0,1,0,0,0), c(0,0,0,1,0), c(0,0,0,0,1), c(0,0,1,0,0), c(0,0,0,0,1)) e112.linwls1<-funlinWLS(model="lin",obj=e112.catdata,U=e112.U) #simetria e112.linwls2<-funlinWLS(model="lin",obj=e112.catdata,X=e112.X) #simetria #Exemplo 11.3 (p.383) / 8.2 (p.233) / 3.2 (p.49) / 1.2 (p.4): Problema do risco de cárie dentária e113.TF<-c(11,5,0,14,34,7,2,13,11) e113.catdata<-readCatdata(TF=e113.TF) e113.U<-rbind(c(0, 1,1,-1,0,0,-1, 0), c(0,-1,0, 1,0,1, 0,-1)) e113.X<-rbind(c(1, 0, 0,0,0,0), c(0, 1, 0,0,0,0), c(0,-1, 1,0,1,0), c(0, 0, 1,0,0,0), c(0, 0, 0,1,0,0), c(0, 1,-1,0,0,1), c(0, 0, 0,0,1,0), c(0, 0, 0,0,0,1)) e113.linwls1<-funlinWLS(model="lin",obj=e113.catdata,U=e113.U) #homogeneidade marginal (HM) e113.linwls2<-funlinWLS(model="lin",obj=e113.catdata,X=e113.X) #homogeneidade marginal (HM) e113.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) ) e113.U2<-rbind(c(1,0,-1, 0), c(0,1, 0,-1)) e113.X2<-rbind(c(1,0), c(0,1), c(1,0), c(0,1) ) e113.linwls3<-funlinWLS(model="lin",obj=e113.catdata,A1=e113.A,U=e113.U2) #HM e113.linwls4<-funlinWLS(model="lin",obj=e113.catdata,A1=e113.A,X=e113.X2) #HM #Exemplo 11.4 (p.384) / 8.3 (p.236) / 3.3 (p.50) / 1.9 (p.12): Problema do tamanho da ninhada e114.TF<-rbind(c(10,21, 96,23), c( 4, 6, 28, 8), c( 9, 7, 58, 7), c( 8,19, 44, 1), c( 5,17, 56, 1), c( 1, 5, 20, 2), c(22,95,103, 4), c(18,49, 62, 0), c( 4,12, 16, 2)) e114.catdata<-readCatdata(TF=e114.TF) e114.A<-kronecker(diag(9),t(c(0,1,2,3))) e114.X<-rbind(c(1, 1, 0, 1, 0, 1, 0, 0, 0), c(1, 1, 0, 0, 1, 0, 1, 0, 0), c(1, 1, 0,-1,-1, -1,-1, 0, 0), c(1, 0, 1, 1, 0, 0, 0, 1, 0), c(1, 0, 1, 0, 1, 0, 0, 0, 1), c(1, 0, 1,-1,-1, 0, 0,-1,-1), c(1,-1,-1, 1, 0, -1, 0,-1, 0), c(1,-1,-1, 0, 1, 0,-1, 0,-1), c(1,-1,-1,-1,-1, 1, 1, 1, 1)) e114.linwls1<-funlinWLS(model="lin",obj=e114.catdata,A1=e114.A,X=e114.X) waldTest(e114.linwls1,cbind(matrix(0,4,5),diag(4))) waldTest(e114.linwls1,cbind(rep(0,2),diag(2),matrix(0,2,6))) waldTest(e114.linwls1,cbind(matrix(0,2,3),diag(2),matrix(0,2,4))) e114.linwls2<-funlinWLS(model="lin",obj=e114.catdata,A1=e114.A,X=e114.X[,1:5]) waldTest(e114.linwls2,cbind(rep(0,2),diag(2),matrix(0,2,2))) waldTest(e114.linwls2,cbind(matrix(0,2,3),diag(2))) e114.linwls3<-funlinWLS(model="lin",obj=e114.catdata,A1=e114.A,X=e114.X[,1:3]) #Exemplo 11.5 (p.388) / 9.1 (p.263): Problema da anemia e115.TF<-c(3,25,32,68) e115.catdata<-readCatdata(TF=e115.TF) e115.U<-c(1,-1,-1,1) e115.X<-rbind(c(0,0), c(0,1), c(1,0), c(1,1)) e115.X2<-rbind(c(0,0,0), c(0,1,0), c(1,0,0), c(1,1,1)) e115.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e115.catdata,U=e115.U) #independência e115.loglinwls2<-funlinWLS(model=c("lin","log"),obj=e115.catdata,X=e115.X) #independência e115.loglinwls3<-funlinWLS(model=c("lin","log"),obj=e115.catdata,X=e115.X2) #modelo saturado e115.loglinwls4<-funlinWLS(model=c("lin","log"),obj=e115.catdata,A1=c(1,-1,-1,1), XL=1) #modelo saturado round(e115.loglinwls4$beta+c(-1,1)*qnorm(0.975)*sqrt(e115.loglinwls4$Vbeta),3) round(exp(e115.loglinwls4$beta),3) #razão de chances round(exp(e115.loglinwls4$beta+c(-1,1)*qnorm(0.975)*sqrt(e115.loglinwls4$Vbeta)),3) #Exemplo 11.6 (p.388) / 9.2 (p.267): Problema da acuidade visual e116.TF<-c(1520,266,124,66, 234,1512,432,78, 117,362,1772,205, 36,82,179,492) e116.catdata<-readCatdata(TF=e116.TF) e116.X1<-rbind(c(1,0,0,0,0,0,0,0,0), c(0,1,0,0,0,0,0,0,0), c(0,0,1,0,0,0,0,0,0), c(0,0,0,1,0,0,0,0,0), c(0,1,0,0,0,0,0,0,0), c(0,0,0,0,1,0,0,0,0), c(0,0,0,0,0,1,0,0,0), c(0,0,0,0,0,0,1,0,0), c(0,0,1,0,0,0,0,0,0), c(0,0,0,0,0,1,0,0,0), c(0,0,0,0,0,0,0,1,0), c(0,0,0,0,0,0,0,0,1), c(0,0,0,1,0,0,0,0,0), c(0,0,0,0,0,0,1,0,0), c(0,0,0,0,0,0,0,0,1)) e116.linwls1<-funlinWLS(model="lin",obj=e116.catdata,X=e116.X1) #simetria em formulação linear e116.A1<-rbind(c(0,1,0,0, -1,0,0,0, 0, 0,0,0, 0, 0, 0,0), c(0,0,1,0, 0,0,0,0, -1, 0,0,0, 0, 0, 0,0), c(0,0,0,1, 0,0,0,0, 0, 0,0,0, -1, 0, 0,0), c(0,0,0,0, 0,0,1,0, 0,-1,0,0, 0, 0, 0,0), c(0,0,0,0, 0,0,0,1, 0, 0,0,0, 0,-1, 0,0), c(0,0,0,0, 0,0,0,0, 0, 0,0,1, 0, 0,-1,0)) e116.linwls2<-funlinWLS(model="lin",obj=e116.catdata,U=e116.A1[,1:15]) #simetria em form.linear #u_1,u_2,u_3, u_{11},u_{12},u_{13},u_{22},u_{23},u_{33} e116.X2<-rbind(c( 2, 0, 0, 1, 0, 0, 0, 0, 0), c( 1, 1, 0, 0, 1, 0, 0, 0, 0), c( 1, 0, 1, 0, 0, 1, 0, 0, 0), c( 0,-1,-1, -1,-1,-1, 0, 0, 0), c( 1, 1, 0, 0, 1, 0, 0, 0, 0), c( 0, 2, 0, 0, 0, 0, 1, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0), c(-1, 0,-1, 0,-1, 0,-1,-1, 0), c( 1, 0, 1, 0, 0, 1, 0, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0), c( 0, 0, 2, 0, 0, 0, 0, 0, 1), c(-1,-1, 0, 0, 0,-1, 0,-1,-1), c( 0,-1,-1, -1,-1,-1, 0, 0, 0), c(-1, 0,-1, 0,-1, 0,-1,-1, 0), c(-1,-1, 0, 0, 0,-1, 0,-1,-1), c(-2,-2,-2, 1, 2, 2, 1, 2, 1)) #análogo à matriz da pág.71 e116.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e116.catdata,X=e116.X2) #simetr.form.log-lin. e116.loglinwls2<-funlinWLS(model=c("lin","log"),obj=e116.catdata,U=e116.A1) #simetr.form.log-lin. e116.A2<-rbind(cbind(kronecker(diag(3),t(rep(1,4))),matrix(0,3,4)), kronecker(t(rep(1,4)),cbind(diag(3),rep(0,3)))) e116.linwls3<-funlinWLS(model="lin",obj=e116.catdata,A1=e116.A2,X=kronecker(rep(1,2),diag(3)))#HM #u_1^A,u_2^A,u_3^A, u_1^B,u_2^B,u_3^B, u_{11},u_{12},u_{13},u_{22},u_{23},u_{33} e116.X3<-rbind(c( 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0), c( 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0), c( 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0), c( 1, 0, 0, -1,-1,-1, -1,-1,-1, 0, 0, 0), c( 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0), c( 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0), c( 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0), c( 0, 1, 0, -1,-1,-1, 0,-1, 0,-1,-1, 0), c( 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0), c( 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0), c( 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1), c( 0, 0, 1, -1,-1,-1, 0, 0,-1, 0,-1,-1), c(-1,-1,-1, 1, 0, 0, -1,-1,-1, 0, 0, 0), c(-1,-1,-1, 0, 1, 0, 0,-1, 0,-1,-1, 0), c(-1,-1,-1, 0, 0, 1, 0, 0,-1, 0,-1,-1), c(-1,-1,-1, -1,-1,-1, 1, 2, 2, 1, 2, 1)) #e92.X2<-cbind(e92.X3[,1]+e92.X3[,4],e92.X3[,2]+e92.X3[,5],e92.X3[,3]+e92.X3[,6],e92.X3[,7:12]) cbind(diag(15),rep(-1,15))%*%e116.X3 #matriz X_G do livro e116.loglinwls3<-funlinWLS(model=c("lin","log"),obj=e116.catdata,X=e116.X3) #quasi-simetria #u_1,u_2,u_3, u_{11},u_{12},u_{13},u_{22},u_{23},u_{33},gama e116.X4<-rbind(c( 2, 0, 0, 1, 0, 0, 0, 0, 0, 0), c( 1, 1, 0, 0, 1, 0, 0, 0, 0, 1), c( 1, 0, 1, 0, 0, 1, 0, 0, 0, 1), c( 0,-1,-1, -1,-1,-1, 0, 0, 0, 1), c( 1, 1, 0, 0, 1, 0, 0, 0, 0, 0), c( 0, 2, 0, 0, 0, 0, 1, 0, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0, 1), c(-1, 0,-1, 0,-1, 0,-1,-1, 0, 1), c( 1, 0, 1, 0, 0, 1, 0, 0, 0, 0), c( 0, 1, 1, 0, 0, 0, 0, 1, 0, 0), c( 0, 0, 2, 0, 0, 0, 0, 0, 1, 0), c(-1,-1, 0, 0, 0,-1, 0,-1,-1, 1), c( 0,-1,-1, -1,-1,-1, 0, 0, 0, 0), c(-1, 0,-1, 0,-1, 0,-1,-1, 0, 0), c(-1,-1, 0, 0, 0,-1, 0,-1,-1, 0), c(-2,-2,-2, 1, 2, 2, 1, 2, 1, 0)) e116.loglinwls4<-funlinWLS(model=c("lin","log"),obj=e116.catdata,X=e116.X4) #simetr.condicional e116.loglinwls5<-funlinWLS(model=c("lin","log"),obj=e116.catdata,A1=e116.A1,XL=rep(1,6))#idem #Exemplo 11.7 (p.390) / 9.3 (p.269) / 1.6 (p.11): Problema dos defeitos de fibras têxteis e117.TF<-rbind(c(28,40,68), c( 5,21,49), c( 1, 4,15)) e117.catdata<-readCatdata(TF=e117.TF) e117.A<-kronecker(diag(3),cbind(diag(2),rep(-1,2))) e117.X1<-rbind(c(1,0,0,0), c(0,1,0,0), c(1,0,2,0), c(0,1,1,0), c(1,0,0,2), c(0,1,0,1)) e117.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e117.catdata,A1=e117.A, XL=e117.X1) #efeito de linha round(exp(-e117.loglinwls1$beta[3:4]),2) round(exp(e117.loglinwls1$beta[3]-e117.loglinwls1$beta[4]),2) waldTest(obj=e117.loglinwls1,C=cbind(0*diag(2),diag(2))) e117.X2<-rbind(c(1,0,0), c(0,1,0), c(1,0,2), c(0,1,1), c(1,0,4), c(0,1,2)) e117.loglinwls2<-funlinWLS(model=c("lin","log"),obj=e117.catdata,A1=e117.A, XL=e117.X2) #associação uniforme round(exp(-e117.loglinwls2$beta[3]),2) round(exp(-e117.loglinwls2$beta[3])*sqrt(e117.loglinwls2$Vbeta[3,3]),2) e117.A2<-rbind(c(1,-1, 0,-1, 1, 0, 0, 0,0), c(0, 1,-1, 0,-1, 1, 0, 0,0), c(0, 0, 0, 1,-1, 0,-1, 1,0), c(0, 0, 0, 0, 1,-1, 0,-1,1)) e117.loglinwls3<-funlinWLS(model=c("lin","log"),obj=e117.catdata,A1=e117.A2, XL=rep(1,4)) #associação uniforme #Exemplo 11.8 (p.392) / 10.2 (p.349) / 6.3 (p.156) / 1.2 (p.4): Problema do risco de cárie dentária e118.TF<-c(11,5,0,14,34,7,2,13,11) e118.catdata<-readCatdata(TF=e118.TF) e118.B<-rbind(c(1,-1,0),c(0,1,-1)) e118.A<-kronecker(e118.B,e118.B) e118.loglinwls1a<-funlinWLS(model=c("lin","log"), obj=e118.catdata,A1=e118.A,XL=rep(1,4)) #subst.padrão de freq=0 por 1/[9*97]=0.001145475 e118.loglinwls1b<-funlinWLS(model=c("lin","log"), obj=readCatdata(TF=c(11,5,1/(9*97),14,34,7,2,13,11)),A1=e118.A,XL=rep(1,4)) #parecido e118.loglinwls2a<-funlinWLS(model=c("lin","log"), obj=e118.catdata,A1=e118.A,XL=rep(1,4),zeroN=1/2) #subst.freq=0 por 1/2 e118.loglinwls2b<-funlinWLS(model=c("lin","log"), obj=readCatdata(TF=c(11,5,1/2,14,34,7,2,13,11)),A1=e118.A,XL=rep(1,4)) #parecido e118.loglinwls3a<-funlinWLS(model=c("lin","log"), obj=e118.catdata,A1=e118.A,XL=rep(1,4),zeroN=1/16) #subst.freq=0 por 1/16 e118.loglinwls3b<-funlinWLS(model=c("lin","log"), obj=readCatdata(TF=c(11,5,1/16,14,34,7,2,13,11)),A1=e118.A,XL=rep(1,4)) #parecido e118.loglinwls4 <-funlinWLS(model=c("lin","log"),obj=readCatdata(TF=1/2+ c(11,5,0,14,34,7,2,13,11)),A1=e118.A,XL=rep(1,4)) #soma 1/2 a todas as freqüências #ao deixar que a rotina substitua os zeros amostrais pelo valor contido em zeroN, a #substituição será realizada apenas onde for necessário (p/obter a matriz de covariâncias das #proporções), enquanto que ao substituir diretamente em TF, todas as quantidades serão afetadas #Exemplo 11.9 (p.393) / 10.3 (p.349) / 6.4 (p.157) / 1.5 (p.5): Problema do uso do fio dental e119.TF<-rbind(c(19,5,4, 2),c( 5,8,0,17),c(11,6,7, 6),c( 2,5,1,22)) e119.catdata<-readCatdata(TF=e119.TF) e119.A<-kronecker(diag(4),t(c(1,-1,-1,1))) e119.XL<-cbind(rep(1,4),c(1,1,0,0),c(0,1,0,1)) e119.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e119.catdata,A1=e119.A,XL=e119.XL) e119.loglinwls2<-funlinWLS(model=c("lin","log"),obj=readCatdata(TF= rbind(c(19,5,4,2),c(5,8,1/2,17),c(11,6,7,6),c(2,5,1,22))),A1=e119.A,XL=e119.XL) #livro e119.loglinwls3<-funlinWLS(model=c("lin","log"),obj=readCatdata(TF= rbind(c(19,5,4,2),c(5,8,1/16,17),c(11,6,7,6),c(2,5,1,22))),A1=e119.A,XL=e119.XL) #livro e119.loglinwls4<-funlinWLS(model=c("lin","log"),obj=readCatdata(TF= 1/2+rbind(c(19,5,4,2),c(5,8,0,17),c(11,6,7,6),c(2,5,1,22))),A1=e119.A,XL=e119.XL) #livro e119.loglinwls5<-funlinWLS(model=c("lin","log"),obj=readCatdata(TF= rbind(c(19,5,4,2),1/2+c(5,8,0,17),c(11,6,7,6),c(2,5,1,22))),A1=e119.A, XL=e119.XL) #adic.1/2 apenas às celas da subpopulação que teve 0 e119.loglinwls6<-funlinWLS(model=c("lin","log"),obj=readCatdata(TF= 1/2+rbind(c(19,5,4,2),c(5,8,0,17),c(11,6,7,6),c(2,5,1,22))),A1=e119.A,XL=e119.XL[,-2]) #livro round(exp(e119.loglinwls6$beta),2) round(exp(e119.loglinwls6$beta-qnorm(0.975)*sqrt(diag(e119.loglinwls6$Vbeta))),2) round(exp(e119.loglinwls6$beta+qnorm(0.975)*sqrt(diag(e119.loglinwls6$Vbeta))),2) #Exemplo 11.10 (p.395) / 10.4 (p.351) / 6.5 (p.160): Problema da complicação pulmonar e1110.catdata<-readCatdata(TF=cbind(c(737,243,39),c(48,74,21))) e1110.A<-rbind(c(0,-1,0,1,0,0),c(0,-1,0,0,0,1)) e1110.loglinwls1<-funlinWLS(model=c("lin","log"),obj=e1110.catdata,A1=e1110.A,X=diag(2)) e1110.loglinwls2<-funlinWLS(model=c("lin","log"),obj=e1110.catdata,A1=e1110.A,X=c(1,1)) round(exp(e1110.loglinwls1$beta),2) round(exp(e1110.loglinwls1$beta-1*qnorm(0.975)*sqrt(diag(e1110.loglinwls1$Vbeta))),2) round(exp(e1110.loglinwls1$beta+1*qnorm(0.975)*sqrt(diag(e1110.loglinwls1$Vbeta))),2) #Exemplo 11.11 (p.395) / 10.5 (p.353) / 6.6 (p.163) / 1.3 (p.4): Problema do peso de recém-nascidos e1111.TF<-rbind(c( 2, 11, 31),c( 5, 24, 95), c( 3, 32, 91),c( 11, 57, 238), c( 15, 58,134),c( 25,105, 445), c(130,362,695),c(231,694,2485), c( 94,225,340),c(105,339,1053)) e1111.catdata<-readCatdata(TF=e1111.TF) e1111.A1<-kronecker(diag(10),rbind(c(1,0,0),c(0,1,1),c(0,1,0),c(0,0,1))) e1111.A2<-kronecker(diag(10),kronecker(diag(2),t(c(1,-1)))) e1111.X1<-kronecker(rbind(c(1,0,0,0,0,0), c(1,0,0,0,0,1), c(1,1,0,0,0,0), c(1,1,0,0,0,1), c(1,0,1,0,0,0), c(1,0,1,0,0,1), c(1,0,0,1,0,0), c(1,0,0,1,0,1), c(1,0,0,0,1,0), c(1,0,0,0,1,1)),diag(2)) e1111.funlinwls1<-funlinWLS(model=c("lin","log","lin"),obj=e1111.catdata, A1=e1111.A1,A2=e1111.A2,X=e1111.X1) e1111.X2<-kronecker(rbind(c(1,0,0), c(1,0,1), c(1,1,0), c(1,1,1), c(1,2,0), c(1,2,1), c(1,3,0), c(1,3,1), c(1,4,0), c(1,4,1)),diag(2)) e1111.funlinwls2<-funlinWLS(model=c("lin","log","lin"),obj=e1111.catdata, A1=e1111.A1,A2=e1111.A2,X=e1111.X2) waldTest(e1111.funlinwls2,c(0,0,1,-1,0,0)) waldTest(e1111.funlinwls2,c(0,0,0,0,1,-1)) round(exp(e1111.funlinwls2$beta),2) round(exp(e1111.funlinwls2$beta-qnorm(0.975)*sqrt(diag(e1111.funlinwls2$Vbeta))),2) round(exp(e1111.funlinwls2$beta+qnorm(0.975)*sqrt(diag(e1111.funlinwls2$Vbeta))),2) e1111.funlinwls3<-funlinWLS(model=c("lin","log","lin"),obj=e1111.catdata,A1=e1111.A1, A2=e1111.A2,X=cbind(e1111.X2[,1:2],e1111.X2[,3]+e1111.X2[,4],e1111.X2[,5]+e1111.X2[,6])) round(exp(e1111.funlinwls3$beta),2) round(exp(e1111.funlinwls3$beta-qnorm(0.975)*sqrt(diag(e1111.funlinwls3$Vbeta))),2) round(exp(e1111.funlinwls3$beta+qnorm(0.975)*sqrt(diag(e1111.funlinwls3$Vbeta))),2) e1111.A3<-kronecker(diag(10),rbind(c(1,0,0),c(0,1,1),c(1,1,0),c(0,0,1))) e1111.funlinwls4<-funlinWLS(model=c("lin","log","lin"),obj=e1111.catdata,A1=e1111.A3, A2=e1111.A2,X=cbind(e1111.X2[,1:2],e1111.X2[,3]+e1111.X2[,4],e1111.X2[,5]+e1111.X2[,6])) round(exp(e1111.funlinwls4$beta),2) round(exp(e1111.funlinwls4$beta-qnorm(0.975)*sqrt(diag(e1111.funlinwls4$Vbeta))),2) round(exp(e1111.funlinwls4$beta+qnorm(0.975)*sqrt(diag(e1111.funlinwls4$Vbeta))),2) #Exemplo 11.12 (p.399) / 6.8 (p.169) / 1.2 (p.4): Problema do risco de cárie dentária e1112.TF<-c(11,5,0,14,34,7,2,13,11) e1112.catdata<-readCatdata(TF=e1112.TF) e1112.A1<-rbind( c(rep(c(1,0,0,0),2),1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e1112.A2<-rbind( cbind(diag(2),matrix(0,2,6)), cbind(matrix(0,3,2),kronecker(t(rep(1,2)),diag(3))) ) e1112.A3<-cbind( c(1,0),c(1,1),-c(2,1)%*%t(rep(1,3)) ) e1112.A4<-t(c(1,-1)) e1112.kappa<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e1112.catdata,A1=e1112.A1,A2=e1112.A2,A3=e1112.A3,A4=e1112.A4,PI1=-1,X=1) round(pnorm((e1112.kappa$beta-0.35)/sqrt(e1112.kappa$Vbeta)),3) #equivalente a round((1-pchisq(waldTest(obj=e1112.kappa,C=1,C0=0.35)$QwH,1))/2,3) round(e1112.kappa$beta+c(-1,1)*qnorm(0.975)*sqrt(e1112.kappa$Vbeta),3) #kappa ponderado (Spitzer, Cohen, Fleiss e Endicott, 1967) W1<-c(1,0.75,0,0.75,1,0.75,0,0.75,1) #pesos quadráticos (Fleiss e Cohen, 1973) W2<-c(1,0.5,0,0.5,1,0.5,0,0.5,1) #pesos absolutos (Cicchetti e Allison, 1971) e1112.w1A1<-rbind( t(W1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e1112.w2A1<-rbind( t(W2), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e1112.wA2<-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))) ) ) e1112.w1A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W1)) ) e1112.w2A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W2)) ) e1112.kappaw1<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e1112.catdata,A1=e1112.w1A1,A2=e1112.wA2,A3=e1112.w1A3,A4=e1112.A4,PI1=-1,X=1) e1112.kappaw2<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e1112.catdata,A1=e1112.w2A1,A2=e1112.wA2,A3=e1112.w2A3,A4=e1112.A4,PI1=-1,X=1) #Exemplo 11.13 (p.399): Problema da poluição por petróleo e1113.TF<-rbind( c(2,0,4,0,6,0,3,0,1,0,1,0,15), c(0,0,1,0,0,0,2,0,2,0,5,0,22), c(0,0,2,0,1,0,3,0,0,0,5,0,21), c(0,0,2,0,1,0,2,0,0,0,6,0,21), c(1,0,4,0,2,0,6,0,1,0,0,0,18), c(0,0,1,0,0,0,0,0,0,0,1,0,30), c(0,0,1,0,0,0,3,0,1,0,0,1,26), c(0,0,0,0,2,0,0,0,0,1,1,0,28)) e1113.catdata<-readCatdata(TF=e1113.TF) e1113.A1<-diag(8)%x%rbind(c(rep(c(1,0),6),0),c(rep((0:5)+0.5,rep(2,6)),6)) e1113.A2<-diag(8)%x%t(c(1,-1)) e1113.X<-cbind(rep(1,8),c(0,1)%x%rep(1,4),c(1,1)%x%rbind(rep(0,3),diag(3)), rbind(matrix(0,5,3),diag(3))) e1113.X2<-cbind(e1113.X[,1:2],rep(c(1,0,0,0),2)) #Análises das log-taxas de mortalidades e1113.ltm1<-funlinWLS(model=c("lin","log","lin"),obj=e1113.catdata, A1=e1113.A1,A2=e1113.A2,X=e1113.X) #modelo saturado waldTest(e1113.ltm1,cbind(matrix(0,3,5),diag(3))) e1113.ltm2<-funlinWLS(model=c("lin","log","lin"),obj=e1113.catdata, A1=e1113.A1,A2=e1113.A2,X=e1113.X[,-(6:8)]) #modelo reduzido e1113.ltm3<-funlinWLS(model=c("lin","log","lin"),obj=e1113.catdata, A1=e1113.A1,A2=e1113.A2,X=e1113.X2) #modelo reduzido2 #Análises das taxas de mortalidades e1113.tm1<-funlinWLS(model=c("exp","lin","log","lin"),obj=e1113.catdata, A1=e1113.A1,A2=e1113.A2,X=e1113.X) #modelo saturado waldTest(e1113.tm1,cbind(matrix(0,3,5),diag(3))) e1113.tm2<-funlinWLS(model=c("exp","lin","log","lin"),obj=e1113.catdata, A1=e1113.A1,A2=e1113.A2,X=e1113.X[,-(6:8)]) #modelo reduzido e1113.tm3<-funlinWLS(model=c("exp","lin","log","lin"),obj=e1113.catdata, A1=e1113.A1,A2=e1113.A2,X=e1113.X2) #modelo reduzido2 round(e1113.tm3$FH[6],3) round(e1113.tm3$FH[6]+c(-1,1)*qnorm(0.975)*sqrt(e1113.tm3$VFH[6,6]),3) round(-e1113.tm3$beta[2],3) round(-e1113.tm3$beta[2]+c(-1,1)*qnorm(0.975)*sqrt(e1113.tm3$Vbeta[2,2]),3) round(e1113.tm3$beta[3],3) round(e1113.tm3$beta[3]+c(-1,1)*qnorm(0.975)*sqrt(e1113.tm3$Vbeta[3,3]),3) #Exemplo 12.1 (pp.419, 427, 442): Problema da infecção urinária e121.raw<-data.frame( inicial=c(1,2,1,2,2,2,2,1,3,2,2,1,2,2,2,2,2,1,1,2,1,3,3,2,2, 2,2,3,2,2,3,0,1,1,1,1,0,0,1,2,1,1,2,2,2,2,3,3,2,3), dias14 =c(0,0,0,0,1,2,2,1,0,1,1,1,0,0,1,1,1,0,0,0,1,1,0,1,0, 0,3,0,2,0,2,0,1,0,0,1,0,0,0,0,1,0,2,2,0,2,1,0,1,0), dias21 =c(0,0,0,0,1,2,2,1,0,2,3,0,0,0,1,1,0,0,0,0,0,0,0,1,0, 0,3,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,2,2,1,0,0,0,1,0)) table(e121.raw[,1]) table(e121.raw[,2]) table(e121.raw[,3]) e121.TF<-c(table(e121.raw[,3:1])) e121.catdata<-readCatdata(TF=e121.TF) e121.v1<-c(0,0,1,1) e121.A1<-rbind(e121.v1%x%rep(1,16) , rep(1,4)%x%e121.v1%x%rep(1,4) , rep(1,16)%x%e121.v1) e121.v2<-c(0,1,2,3) e121.A2<-rbind(e121.v2%x%rep(1,16) , rep(1,4)%x%e121.v2%x%rep(1,4) , rep(1,16)%x%e121.v2) e121.X<-rbind(c(1,0,0),c(1,1,0),c(1,0,1)) e121.propwls<-funlinWLS(model="lin",obj=e121.catdata,A1=e121.A1,X=e121.X) waldTest(e121.propwls,rbind(c(0,1,0),c(0,0,1))) waldTest(e121.propwls,c(0,1,-1)) e121.scorwls<-funlinWLS(model="lin",obj=e121.catdata,A1=e121.A2,X=e121.X) waldTest(e121.scorwls,rbind(c(0,1,0),c(0,0,1))) waldTest(e121.scorwls,c(0,1,-1)) e121.X2<-rbind(c(1,0),c(1,1),c(1,1)) e121.propwls2<-funlinWLS(model="lin",obj=e121.catdata,A1=e121.A1,X=e121.X2) e121.scorwls2<-funlinWLS(model="lin",obj=e121.catdata,A1=e121.A2,X=e121.X2) round(e121.propwls2$beta[1],2) round(e121.propwls2$beta[1]+c(-1,1)*qnorm(0.975)*sqrt(e121.propwls2$Vbeta[1,1]),2) round(-e121.propwls2$beta[2],2) round(-e121.propwls2$beta[2]+c(-1,1)*qnorm(0.975)*sqrt(e121.propwls2$Vbeta[2,2]),2) round(e121.scorwls2$beta[1],2) round(e121.scorwls2$beta[1]+c(-1,1)*qnorm(0.975)*sqrt(e121.scorwls2$Vbeta[1,1]),2) round(-e121.scorwls2$beta[2],2) round(-e121.scorwls2$beta[2]+c(-1,1)*qnorm(0.975)*sqrt(e121.scorwls2$Vbeta[2,2]),2) e121.v3<-rbind(c(1,0,0,0),c(0,1,1,1),c(1,1,0,0),c(0,0,1,1),c(1,1,1,0),c(0,0,0,1)) e121.A31<-rbind(e121.v3%x%t(rep(1,16)) , t(rep(1,4))%x%e121.v3%x%t(rep(1,4)) , t(rep(1,16))%x%e121.v3) e121.A32<-diag(9)%x%t(c(1,-1)) e121.X3<-cbind(rep(1,3)%x%diag(3),cbind(c(0,1,0),c(0,0,1))%x%rep(1,3)) e121.chprwls<-funlinWLS(model=c("lin","log","lin"),obj=e121.catdata, A1=e121.A31,A2=e121.A32,X=e121.X3) round(e121.chprwls$beta,2) round(exp(e121.chprwls$beta),2) round(exp(e121.chprwls$beta-qnorm(0.975)*sqrt(diag(e121.chprwls$Vbeta))),2) round(exp(e121.chprwls$beta+qnorm(0.975)*sqrt(diag(e121.chprwls$Vbeta))),2) waldTest(e121.chprwls,c(0,0,0,1,-1)) e121.X4<-cbind(e121.X3[,1:3],e121.X3[,4]+e121.X3[,5]) e121.chprwls2<-funlinWLS(model=c("lin","log","lin"),obj=e121.catdata, A1=e121.A31,A2=e121.A32,X=e121.X4) round(e121.chprwls2$beta[4],2) round(exp(e121.chprwls2$beta[4]),2) round(exp(e121.chprwls2$beta[4]+c(-1,1)*qnorm(0.975)*sqrt(e121.chprwls$Vbeta[4,4])),2) table(e121.raw[,2:1]) table(e121.raw[,c(3,1)]) table(e121.raw[,3:2]) e121.A41<-rbind(diag(16)%x%t(rep(1,4)),diag(4)%x%t(rep(1,4))%x%diag(4),t(rep(1,4))%x%diag(16)) e121.A42<-rbind(c(1,0,0),c(0,0,1))%x%rbind(c(0,0,0,0,1,0,0,0,1,1,0,0,1,1,1,0)) e121.probsbiv<-funlinWLS(model=c("lin","lin"),obj=e121.catdata,A1=e121.A41,A2=e121.A42, X=diag(2)) waldTest(e121.probsbiv,c(1,-1)) round(e121.probsbiv$beta,2) round(e121.probsbiv$beta-qnorm(0.975)*sqrt(diag(e121.probsbiv$Vbeta)),2) round(e121.probsbiv$beta+qnorm(0.975)*sqrt(diag(e121.probsbiv$Vbeta)),2) #pp.442-444 e121.raw2<-data.frame(paciente=rep(1:50,3),corrimento=with(e121.raw,c(inicial,dias14,dias21)), avaliacao=c(rep("inicial",50),rep("dias14",50),rep("dias21",50))) e121.raw2<-e121.raw2[order(e121.raw2$paciente),] require(gee) e121.gee<-with(e121.raw2,gee(I(ifelse(corrimento>=2,1,0))~C(avaliacao,base=3),id=paciente, family=quasi(link="identity",variance="mu(1-mu)"),corstr="exchangeable",scale.value=1, scale.fix=TRUE)) e121.gee2<-with(e121.raw2,gee(I(ifelse(corrimento>=2,1,0))~C(avaliacao,base=3),id=paciente, family=quasi(link="identity",variance="mu(1-mu)"),corstr="unstructured",scale.value=1, scale.fix=TRUE)) round(matrix(c(e121.propwls$beta,sqrt(diag(e121.propwls$Vbeta)), e121.gee$coef,sqrt(diag(e121.gee$robust)),sqrt(diag(e121.gee$naive)), e121.gee2$coef,sqrt(diag(e121.gee2$robust)),sqrt(diag(e121.gee2$naive))),3),3) #Tabela 12.10 cov2cor(e121.propwls$VFH) #compare a estimativa da matriz de correlações da análise por MQG... e121.gee2$work #...com a estimativa da matriz de de correlação de trabalho não estruturada e121.gee$work #Exemplo 12.2 (pp.421, 431): Problema da sensibilidade dentinária e122.TF<-rbind(c(22,1,3,6),c(12,10,7,4),c(10,6,12,3),c(5,13,11,3)) e122.catdata<-readCatdata(TF=e122.TF) e122.A1<-diag(4)%x%(1-rbind(diag(2)%x%t(c(1,1)),t(c(1,1))%x%diag(2))) e122.A2<-diag(4)%x%t(c(-1,1,1,-1)) e122.X<-cbind(rep(1,4),c(1,1,-1,-1),rep(c(1,-1),2),c(1,-1,-1,1)) e122.wls<-funlinWLS(model=c("lin","log","lin"),obj=e122.catdata,A1=e122.A1,A2=e122.A2, X=e122.X) e122.wls2<-funlinWLS(model=c("lin","log","lin"),obj=e122.catdata,A1=e122.A1,A2=e122.A2, X=rep(1,4)) round(exp(e122.wls2$beta),2) round(exp(e122.wls2$beta+c(-1,1)*qnorm(0.975)*sqrt(e122.wls2$Vbeta)),2) #pp.444-445 e122.A3<-diag(4)%x%diag(2)%x%t(c(-1,1)) e122.X2<-e122.X %x% diag(2) e122.wls3<-funlinWLS(model=c("lin","log","lin"),obj=e122.catdata,A1=e122.A1,A2=e122.A3, X=e122.X2) waldTest(e122.wls3,c(0,0,1,1,0,0,0,0)) waldTest(e122.wls3,c(0,0,0,0,1,1,0,0)) waldTest(e122.wls3,c(0,0,0,0,0,0,1,1)) waldTest(e122.wls3,rep(1,4)%x%c(1,-1)) waldTest(e122.wls3,c(0,0,1,-1,0,0,0,0)) waldTest(e122.wls3,c(0,0,0,0,1,-1,0,0)) waldTest(e122.wls3,c(0,0,0,0,0,0,1,-1)) cov2cor(e122.wls3$VFu) #Exemplo 12.3 (pp.422, 432): Problema da maturação do colo do útero e123.raw<-data.frame( paridade=c(rep("N",46),rep("M",37)), cons00=c(rep(2,28),1,1,2,2,1,2,2,2,2,1,1,1,1,1,2,1,1,1, 2,2,2,1,2,2,2,2,rep(1,11),2,2,1,2,1,1,1,2,2,1,1,2,1,2,2,2,2,2), cons24=c(2,0,2,0,1,0,0,0,2,0,0,0,0,2,2,2,0,2,2,1,1,1,0,2,1,2,0,2,1,1,2,0,0,1,0,0,0,0,0,1, 0,1,1,1,1,0, 2,2,2,1,1,2,2,2,1,1,1,0,1,0,1,1,0,0,0,1,0,0,2,1,1,0,0,2,0,1,1,1,0,2,0,2,0), cons48=c(1,0,2,0,1,rep(0,8),1,2,1,0,0,0,1,1,1,0,1,0,2,0,0,1,1,2,0,0,1,rep(0,9),1,1,0, 2,2,2,1,1,2,1,1,1,0,1,0,1,0,1,1,rep(0,7),1,1,0,0,2,0,1,1,0,0,2,0,2,0), cons72=c(rep(0,14),1,rep(0,10),2,0,0,0,1,2,rep(0,12),1,1,0, 2,1,2,0,0,0,1,0,0,0,1,0,1,0,0,1,rep(0,7),1,rep(0,9),2,0,2,0), cons96=c(rep(0,14),1,rep(0,10),2,0,0,0,1,2,rep(0,12),1,1,0, 2,0,2,rep(0,7),1,0,1,0,0,1,rep(0,7),1,rep(0,9),2,0,2,0)) e123.agr1<-with(e123.raw,aggregate(rep(1,nrow(e123.raw)),by=list( cons00=cons00,cons24=cons24,cons48=cons48,cons72=cons72,cons96=cons96,paridade=paridade),sum)) e123.agr2<-with(e123.raw,aggregate(rep(1,nrow(e123.raw)),by=list( cons00=cons00,cons24=cons24,cons48=cons48,cons72=cons72,paridade=paridade),sum)) i1<-0;e123.perfil1<-numeric(20);e123.TF1<-matrix(0,2,20) i2<-0;e123.perfil2<-numeric(14);e123.TF2<-matrix(0,2,14) for(i00 in 2:1) { for(i24 in i00:0) { for(i48 in i24:0) { for(i72 in i48:0) { for(i96 in i72:0) { i1<-i1+1 e123.perfil1[i1]<-10000*i00+1000*i24+100*i48+10*i72+i96 acess<-with(e123.agr1,x[cons00==i00 & cons24==i24 & cons48==i48 & cons72==i72 & cons96==i96 & paridade=="M"]) e123.TF1[1,i1]<-ifelse(length(acess)==0,0,acess) acess<-with(e123.agr1,x[cons00==i00 & cons24==i24 & cons48==i48 & cons72==i72 & cons96==i96 & paridade=="N"]) e123.TF1[2,i1]<-ifelse(length(acess)==0,0,acess) } i2<-i2+1 e123.perfil2[i2]<-1000*i00+100*i24+10*i48+i72 acess<-with(e123.agr1,x[cons00==i00 & cons24==i24 & cons48==i48 & cons72==i72 & paridade=="M"]) e123.TF2[1,i2]<-ifelse(length(acess)==0,0,acess) acess<-with(e123.agr1,x[cons00==i00 & cons24==i24 & cons48==i48 & cons72==i72 & paridade=="N"]) e123.TF2[2,i2]<-ifelse(length(acess)==0,0,acess) } } } } e123.catdata<-readCatdata(TF=e123.TF2) e123.perfil2 # 2222 2221 2220 2211 2210 2200 2111 2110 2100 2000 1111 1110 1100 1000 e123.A0<-diag(2)%x%rbind( #psi_{i12} c(0,0,0,0,0,0,0,0,0,0,0,0,0,1), c(0,0,0,0,0,0,0,0,0,0,1,1,1,0), c(0,0,0,0,0,0,0,0,0,1,0,0,0,0), c(0,0,0,0,0,0,1,1,1,0,0,0,0,0), c(1,1,1,1,1,1,0,0,0,0,0,0,0,0), #psi_{i23} c(0,0,0,0,0,0,0,0,0,1,0,0,0,1), c(0,0,0,0,0,0,0,0,1,0,0,0,1,0), c(0,0,0,0,0,0,1,1,0,0,1,1,0,0), c(0,0,0,0,0,1,0,0,0,0,0,0,0,0), c(0,0,0,1,1,0,0,0,0,0,0,0,0,0), c(1,1,1,0,0,0,0,0,0,0,0,0,0,0), #psi_{i34} c(0,0,0,0,0,1,0,0,1,1,0,0,1,1), c(0,0,0,0,1,0,0,1,0,0,0,1,0,0), c(0,0,0,1,0,0,1,0,0,0,1,0,0,0), c(0,0,1,0,0,0,0,0,0,0,0,0,0,0), c(0,1,0,0,0,0,0,0,0,0,0,0,0,0), c(1,0,0,0,0,0,0,0,0,0,0,0,0,0)) e123.A1<-diag(2)%x%(diag(3)%x%rbind(c(0,1,0,1,0,0),c(0,1,1,1,1,1)))[,-1] e123.A2<-diag(6)%x%t(c(1,-1)) e123.probtrans<-funlinWLS(model=c("exp","lin","log","lin","lin"),obj=e123.catdata,A1=e123.A0, A2=e123.A1,A3=e123.A2,X=diag(6)) waldTest(e123.probtrans,cbind(diag(3),-diag(3))) e123.probtrans2<-funlinWLS(model=c("exp","lin","log","lin","lin"),obj=e123.catdata,A1=e123.A0, A2=e123.A1,A3=e123.A2,X=rep(1,2)%x%diag(3)) #Exemplo 13.2 (p.466) e132.TF<-c(7,11,2,3,9,5,0,10,4, 8,7,3,0, 0,7,14,7) e132.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))) ) ) e132.Rp<-c(4,4) e132.catdata<-readCatdata(TF=e132.TF,Zp=e132.Zp,Rp=e132.Rp) #p_{ij} e \hat{\sigma}(p_{ij}) e132.satmarml<-satMarML(e132.catdata) #\hat{\theta}_{ij}, \hat{\sigma}(...), Q_V e Q_P(M_2|M1) e132.satmarml$alphast #Tabela 13.3 - EMV das prob.condicionais de omissão e132.TF2<-c(7,11,2,3,9,5,1e-5,10,4, 8,7,3,0, 0,7,14,7) #subst.zero por valor peq. e132.catdata2<-readCatdata(TF=e132.TF2,Zp=e132.Zp,Rp=e132.Rp) e132.satmarml2<-satMarML(e132.catdata2) e132.U<-rbind(c(0, 1,1,-1,0,0,-1, 0), c(0,-1,0, 1,0,1, 0,-1) ) e132.linml<-linML(e132.satmarml2,U=e132.U) #\hat{\theta}_{ij}(H), \hat{\sigma}, Q_V(H|M1) e132.linwls<-funlinWLS(model="lin",obj=e132.satmarml2,U=e132.U) #abordagem híbrida e132.kA1<-rbind( c(rep(c(1,0,0,0),2),1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e132.kA2<-rbind( cbind(diag(2),matrix(0,2,6)), cbind(matrix(0,3,2),kronecker(t(rep(1,2)),diag(3))) ) e132.kA3<-cbind( c(1,0),c(1,1),-c(2,1)%*%t(rep(1,3)) ) e132.kA4<-t(c(1,-1)) e132.kappa<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e132.satmarml,A1=e132.kA1,A2=e132.kA2,A3=e132.kA3,A4=e132.kA4,PI1=-1,X=1) #Estes resultados não estão no livro, mas ilustram o cálculo do kappa ponderado W1<-c(1,0.75,0,0.75,1,0.75,0,0.75,1) #pesos quadráticos Fleiss e Cohen (1973) W2<-c(1,0.5,0,0.5,1,0.5,0,0.5,1) #pesos absolutos Agresti (2002) e132.kw1A1<-rbind( t(W1), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e132.kw2A1<-rbind( t(W2), rep(1,9), kronecker(diag(3),t(rep(1,3))), kronecker(t(rep(1,3)),diag(3)) ) e132.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))))) e132.kw1A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W1)) ) e132.kw2A3<-cbind( c(1,0),c(1,1),kronecker(-c(2,1),t(W2)) ) e132.kappaw1<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e132.satmarml,A1=e132.kw1A1,A2=e132.kwA2,A3=e132.kw1A3,A4=e132.kA4,PI1=-1,X=1) e132.kappaw2<-funlinWLS(model=c("add","exp","lin","log","lin","exp","lin","log","lin"), obj=e132.satmarml,A1=e132.kw2A1,A2=e132.kwA2,A3=e132.kw2A3,A4=e132.kA4,PI1=-1,X=1) #Exemplo 13.3 (p.469) e133.TF<-c(77,87,94,70,67,36,143,78, 14,8,3,9, 25,18,43,16, 14,12) e133.Zp<-cbind(kronecker(diag(4),rep(1,2)), kronecker(diag(2),kronecker(rep(1,2),diag(2))), kronecker(diag(2),rep(1,4)) ) e133.Rp<-c(4,4,2) e133.catdata<-readCatdata(TF=e133.TF,Zp=e133.Zp,Rp=e133.Rp) #p_{ij}, \hat{\sigma}(p_{ij}) e133.satmcarml<-satMarML(e133.catdata,missing="MCAR")#\hat{\theta}_{ij}, Q_V e Q_P(M_2|M1) e133.satmcarwls<-satMcarWLS(e133.catdata) #\tilde{\theta}_{ij},\hat{\sigma}(...), QN_(H_0) e133.X<-rbind(c( 1, 1, 1, 1, 1, 1), c( 1, 1,-1, 1,-1,-1), c( 1,-1, 1, -1, 1,-1), c( 1,-1,-1, -1,-1, 1), c(-1, 1, 1, -1,-1, 1), c(-1, 1,-1, -1, 1,-1), c(-1,-1, 1, 1,-1,-1), c(-1,-1,-1, 1, 1, 1)) e133.loglinml<-loglinML(obj=e133.satmcarml,X=e133.X) e133.loglinwls<-funlinWLS(model=c("lin","log"),obj=e133.satmcarwls,X=e133.X) e133.loglinhib<-funlinWLS(model=c("lin","log"),obj=e133.satmcarml,X=e133.X) e133.loglinml2<-loglinML(obj=e133.satmcarml,X=e133.X[,-4]) e133.loglinml3<-loglinML(obj=e133.satmcarml,X=e133.X[,-5]) e133.loglinml4<-loglinML(obj=e133.satmcarml,X=e133.X[,-6]) e133.loglinml2$QvH-e133.loglinml$QvH e133.loglinml3$QvH-e133.loglinml$QvH e133.loglinml4$QvH-e133.loglinml$QvH 1-pchisq(e133.loglinml2$QvH-e133.loglinml$QvH,1) 1-pchisq(e133.loglinml3$QvH-e133.loglinml$QvH,1) 1-pchisq(e133.loglinml4$QvH-e133.loglinml$QvH,1) #note delta^{A(1)}=delta^{A(2)} sob ausência de interação de 2a.ordem c(1,-1,-1,1,0,0,0,0)%*%log(e133.loglinml$thetaH) c(0,0,0,0,1,-1,-1,1)%*%log(e133.loglinml$thetaH) c(1,-1,-1,1,0,0,0,0)%*%e133.X c(0,0,0,0,1,-1,-1,1)%*%e133.X round(4*e133.loglinml$beta[6:4],3) round(4*sqrt(diag(e133.loglinml$Vbeta))[6:4],3) round(4*e133.loglinwls$beta[6:4],3) round(4*sqrt(diag(e133.loglinwls$Vbeta))[6:4],3) round(4*e133.loglinhib$beta[6:4],3) round(4*sqrt(diag(e133.loglinhib$Vbeta))[6:4],3) #Exemplo 13.4 (p.472) / 13.1 (p.454) e134.TF<-c(12,4,5,2, 50,31, 27,12) e134.Zp<-cbind(kronecker(diag(2),rep(1,2)),kronecker(rep(1,2),diag(2))) e134.Rp<-c(2,2) e134.catdata<-readCatdata(TF=e134.TF,Zp=e134.Zp,Rp=e134.Rp) e134.satmcarml<-satMarML(e134.catdata,miss="MCAR")#MV MCAR Tabs.13.10 (Inf.Fisher),13.11,13.12 e134.satmarml<-satMarML(e134.catdata)#Results.MV MAR_{sat} Tabs.13.11-13.12 e134.satmcarwls<-satMcarWLS(e134.catdata) e134.A<-rbind(c(1,1,0,0),c(1,0,1,0)) e134.hmmcarml<-linML(e134.satmcarml,A=e134.A,X=rep(1,2)) #Tab.13.12/13.13 e134.hmmcarwls<-funlinWLS(model="lin",obj=e134.satmcarwls,A1=e134.A,X=rep(1,2)) #Tab.13.12 e134.hmmcarhib<-funlinWLS(model="lin",obj=e134.satmcarml,A1=e134.A,X=rep(1,2)) e134.TF2<-c(e134.TF,24) #para mecanismos MNAR, cenários de omissão total trazem inf.na estim. mnarsat.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1,\alpha_2 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7];a2<-p[8] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))) )+ n121*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (exp(a20+a1+a2)/(1+exp(a20+a1+a2))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))) )+ n22*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (1/(1+exp(a20+a1+a2))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))* (exp(a30+a1+a2)/(1+exp(a30+a1+a2))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))) + t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) ) value } mnarsat.der<-deriv3(~-( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))) )+ n121*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (exp(a20+a1+a2)/(1+exp(a20+a1+a2))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))) )+ n22*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (1/(1+exp(a20+a1+a2))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))* (exp(a30+a1+a2)/(1+exp(a30+a1+a2))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))) + t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) ),c("t11","t12","t21","a10","a20","a30","a1","a2"), c("t11","t12","t21","a10","a20","a30","a1","a2", "n111","n112","n121","n122","n21","n22","n31","n32","N4") ) #obtém o gradiente e a hessiana analiticamente mnarsat.esp<-function(p,N){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1,\alpha_2 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7];a2<-p[8] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))), t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))*(exp(a20+a1+a2)/(1+exp(a20+a1+a2))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))), t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))*(1/(1+exp(a20+a1+a2))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))), t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(exp(a30+a1+a2)/(1+exp(a30+a1+a2))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))), t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))), (1-t11-t12-t21)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) value } mnarred.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10+a1)/(1+exp(a10+a1)))* (exp(a20+a1)/(1+exp(a20+a1))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10+a1)/(1+exp(a10+a1)))* (1/(1+exp(a20+a1))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1)))* (exp(a30+a1)/(1+exp(a30+a1))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) ) ) value } mnarred.der<-deriv3(~-( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10+a1)/(1+exp(a10+a1)))* (exp(a20+a1)/(1+exp(a20+a1))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10+a1)/(1+exp(a10+a1)))* (1/(1+exp(a20+a1))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1)))* (exp(a30+a1)/(1+exp(a30+a1))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t21)*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) ) ),c("t11","t12","t21","a10","a20","a30","a1"), c("t11","t12","t21","a10","a20","a30","a1", "n111","n112","n121","n122","n21","n22","n31","n32","N4") ) mnarred.esp<-function(p,N){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t21*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t21*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t21*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), (1-t11-t12-t21)*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), t21*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))), (1-t11-t12-t21)*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) ) value } marred.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(1/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) ) value } marred.der<-deriv3(~-( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(1/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) ),c("t11","t12","t21","a10","a20","a30","a1"), c("t11","t12","t21","a10","a20","a30","a1", "n111","n112","n121","n122","n21","n22","n31","n32","N4") ) marred.esp<-function(p,N){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6];a1<-p[7] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t21*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))), (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t21*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), t21*(1/(1+exp(a10)))*(1/(1+exp(a30))), (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) value } mcar.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30} t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n122*log( (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(1/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) ) value } mcar.der<-deriv3(~-( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n122*log( (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t21*(1/(1+exp(a10)))*(1/(1+exp(a30))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) ),c("t11","t12","t21","a10","a20","a30"), c("t11","t12","t21","a10","a20","a30", "n111","n112","n121","n122","n21","n22","n31","n32","N4") ) mcar.esp<-function(p,N){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30} t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5];a30<-p[6] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t21*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t21*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), t21*(1/(1+exp(a10)))*(1/(1+exp(a30))), (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) value } mcarred.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10}=\alpha_{30},\alpha_{20} t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n122*log( (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) + t21*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a10))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a10))) + t21*(1/(1+exp(a10)))*(1/(1+exp(a10))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a10))) ) ) value } mcarred.der<-deriv3(~-( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t21*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n122*log( (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) + t21*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a10))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a10))) + t21*(1/(1+exp(a10)))*(1/(1+exp(a10))) + (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a10))) ) ),c("t11","t12","t21","a10","a20"), c("t11","t12","t21","a10","a20", "n111","n112","n121","n122","n21","n22","n31","n32","N4") ) mcarred.esp<-function(p,N){ #p=\theta_{11},\theta_{12},\theta_{21},\alpha_{10}=\alpha_{30},\alpha_{20} t11<-p[1];t12<-p[2];t21<-p[3] a10<-p[4];a20<-p[5] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t21*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t21*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), (1-t11-t12-t21)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t11*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), t21*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), (1-t11-t12-t21)*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), t11*(1/(1+exp(a10)))*(1/(1+exp(a10))), t12*(1/(1+exp(a10)))*(1/(1+exp(a10))), t21*(1/(1+exp(a10)))*(1/(1+exp(a10))), (1-t11-t12-t21)*(1/(1+exp(a10)))*(1/(1+exp(a10))) ) value } require(geoR) #.nlmP adapta nlm p/restringir o espaço paramétrico. Isso é importante, inipars<-c(0.25,0.25,0.25,0,0,0,0,0) #pois mecanismos MNAR resultam facilmente em estims. minpars<-c(0,0,0,-Inf,-Inf,-Inf,-Inf,-Inf) #p/probs. >1 ou <0 quando não se usa o EM ou maxpars<-c(1,1,1,Inf,Inf,Inf,Inf,Inf) #funções ligações próprias para probs.(e.g.,logito) mnarsat<-.nlmP(objfunc=mnarsat.mlv,params=inipars , lower=minpars ,upper=maxpars ,hessian=T) mnarred<-.nlmP(objfunc=mnarred.mlv,params=inipars[-8] , lower=minpars[-8] ,upper=maxpars[-8] ,hessian=T) marred <-.nlmP(objfunc=marred.mlv ,params=inipars[-8] , lower=minpars[-8] ,upper=maxpars[-8] ,hessian=T) mcar <-.nlmP(objfunc=mcar.mlv ,params=inipars[-(7:8)], lower=minpars[-(7:8)],upper=maxpars[-(7:8)],hessian=T) mcarred<-.nlmP(objfunc=mcarred.mlv,params=inipars[-(6:8)], lower=minpars[-(6:8)],upper=maxpars[-(6:8)],hessian=T) p<-mnarsat$est mnarsat.infobs<-attr(mnarsat.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],e134.TF2[1],e134.TF2[2], e134.TF2[3],e134.TF2[4],e134.TF2[5],e134.TF2[6],e134.TF2[7],e134.TF2[8],e134.TF2[9]), "hessian")[1,,] #informação observada obtida analiticamente p<-mnarred$est mnarred.infobs<-attr(mnarred.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],e134.TF2[1],e134.TF2[2], e134.TF2[3],e134.TF2[4],e134.TF2[5],e134.TF2[6],e134.TF2[7],e134.TF2[8],e134.TF2[9]), "hessian")[1,,] p<-marred$est marred.infobs<-attr(marred.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],e134.TF2[1],e134.TF2[2], e134.TF2[3],e134.TF2[4],e134.TF2[5],e134.TF2[6],e134.TF2[7],e134.TF2[8],e134.TF2[9]), "hessian")[1,,] p<-mcar$est mcar.infobs<-attr(mcar.der(p[1],p[2],p[3],p[4],p[5],p[6],e134.TF2[1],e134.TF2[2], e134.TF2[3],e134.TF2[4],e134.TF2[5],e134.TF2[6],e134.TF2[7],e134.TF2[8],e134.TF2[9]), "hessian")[1,,] p<-mcarred$est mcarred.infobs<-attr(mcarred.der(p[1],p[2],p[3],p[4],p[5],e134.TF2[1],e134.TF2[2], e134.TF2[3],e134.TF2[4],e134.TF2[5],e134.TF2[6],e134.TF2[7],e134.TF2[8],e134.TF2[9]), "hessian")[1,,] mnarsat.infobs2<-mnarsat$hess #informação observada obtida numericamente mnarred.infobs2<-mnarred$hess marred.infobs2<-marred$hess mcar.infobs2<-mcar$hess mcarred.infobs2<-mcarred$hess round(mnarsat.infobs-mnarsat.infobs2,3) #note as diferenças round(mnarred.infobs-mnarred.infobs2,3) round(marred.infobs-marred.infobs2,3) round(mcarred.infobs-mcarred.infobs2,3) round(mcar.infobs-mcar.infobs2,3) mnarsat.cov<-solve(mnarsat.infobs) mnarred.cov<-solve(mnarred.infobs) marred.cov<-solve(marred.infobs) mcar.cov<-solve(mcar.infobs) mcarred.cov<-solve(mcarred.infobs) mnarsat.cov2<-solve(mnarsat.infobs2) mnarred.cov2<-solve(mnarred.infobs2) marred.cov2<-solve(marred.infobs2) mcar.cov2<-solve(mcar.infobs2) mcarred.cov2<-solve(mcarred.infobs2) round(mnarsat.cov-mnarsat.cov2,6) #as diferenças são maiores justamente round(mnarred.cov-mnarred.cov2,6) #para os parâmetros de interesse round(marred.cov-marred.cov2,6) round(mcarred.cov-mcarred.cov2,6) round(mcar.cov-mcar.cov2,6) b<-c(rep(0,3),1) B<-rbind(diag(3),rep(-1,3)) mnarsat.estp<-c(b+B%*%mnarsat$est[1:3]) mnarred.estp<-c(b+B%*%mnarred$est[1:3]) marred.estp<-c(b+B%*%marred$est[1:3]) mcar.estp<-c(b+B%*%mcar$est[1:3]) mcarred.estp<-c(b+B%*%mcarred$est[1:3]) mnarsat.covp<-B%*%mnarsat.cov[1:3,1:3]%*%t(B) mnarred.covp<-B%*%mnarred.cov[1:3,1:3]%*%t(B) marred.covp<-B%*%marred.cov[1:3,1:3]%*%t(B) mcar.covp<-B%*%mcar.cov[1:3,1:3]%*%t(B) mcarred.covp<-B%*%mcarred.cov[1:3,1:3]%*%t(B) #Recorde, no Exercício~13.5, item (c), que a estimativa da matriz de informação de Fisher #relativa a \theta sob o mecanismo MAR é igual à estimativa da matriz de inf.observada #relativa a \theta sob os mecanismos MAR e MCAR #Como a rotina satMarML utiliza a informação de Fisher, devemos comparar a estimativa #da matriz de covariâncias obtida sob o mecanismo MAR com as dos mecanismos MAR_red, #MCAR, MCAR_red baseadas na matriz de informação observada round(e134.satmarml$Vtheta-marred.covp,6) round(e134.satmarml$Vtheta-mcar.covp,6) round(e134.satmarml$Vtheta-mcarred.covp,6) round(e134.satmarml$Vtheta-B%*%marred.cov2[1:3,1:3]%*%t(B),6) round(e134.satmarml$Vtheta-B%*%mcar.cov2[1:3,1:3]%*%t(B),6) round(e134.satmarml$Vtheta-B%*%mcarred.cov2[1:3,1:3]%*%t(B),6) #Portanto, sugere-se calcular sempre a inf.obs.analítica! e134.satmarml$theta-marred.estp e134.satmarml$theta-mcar.estp e134.satmarml$theta-mcarred.estp #Recorde tb que as estimativas de MV dos \theta são sempre iguais para todos os #mecanismo MAR, MCAR e estruturas mais reduzidas destes. cbind(mnarsat.estp,mnarred.estp,marred.estp,mcar.estp,mcarred.estp) #Tab.13.10 mnarsat$est[-(1:3)];mnarred$est[-(1:3)] #EMV dos \alphas da Tabela 13.10 marred$est[-(1:3)];mcar$est[-(1:3)];mcarred$est[-(1:3)] cbind(sqrt(diag(mnarsat.covp)),sqrt(diag(mnarred.covp)),sqrt(diag(marred.covp)), sqrt(diag(mcar.covp)),sqrt(diag(mcarred.covp))) #erros padrões dos \thetas Tab.13.10 sqrt(diag(mnarsat.cov))[4:8];sqrt(diag(mnarred.cov))[4:7];sqrt(diag(marred.cov))[4:7] sqrt(diag(mcar.cov))[4:6];sqrt(diag(mcarred.cov))[4:5] #e.p.s dos \alphas Tab.13.10 mnarsat.wls<-funlinWLS(model="lin",theta=mnarsat.estp,Vtheta=mnarsat.covp,A1=t(c(0,1,-1,0)),X=1) mnarred.wls<-funlinWLS(model="lin",theta=mnarred.estp,Vtheta=mnarred.covp,A1=t(c(0,1,-1,0)),X=1) marred.wls<-funlinWLS(model="lin",theta=marred.estp,Vtheta=marred.covp,A1=t(c(0,1,-1,0)),X=1) mcar.wls<-funlinWLS(model="lin",theta=mcar.estp,Vtheta=mcar.covp,A1=t(c(0,1,-1,0)),X=1) mcarred.wls<-funlinWLS(model="lin",theta=mcarred.estp,Vtheta=mcarred.covp,A1=t(c(0,1,-1,0)),X=1) rbind(c(mnarsat.wls$beta,mnarred.wls$beta,marred.wls$beta,mcar.wls$beta,mcarred.wls$beta), c(mnarsat.wls$Vbeta,mnarred.wls$Vbeta,marred.wls$Vbeta,mcar.wls$Vbeta,mcarred.wls$Vbeta)) -c(mnarsat$min,mnarred$min,marred$min,mcar$min,mcarred$min)#log-veros. Tab.13.10 sat.lv<-sum(e134.TF2*log(e134.TF2/sum(e134.TF2))) #vlr.máx.da log-veros.de um mod.sat. #Modelos MNAR saturados podem não ter um ajuste perfeito, veja Baker e Laird (1988) ou #Poleto (2006, pp.21-26). Além disso, modelos MNAR podem ter problemas de #identificabilidade dos parâmetros, veja Poleto (2006, pp.27-30). #Poleto (2006, pp.31-40) realiza um estudo de simulação para avaliar estas 2 patologias. -2*(-c(mnarsat$min,mnarred$min,marred$min,mcar$min,mcarred$min)-sat.lv) 1-pchisq(2*(c(mnarred$min,marred$min,mcar$min,mcarred$min)+sat.lv),c(1,1,2,3)) #Valor-P c(mnarsat$code,mnarred$code,marred$code,mcar$code,mcarred$code) c(mnarsat$it,mnarred$it,marred$it,mcar$it,mcarred$it) #Tabela 13.11 matrix(mnarsat.esp(p=mnarsat$est,sum(e134.TF2))[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) matrix(mnarred.esp(p=mnarred$est,sum(e134.TF2))[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) matrix(c(e134.satmarml$yst$st1.1,e134.satmarml$yst$st1.2,e134.satmarml$yst$st1.3, 24*e134.satmarml$theta)[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) #MAR saturado matrix(marred.esp(p=marred$est,sum(e134.TF2))[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) matrix(mcar.esp(p=mcar$est,sum(e134.TF2))[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) matrix(mcarred.esp(p=mcarred$est,sum(e134.TF2))[rep(c(1,3,2,4),4)+rep(seq(0,12,4),rep(4,4))],2) #Em tabelas 2x2, homogeneidade marginal <=> simetria mnarsatHM.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1,\alpha_2 t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5];a1<-p[6];a2<-p[7] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))) )+ n121*log( t12*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t12)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (exp(a20+a1+a2)/(1+exp(a20+a1+a2))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))) )+ n22*log( t12*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t12)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))* (1/(1+exp(a20+a1+a2))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t12*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))) + (1-t11-t12-t12)*(1/(1+exp(a10+a1+a2)))* (exp(a30+a1+a2)/(1+exp(a30+a1+a2))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))) + t12*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t12)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) ) value } mnarsatHM.esp<-function(p,N){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1,\alpha_2 t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5];a1<-p[6];a2<-p[7] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10+a2)/(1+exp(a10+a2)))*(exp(a20+a2)/(1+exp(a20+a2))), t12*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t12)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))*(exp(a20+a1+a2)/(1+exp(a20+a1+a2))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10+a2)/(1+exp(a10+a2)))*(1/(1+exp(a20+a2))), t12*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), (1-t11-t12-t12)*(exp(a10+a1+a2)/(1+exp(a10+a1+a2)))*(1/(1+exp(a20+a1+a2))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10+a2)))*(exp(a30+a2)/(1+exp(a30+a2))), t12*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), (1-t11-t12-t12)*(1/(1+exp(a10+a1+a2)))*(exp(a30+a1+a2)/(1+exp(a30+a1+a2))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10+a2)))*(1/(1+exp(a30+a2))), t12*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))), (1-t11-t12-t12)*(1/(1+exp(a10+a1+a2)))*(1/(1+exp(a30+a1+a2))) ) value } mnarredHM.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5];a1<-p[6] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t12*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t12)*(exp(a10+a1)/(1+exp(a10+a1)))* (exp(a20+a1)/(1+exp(a20+a1))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t12*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t12)*(exp(a10+a1)/(1+exp(a10+a1)))* (1/(1+exp(a20+a1))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t12*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t12)*(1/(1+exp(a10+a1)))* (exp(a30+a1)/(1+exp(a30+a1))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) + (1-t11-t12-t12)*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) ) ) value } mnarredHM.esp<-function(p,N){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5];a1<-p[6] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t12)*(exp(a10+a1)/(1+exp(a10+a1)))*(exp(a20+a1)/(1+exp(a20+a1))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), (1-t11-t12-t12)*(exp(a10+a1)/(1+exp(a10+a1)))*(1/(1+exp(a20+a1))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t12*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), (1-t11-t12-t12)*(1/(1+exp(a10+a1)))*(exp(a30+a1)/(1+exp(a30+a1))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))), (1-t11-t12-t12)*(1/(1+exp(a10+a1)))*(1/(1+exp(a30+a1))) ) value } marredHM.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5];a1<-p[6] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n122*log( (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))) + (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t12)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + (1-t11-t12-t12)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) ) value } marredHM.esp<-function(p,N){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30},\alpha_1 t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5];a1<-p[6] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))), (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(exp(a20+a1)/(1+exp(a20+a1))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))), (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20+a1))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), (1-t11-t12-t12)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), (1-t11-t12-t12)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) value } mcarHM.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30} t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n122*log( (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) + (1-t11-t12-t12)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a30))) + (1-t11-t12-t12)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) ) value } mcarHM.esp<-function(p,N){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10},\alpha_{20},\alpha_{30} t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4];a30<-p[5] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t11*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t12*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), (1-t11-t12-t12)*(1/(1+exp(a10)))*(exp(a30)/(1+exp(a30))), t11*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), t12*(1/(1+exp(a10)))*(1/(1+exp(a30))), (1-t11-t12-t12)*(1/(1+exp(a10)))*(1/(1+exp(a30))) ) value } mcarredHM.mlv<-function(p,n111=e134.TF2[1],n112=e134.TF2[2],n121=e134.TF2[3],n122=e134.TF2[4], n21=e134.TF2[5],n22=e134.TF2[6],n31=e134.TF2[7],n32=e134.TF2[8],N4=e134.TF2[9]){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10}=\alpha_{30},\alpha_{20} t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4] value<- -( n111*log( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n112*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n121*log( t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n122*log( (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))) )+ n21*log( t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n22*log( t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) + (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))) )+ n31*log( t11*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) + t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) )+ n32*log( t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) + (1-t11-t12-t12)*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))) )+ N4*log( t11*(1/(1+exp(a10)))*(1/(1+exp(a10))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a10))) + t12*(1/(1+exp(a10)))*(1/(1+exp(a10))) + (1-t11-t12-t12)*(1/(1+exp(a10)))*(1/(1+exp(a10))) ) ) value } mcarredHM.esp<-function(p,N){ #p=\theta_{11},\theta_{12}=\theta_{21},\alpha_{10}=\alpha_{30},\alpha_{20} t11<-p[1];t12<-p[2] a10<-p[3];a20<-p[4] value<-N*c( t11*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(exp(a20)/(1+exp(a20))), t11*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t12*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), (1-t11-t12-t12)*(exp(a10)/(1+exp(a10)))*(1/(1+exp(a20))), t11*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), #acertar a ordem dos thetas t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), t12*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), (1-t11-t12-t12)*(1/(1+exp(a10)))*(exp(a10)/(1+exp(a10))), t11*(1/(1+exp(a10)))*(1/(1+exp(a10))), t12*(1/(1+exp(a10)))*(1/(1+exp(a10))), t12*(1/(1+exp(a10)))*(1/(1+exp(a10))), (1-t11-t12-t12)*(1/(1+exp(a10)))*(1/(1+exp(a10))) ) value } mnarsatHM<-.nlmP(objfunc=mnarsatHM.mlv,params=inipars[-3] , lower=minpars[-3] ,upper=maxpars[-3]) mnarredHM<-.nlmP(objfunc=mnarredHM.mlv,params=inipars[-c(3,8)] , lower=minpars[-c(3,8)] ,upper=maxpars[-c(3,8)]) marredHM <-.nlmP(objfunc=marredHM.mlv ,params=inipars[-c(3,8)] , lower=minpars[-c(3,8)] ,upper=maxpars[-c(3,8)]) mcarHM <-.nlmP(objfunc=mcarHM.mlv ,params=inipars[-c(3,7:8)], lower=minpars[-c(3,7:8)],upper=maxpars[-c(3,7:8)]) mcarredHM<-.nlmP(objfunc=mcarredHM.mlv,params=inipars[-c(3,6:8)], lower=minpars[-c(3,6:8)],upper=maxpars[-c(3,6:8)]) c(mnarsatHM$code,mnarredHM$code,marredHM$code,mcarHM$code,mcarredHM$code) c(mnarsatHM$it,mnarredHM$it,marredHM$it,mcarHM$it,mcarredHM$it) -c(mnarsatHM$min,mnarredHM$min,marredHM$min,mcarHM$min,mcarredHM$min) 2*(c(mnarsatHM$min,mnarredHM$min,marredHM$min,mcarHM$min,mcarredHM$min)+sat.lv) #Q_V(M,H) 1-pchisq(2*(c(mnarsatHM$min,mnarredHM$min,marredHM$min,mcarHM$min, mcarredHM$min)+sat.lv),c(1,2,2,3,4)) #Valor-P 2*(c(mnarsatHM$min,mnarredHM$min,marredHM$min,mcarHM$min,mcarredHM$min)- c(mnarsat$min,mnarred$min,marred$min,mcar$min,mcarred$min)) #Q_V(H|M) 1-pchisq(2*(c(mnarsatHM$min,mnarredHM$min,marredHM$min,mcarHM$min,mcarredHM$min)- c(mnarsat$min,mnarred$min,marred$min,mcar$min,mcarred$min)),1) #Valor-P QP<-function(esp){ esp2<-c(esp[1:4],sum(esp[5:6]),sum(esp[7:8]),sum(esp[c(9,11)]),sum(esp[c(10,12)]), sum(esp[13:16])) c(t(e134.TF2-esp2)%*%solve(diag(esp2))%*%(e134.TF2-esp2)) } QPs<-c(QP(mnarsatHM.esp(p=mnarsatHM$est,sum(e134.TF2))), QP(mnarredHM.esp(p=mnarredHM$est,sum(e134.TF2))), QP(marredHM.esp(p=marredHM$est,sum(e134.TF2))), QP(mcarHM.esp(p=mcarHM$est,sum(e134.TF2))), QP(mcarredHM.esp(p=mcarHM$est,sum(e134.TF2)))) rbind(QPs,1-pchisq(QPs,c(1,2,2,3,4))) #Q_P(M,H) / Valor-P