Description Usage Arguments Details Value Author(s) References Examples
funlinWLS
fits functional linear models by WLS (weighted least squares).
For complete data, it is based on a object of the class readCatdata
.
For missing data, it is based on a object of the class satMarML
(under MAR or MCAR) or satMcarWLS
(under MCAR).
For both complete and missing data, another alternative is to use as inputs the
maximum likelihood (ML) or any other best asymptotically normal (BAN) estimate
for the product-multinomial probabilities under a fitted model (which may encompass
MAR, MCAR or MNAR assumption) and a consistent estimate of its asymptotic covariance
matrix.
Depending on the formulation (freedom equations or constraints) and on the model
specification, different arguments must be informed.
1 2 | funlinWLS(model, obj, theta, Vtheta, A1, A2, A3, A4, A5, A6, A7, A8,
A9, X, U, XL, UL, zeroN, PI1, PI2, PI3, PI4, PI5, PI6, PI7, PI8, PI9)
|
model |
vector of functions to be modeled linearly; they need to be specified
in the order with which they are applied to the vector of proportions; the supported functions are:
linear ( |
obj |
an object of class |
theta |
BAN estimate of the product-multinomial probabilities. |
Vtheta |
consistent estimate of asymptotic covariance matrix of the estimators of theta. |
A1 |
a matrix for the 1st linear function; for linear model
( |
A2 |
a matrix for the 2nd linear function. |
A3 |
a matrix for the 3rd linear function. |
A4 |
a matrix for the 4th linear function. |
A5 |
a matrix for the 5th linear function. |
A6 |
a matrix for the 6th linear function. |
A7 |
a matrix for the 7th linear function. |
A8 |
a matrix for the 8th linear function. |
A9 |
a matrix for the 9th linear function. |
X |
a model specification matrix for freedom equation formulation for all
functional linear models; for log-linear model
( |
U |
a matrix for constraint formulation for all functional linear models;
for log-linear model ( |
XL |
a model specification matrix for freedom equation formulation of
generalized log-linear models ( |
UL |
a matrix for constraint formulation of generalized log-linear models ( |
zeroN |
for complete data, the vector has S values that are used to replace null
frequencies in each subpopulation in order to compute the proportions used in the estimate
for the covariance matrix; the default value is |
PI1 |
the 1st vector of constants to be added. |
PI2 |
the 2nd vector of constants to be added. |
PI3 |
the 3rd vector of constants to be added. |
PI4 |
the 4th vector of constants to be added. |
PI5 |
the 5th vector of constants to be added. |
PI6 |
the 6th vector of constants to be added. |
PI7 |
the 7th vector of constants to be added. |
PI8 |
the 8th vector of constants to be added. |
PI9 |
the 9th vector of constants to be added. |
Every linear function demands the specification of a matrix in Ai, where i may vary from 1 to 9; such matrices must be numbered from right to left, in the order of which the operations are applied. Similarly, the user needs to specify a vector PIi for every addition of constants.
Examples of functions (for simplicity, consider "*" as a matrix multiplication in the following functions)
Function F(Theta) | model |
A1*Theta | "lin" |
log(Theta) | "log" |
exp(Theta) | "exp" |
PI1+Theta | "add" |
A1*log(Theta) | c("lin","log") |
exp[A1*log(Theta)] | c("exp","lin","log") |
PI3+exp[PI2+A1*log(PI1+Theta)] | c("add","exp","add","lin","log","add") |
PI1+exp(A4*logA3*exp[A2*log(A1*Theta)]) | c("add","exp","lin","log","lin","exp","lin","log","lin") |
Function F(Theta) | Arguments that must be supplied |
A1*Theta | A1 |
log(Theta) | (none) |
exp(Theta) | (none) |
PI1+Theta | PI1 |
A1*log(Theta) | A1 |
exp[A1*log(Theta)] | A1 |
PI3+exp[PI2+A1*log(PI1+Theta)] | A1, PI1, PI2, PI3 |
PI1+exp(A4*logA3*exp[A2*log(A1*Theta)]) | A1, A2, A3, A4, PI1 |
Functional linear models may be fitted to the functions F(Theta) using a freedom equation formulation F(Theta)=X%*%Beta, where the elements of Beta are the parameters to be estimated, or using a constraint formulation U%*%F(Theta)=0. Both formulations lead to an equivalent model fit if U%*%X=0.
Specifically for log-linear models (model=c("lin","log")
), X and U are used
for ordinary log-linear models, and XL and UL are used for generalized log-linear models,
namely log(Theta) = nu+X%*%Beta, U%*%log(Theta) = 0, A1%*%log(Theta) = XL%*%Beta,
UL%*%A1%*%log(Theta) = 0, where nu are non-estimated parameters included only to satisfy
the natural constraints of the product-multinomial distribution.
The generic functions print
and summary
are used to print the results and to obtain a
summary thereof.
An object of the class funlinWLS
is a list containing most of the components of the
argument obj
as well as the following components:
thetaH |
vector of WLS estimates for all product-multinomial probabilities under the functional linear model (in the case of missing data, conditional on the previously assumed model for the missingness mechanism. |
VthetaH |
corresponding estimated covariance matrix. |
beta |
vector of WLS estimates for parameters of the functional linear model (only for the freedom equation formulation). |
Vbeta |
corresponding estimated covariance matrix (only for the freedom equation formulation). |
Fu |
observed functions, without model constraints. |
VFu |
corresponding estimated covariance matrix. |
FH |
WLS estimates for the functions under the fitted model. |
VFH |
corresponding estimated covariance matrix. |
QwH |
Wald statistic for testing the goodness of fit of the functional linear model (for missing data, conditional on the assumed missingness mechanism). |
glH |
degrees of freedom for testing the goodness of fit of the functional linear model (for missing data, conditional on the assumed missingness mechanism). |
ystH |
for complete data, it contains the WLS estimates for the frequencies under the functional linear model; for missing data, it contains the WLS estimates for the augmented frequencies under both the linear model and the assumed missingness mechanism; for both missing and complete data, this is computed only for linear and certain log-linear models wherein it is possible to estimate the marginal probabilities of categorization. |
Frederico Zanqueta Poleto(frederico@poleto.com)
Julio da Motta Singer (jmsinger@ime.usp.br)
Carlos Daniel Paulino (daniel.paulino@math.ist.utl.pt)
with the collaboration of
Fabio Mathias Correa (fmcorrea@uesc.br)
Enio Galinkin Jelihovschi (eniojelihovs@gmail.com)
Paulino, C.D. e Singer, J.M. (2006). Analise de dados categorizados (in Portuguese). Sao Paulo: Edgard Blucher.
Poleto, F.Z. (2006). Analise de dados categorizados com omissao (in Portuguese). Dissertacao de mestrado. IME-USP. http://www.poleto.com/missing.html.
Poleto, F.Z., Singer, J.M. e Paulino, C.D. (2007). Analyzing categorical data with complete or missing responses using the Catdata package. Unpublished vignette. http://www.poleto.com/missing.html.
Poleto, F.Z., Singer, J.M. e Paulino, C.D. (2012). A product-multinomial framework for categorical data analysis with missing responses. To appear in Brazilian Journal of Probability and Statistics. http://imstat.org/bjps/papers/BJPS198.pdf.
Singer, J. M., Poleto, F. Z. and Paulino, C. D. (2007). Catdata: software for analysis of categorical data with complete or missing responses. Actas de la XII Reunion Cientifica del Grupo Argentino de Biometria y I Encuentro Argentino-Chileno de Biometria. http://www.poleto.com/SingerPoletoPaulino2007GAB.pdf.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 | #Example 11.2 of Paulino and Singer (2006)
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))
#Two equivalent ways of fitting the same symmetry model
e112.linwls1<-funlinWLS(model="lin",obj=e112.catdata,U=e112.U)
e112.linwls2<-funlinWLS(model="lin",obj=e112.catdata,X=e112.X)
e112.linwls1 #constraint formulation
e112.linwls2 #freedom equation formulation
summary(e112.linwls1)
#Example 11.5 of Paulino and Singer (2006)
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)
e115.loglinwls2<-funlinWLS(model=c("lin", "log"), obj=e115.catdata,
X=e115.X)
e115.loglinwls3<-funlinWLS(model=c("lin", "log"), obj=e115.catdata,
X=e115.X2)
e115.loglinwls4<-funlinWLS(model=c("lin", "log"), obj=e115.catdata,
A1=c(1,-1,-1,1), XL=1)
#Independence ordinary log-linear model, constraint formulation
e115.loglinwls1
#Independence ordinary log-linear model, freedom equation formulation
e115.loglinwls2
#Saturated ordinary log-linear model, freedom equation formulation
e115.loglinwls3
#Saturated generalized log-linear model, freedom equation formulation
e115.loglinwls4
#95% confidence interval for log-odds ratio and for odds ratio
round(e115.loglinwls4$beta+c(-1,1)*qnorm(0.975)*sqrt(e115.loglinwls4$Vbeta),3)
round(exp(e115.loglinwls4$beta),3)
round(exp(e115.loglinwls4$beta+c(-1,1)*qnorm(0.975)*sqrt(e115.loglinwls4$Vbeta)),3)
#Example 11.3 of Paulino and Singer (2006)
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)
e113.linwls2<-funlinWLS(model="lin",obj=e113.catdata,X=e113.X)
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)
e113.linwls4<-funlinWLS(model="lin",obj=e113.catdata,A1=e113.A,X=e113.X2)
#Four equivalent ways of fitting the same marginal homogeneity model
e113.linwls1
e113.linwls2
e113.linwls3
e113.linwls4
#Example 11.12 of Paulino and Singer (2006)
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),
tcrossprod(-c(2,1),(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)
# confidence interval
round(e1112.kappa$beta+c(-1,1)*qnorm(0.975)*sqrt(e1112.kappa$Vbeta),3)
#weighted kappa (Spitzer, Cohen, Fleiss e Endicott, 1967)
#squared weights (Fleiss e Cohen, 1973)
W1<-c(1,0.75,0,0.75,1,0.75,0,0.75,1)
#absolute weights (Cicchetti e Allison, 1971)
W2<-c(1,0.5,0,0.5,1,0.5,0,0.5,1)
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)
#Example 1 of Poleto et al (2012)
smoking.TF<-rbind(c(167,17,19,10,1,3,52,10,11, 176,24,121, 28,10,12),
c(120,22,19, 8,5,1,39,12,12, 103, 3, 80, 31, 8,14))
smoking.Zp<-kronecker(t(rep(1,2)),cbind(kronecker(diag(3),rep(1,3)),
kronecker(rep(1,3),diag(3))))
smoking.Rp<-rbind(c(3,3),c(3,3))
smoking.catdata<-readCatdata(TF=smoking.TF,Zp=smoking.Zp,Rp=smoking.Rp)
smoking.catdata #Proportions of the complete data
smoking.satmarml<-satMarML(smoking.catdata)
smoking.satmcarml<-satMarML(smoking.catdata,missing="MCAR")
smoking.satmcarwls<-satMcarWLS(smoking.catdata)
smoking.E<-rbind(c(1,-1,0),c(0,1,-1))
smoking.A<-kronecker(kronecker(diag(2),smoking.E),smoking.E)
smoking.loglin2.marhybrid<-funlinWLS(model=c("lin","log"),
obj=smoking.satmarml, A1=smoking.A, XL=rep(1,8))
smoking.loglin2.mcarhybrid<-funlinWLS(model=c("lin","log"),
obj=smoking.satmcarml, A1=smoking.A, XL=rep(1,8))
smoking.loglin2.mcarwls<-funlinWLS(model=c("lin","log"),
obj=smoking.satmcarwls, A1=smoking.A, XL=rep(1,8))
#MNAR example
#Minus log-likelihood of the MNAR described in the last paragraph of Section 3
mnar.mll<-function(p,fr){
# p=(pi11(1),...,pi33(2),a2(11),a2(21),a2(31),a3(11),
# a3(21),a3(31),a2(12),a2(22),
# a2(32),a3(12),a3(22),a3(32))
pi11.1<-p[1]; pi12.1<-p[2]; pi13.1<-p[3]
pi21.1<-p[4]; pi22.1<-p[5]; pi23.1<-p[6]
pi31.1<-p[7]; pi32.1<-p[8]
pi33.1=1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1
pi11.2<-p[9]; pi12.2<-p[10];pi13.2<-p[11]
pi21.2<-p[12];pi22.2<-p[13];pi23.2<-p[14]
pi31.2<-p[15];pi32.2<-p[16]
pi33.2=1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2
a2.1.1<-p[17];a2.2.1<-p[18];a2.3.1<-p[19]
a3.1.1<-p[20];a3.2.1<-p[21];a3.3.1<-p[22]
a2.1.2<-p[23];a2.2.2<-p[24];a2.3.2<-p[25]
a3.1.2<-p[26];a3.2.2<-p[27];a3.3.2<-p[28]
value<- -(
fr[1,1]*log(pi11.1*(1-a2.1.1-a3.1.1))+fr[1,2]*log(pi12.1*(1-a2.2.1-a3.1.1))+
fr[1,3]*log(pi13.1*(1-a2.3.1-a3.1.1))+
fr[1,4]*log(pi21.1*(1-a2.1.1-a3.2.1))+fr[1,5]*log(pi22.1*(1-a2.2.1-a3.2.1))+
fr[1,6]*log(pi23.1*(1-a2.3.1-a3.2.1))+
fr[1,7]*log(pi31.1*(1-a2.1.1-a3.3.1))+fr[1,8]*log(pi32.1*(1-a2.2.1-a3.3.1))+
fr[1,9]*log(pi33.1*(1-a2.3.1-a3.3.1))+
fr[1,10]*log(pi11.1*a2.1.1 + pi12.1*a2.2.1 + pi13.1*a2.3.1)+
fr[1,11]*log(pi21.1*a2.1.1 + pi22.1*a2.2.1 + pi23.1*a2.3.1)+
fr[1,12]*log(pi31.1*a2.1.1 + pi32.1*a2.2.1 + pi33.1*a2.3.1)+
fr[1,13]*log(pi11.1*a3.1.1 + pi21.1*a3.2.1 + pi31.1*a3.3.1)+
fr[1,14]*log(pi12.1*a3.1.1 + pi22.1*a3.2.1 + pi32.1*a3.3.1)+
fr[1,15]*log(pi13.1*a3.1.1 + pi23.1*a3.2.1 + pi33.1*a3.3.1)+
fr[2,1]*log(pi11.2*(1-a2.1.2-a3.1.2))+fr[2,2]*log(pi12.2*(1-a2.2.2-a3.1.2))+
fr[2,3]*log(pi13.2*(1-a2.3.2-a3.1.2))+
fr[2,4]*log(pi21.2*(1-a2.1.2-a3.2.2))+fr[2,5]*log(pi22.2*(1-a2.2.2-a3.2.2))+
fr[2,6]*log(pi23.2*(1-a2.3.2-a3.2.2))+
fr[2,7]*log(pi31.2*(1-a2.1.2-a3.3.2))+fr[2,8]*log(pi32.2*(1-a2.2.2-a3.3.2))+
fr[2,9]*log(pi33.2*(1-a2.3.2-a3.3.2))+
fr[2,10]*log(pi11.2*a2.1.2 + pi12.2*a2.2.2 + pi13.2*a2.3.2)+
fr[2,11]*log(pi21.2*a2.1.2 + pi22.2*a2.2.2 + pi23.2*a2.3.2)+
fr[2,12]*log(pi31.2*a2.1.2 + pi32.2*a2.2.2 + pi33.2*a2.3.2)+
fr[2,13]*log(pi11.2*a3.1.2 + pi21.2*a3.2.2 + pi31.2*a3.3.2)+
fr[2,14]*log(pi12.2*a3.1.2 + pi22.2*a3.2.2 + pi32.2*a3.3.2)+
fr[2,15]*log(pi13.2*a3.1.2 + pi23.2*a3.2.2 + pi33.2*a3.3.2)
)
value
}
mnar.fit<-constrOptim(theta=c(rep(1/9,16), rep(1/3,12)), f=mnar.mll,
method="Nelder-Mead", ui=diag(28), ci=rep(0,28),
control=list(maxit=10000), outer.iterations=1000, fr=smoking.TF)
#hessian matrix
mnar.der<-deriv3(~-(
o1.1*log(pi11.1*(1-a2.1.1-a3.1.1))+o1.2*log(pi12.1*(1-a2.2.1-a3.1.1))+
o1.3*log(pi13.1*(1-a2.3.1-a3.1.1))+
o1.4*log(pi21.1*(1-a2.1.1-a3.2.1))+o1.5*log(pi22.1*(1-a2.2.1-a3.2.1))+
o1.6*log(pi23.1*(1-a2.3.1-a3.2.1))+
o1.7*log(pi31.1*(1-a2.1.1-a3.3.1))+o1.8*log(pi32.1*(1-a2.2.1-a3.3.1))+
o1.9*log((1-pi11.1-pi12.1-pi13.1-pi21.1-
pi22.1-pi23.1-pi31.1-pi32.1)*(1-a2.3.1-a3.3.1))+
o1.10*log(pi11.1*a2.1.1 + pi12.1*a2.2.1 + pi13.1*a2.3.1)+
o1.11*log(pi21.1*a2.1.1 + pi22.1*a2.2.1 + pi23.1*a2.3.1)+
o1.12*log(pi31.1*a2.1.1 + pi32.1*a2.2.1 +
(1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*a2.3.1)+
o1.13*log(pi11.1*a3.1.1 + pi21.1*a3.2.1 + pi31.1*a3.3.1)+
o1.14*log(pi12.1*a3.1.1 + pi22.1*a3.2.1 + pi32.1*a3.3.1)+
o1.15*log(pi13.1*a3.1.1 + pi23.1*a3.2.1 +
(1-pi11.1-pi12.1-pi13.1-pi21.1-pi22.1-pi23.1-pi31.1-pi32.1)*a3.3.1)+
o2.1*log(pi11.2*(1-a2.1.2-a3.1.2))+o2.2*log(pi12.2*(1-a2.2.2-a3.1.2))+
o2.3*log(pi13.2*(1-a2.3.2-a3.1.2))+
o2.4*log(pi21.2*(1-a2.1.2-a3.2.2))+o2.5*log(pi22.2*(1-a2.2.2-a3.2.2))+
o2.6*log(pi23.2*(1-a2.3.2-a3.2.2))+
o2.7*log(pi31.2*(1-a2.1.2-a3.3.2))+o2.8*log(pi32.2*(1-a2.2.2-a3.3.2)) +
o2.9*log((1-pi11.2-pi12.2-pi13.2-pi21.2-
pi22.2-pi23.2-pi31.2-pi32.2)*(1-a2.3.2-a3.3.2))+
o2.10*log(pi11.2*a2.1.2 + pi12.2*a2.2.2 + pi13.2*a2.3.2)+
o2.11*log(pi21.2*a2.1.2 + pi22.2*a2.2.2 + pi23.2*a2.3.2)+
o2.12*log(pi31.2*a2.1.2 + pi32.2*a2.2.2 +
(1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*a2.3.2)+
o2.13*log(pi11.2*a3.1.2 + pi21.2*a3.2.2 + pi31.2*a3.3.2)+
o2.14*log(pi12.2*a3.1.2 + pi22.2*a3.2.2 + pi32.2*a3.3.2)+
o2.15*log(pi13.2*a3.1.2 + pi23.2*a3.2.2 +
(1-pi11.2-pi12.2-pi13.2-pi21.2-pi22.2-pi23.2-pi31.2-pi32.2)*a3.3.2)
),c("pi11.1","pi12.1","pi13.1","pi21.1","pi22.1","pi23.1","pi31.1","pi32.1",
"pi11.2","pi12.2","pi13.2","pi21.2","pi22.2","pi23.2","pi31.2","pi32.2",
"a2.1.1","a2.2.1","a2.3.1","a3.1.1","a3.2.1","a3.3.1","a2.1.2","a2.2.2",
"a2.3.2","a3.1.2","a3.2.2","a3.3.2"),
c("pi11.1","pi12.1","pi13.1","pi21.1","pi22.1","pi23.1","pi31.1","pi32.1",
"pi11.2","pi12.2","pi13.2","pi21.2","pi22.2","pi23.2","pi31.2","pi32.2",
"a2.1.1","a2.2.1","a2.3.1","a3.1.1","a3.2.1","a3.3.1","a2.1.2","a2.2.2",
"a2.3.2","a3.1.2","a3.2.2","a3.3.2",
"o1.1","o1.2","o1.3","o1.4","o1.5","o1.6","o1.7","o1.8","o1.9","o1.10",
"o1.11","o1.12","o1.13","o1.14","o1.15",
"o2.1","o2.2","o2.3","o2.4","o2.5","o2.6","o2.7","o2.8","o2.9","o2.10",
"o2.11","o2.12","o2.13","o2.14","o2.15")
)
p<-mnar.fit$par;TF<-smoking.TF
mnar.InfObs<-mnar.der(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9],p[10],
p[11],p[12],p[13],p[14], p[15],p[16],p[17],p[18],p[19],p[20],p[21],
p[22],p[23],p[24],p[25],p[26],p[27],p[28], TF[1,1],TF[1,2],TF[1,3],
TF[1,4],TF[1,5],TF[1,6],TF[1,7],TF[1,8],TF[1,9],TF[1,10],
TF[1,11],TF[1,12],TF[1,13],TF[1,14],TF[1,15],
TF[2,1],TF[2,2],TF[2,3],TF[2,4],TF[2,5],TF[2,6],TF[2,7],TF[2,8],TF[2,9],TF[2,10],
TF[2,11],TF[2,12],TF[2,13],TF[2,14],TF[2,15])
b<-smoking.catdata$b #b in (8), i.e., rep(1,2)%x%c(rep(0,8),1)
B<-smoking.catdata$B #B in (8), i.e., diag(2)%x%rbind(diag(8),rep(-1,8))
smoking.loglin2mnar.hybrid<-funlinWLS(model=c("lin","log"),
theta=as.vector(b+c(B%*%mnar.fit$par[1:16])),
Vtheta=B%*%solve(attr(mnar.InfObs,"hessian")[1,,])[1:16,1:16]%*%t(B),
A1=smoking.A,X=rep(1,8))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.