1##logit and invlogit transformation 2logit <- function(X) 3 { 4 Y<-log(X/(1-X)) 5 Y 6 } 7 8invlogit <-function(Y) 9 { 10 X<-exp(Y)/(1+exp(Y)) 11 X 12 } 13 14####assuming theta.em 15##2 d: mu1, mu2, sig1, sig2, r12 16##3 d: mu3, mu1, mu2, sig3, sig1, sig2, r13, r23, r12 17param.pack<-function(theta.em, fix.rho=FALSE,r12=0, dim) 18 { 19 mu<-rep(0, dim) 20 Sig<-matrix(0,dim, dim) 21 22 mu<-theta.em[1:dim] 23 24 for (i in 1:dim) 25 Sig[i,i]<-theta.em[dim+i] 26 27 if (!fix.rho) { 28 Sig[1,2]<-Sig[2,1]<-theta.em[2*dim+1]*sqrt(Sig[1,1]*Sig[2,2]) 29 if (dim==3) { 30 Sig[1,3]<-Sig[3,1]<-theta.em[2*dim+2]*sqrt(Sig[1,1]*Sig[3,3]) 31 Sig[2,3]<-Sig[3,2]<-theta.em[2*dim+3]*sqrt(Sig[2,2]*Sig[3,3]) 32 } 33 } 34 if (fix.rho) { 35 if (dim==2) 36 Sig[1,2]<-Sig[2,1]<-r12*sqrt(Sig[1,1]*Sig[2,2]) 37 if (dim==3) { 38 Sig[1,2]<-Sig[2,1]<-theta.em[2*dim+1]*sqrt(Sig[1,1]*Sig[2,2]) 39 Sig[1,3]<-Sig[3,1]<-theta.em[2*dim+2]*sqrt(Sig[1,1]*Sig[3,3]) 40 Sig[2,3]<-Sig[3,2]<-r12*sqrt(Sig[2,2]*Sig[3,3]) 41 } 42 } 43 return(list(mu=mu, Sigma=Sig)) 44} 45 46## transformation of BVN parameter into 47## Fisher scale or unit scale 48## in 2 D, mu1, mu2, sigma1, sigma2, r12 49## in 3 D, mu3, mu1, mu2, sigma3, sigma1, sigma2, sigma31, sigma32, sigma12 50param.trans <-function(X, transformation="Fisher") { 51 p<-length(X) 52 Y<-rep(0,p) 53 54 if (transformation=="Fisher") { 55 if (p<=5) { 56 Y[1:2]<-X[1:2] 57 Y[3:4]<-log(X[3:4]) 58 if (p==5) 59 Y[5]<-0.5*log((1+X[5])/(1-X[5])) 60 } 61 62 if (p>5) { 63 Y[1:3]<-X[1:3] 64 Y[4:6]<-log(X[4:6]) 65 Y[7:8]<-0.5*log((1+X[7:8])/(1-X[7:8])) 66 if (p==9) 67 Y[9]<-0.5*log((1+X[9])/(1-X[9])) 68 } 69 } 70 71 if (transformation=="unitscale") { 72 if (p<=5) { 73 Y[1:2] <- invlogit(X[1:2]) 74 Y[3:4] <- X[3:4]*exp(2*X[1:2])/(1+exp(X[1:2]))^4 75 if (p==5) 76 Y[5] <- X[5] 77 } 78 79 if (p>5) { 80 Y[1:3]<-invlogit(X[1:3]) 81 Y[4:6]<-X[4:6]*exp(2*X[4:6])/(1+exp(X[4:6]))^4 82 Y[7:8]<-X[7:8] 83 if (p==9) 84 Y[9]<-X[9] 85 } 86 } 87 return(Y) 88} 89 90vec<-function(mat) { 91 v<-as.vector(mat, mode="any") 92 v 93} 94 95tr<-function(mat) { 96 trace<-sum(diag(mat)) 97 trace 98} 99## I_{com} 100## the gradient function for multivariate normal function 101#du.theta and dSig.theta are the first derivative of mu and Sigma 102#with respect to theta 103#du.theta[n.u, n.theta] 104#dSig.theta[n.u, n.u, n.theta] 105 106d1st.mvn<-function(mu,Sigma, fix.rho=FALSE) { 107 #r12, r13,r23 are internal here, 108 # r12 doesn't correspond to cor(w1, w2) in 3d case (intead, r12=>cor(W1,x) 109 d<-length(mu) 110 p<-d+d+d*(d-1)/2 111 u1<-mu[1] 112 u2<-mu[2] 113 s1<-Sigma[1,1] 114 s2<-Sigma[2,2] 115 r12<-Sigma[1,2]/sqrt(s1*s2) 116 117 if (d==3) { 118 u3<-mu[3] 119 s3<-Sigma[3,3] 120 r13<-Sigma[1,3]/sqrt(s1*s3) 121 r23<-Sigma[2,3]/sqrt(s2*s3) 122 } 123 124 if (fix.rho) p<-p-1 125 126 du.theta<-matrix(0,d,p) 127 for (j in 1:d) 128 du.theta[j,j]<-1 129 130 dSig.theta<-array(0,c(d,d,p)) 131 132 for (i in 1:d) 133 dSig.theta[i,i,d+i]<-1 134 135 dSig.theta[1,2,d+1]<-dSig.theta[2,1,d+1]<-1/2*s1^(-1/2)*s2^(1/2)*r12 136 dSig.theta[1,2,d+2]<-dSig.theta[2,1,d+2]<-1/2*s2^(-1/2)*s1^(1/2)*r12 137 if (d==3) { 138 dSig.theta[1,3,d+1]<-dSig.theta[3,1,d+1]<-1/2*s1^(-1/2)*s3^(1/2)*r13 139 dSig.theta[1,3,d+3]<-dSig.theta[3,1,d+3]<-1/2*s3^(-1/2)*s1^(1/2)*r13 140 dSig.theta[2,3,d+2]<-dSig.theta[3,2,d+2]<-1/2*s2^(-1/2)*s3^(1/2)*r23 141 dSig.theta[2,3,d+3]<-dSig.theta[3,2,d+3]<-1/2*s3^(-1/2)*s2^(1/2)*r23 142 } 143 144 if (!fix.rho) { 145 dSig.theta[1,2,2*d+1]<-dSig.theta[2,1,2*d+1]<-sqrt(s1*s2) 146 if (d==3) { 147 dSig.theta[1,3,2*d+2]<-dSig.theta[3,1,2*d+2]<-sqrt(s1*s3) 148 dSig.theta[2,3,2*d+3]<-dSig.theta[3,2,2*d+3]<-sqrt(s2*s3) 149 } 150 } 151 if (fix.rho) { 152 if (d==3) { 153 dSig.theta[1,3,2*d+1]<-dSig.theta[3,1,2*d+1]<-sqrt(s1*s3) 154 dSig.theta[2,3,2*d+2]<-dSig.theta[3,2,2*d+2]<-sqrt(s2*s3) 155 } 156 } 157 158 return(list(du.theta=du.theta, dSig.theta=dSig.theta)) 159} 160 161d2nd.mvn<-function(mu,Sigma, fix.rho=FALSE) { 162 #r12, r13,r23 are internal here, 163 # r12 doesn't correspond to cor(w1, w2) in 3d case (intead, r12=>cor(W1,x) 164 d<-length(mu) 165 p<-d+d+d*(d-1)/2 166 u1<-mu[1] 167 u2<-mu[2] 168 s1<-Sigma[1,1] 169 s2<-Sigma[2,2] 170 r12<-Sigma[1,2]/sqrt(s1*s2) 171 if (d==3) { 172 u3<-mu[3] 173 s3<-Sigma[3,3] 174 r13<-Sigma[1,3]/sqrt(s1*s3) 175 r23<-Sigma[2,3]/sqrt(s2*s3) 176 } 177 178 if (fix.rho) p<-p-1 179 180 ddu.theta<-array(0,c(d,p,p)) 181 182 ddSig.theta<-array(0,c(d,d,p,p)) 183 184 ddSig.theta[1,2,d+1,d+1]<-ddSig.theta[2,1,d+1,d+1]<- -1/4*s1^(-3/2)*s2^(1/2)*r12 185 ddSig.theta[1,2,d+1,d+2]<-ddSig.theta[2,1,d+1,d+2]<- 1/4*s1^(-1/2)*s2^(-1/2)*r12 186 ddSig.theta[1,2,d+2,d+2]<-ddSig.theta[2,1,d+2,d+2]<- -1/4*s1^(1/2)*s2^(-3/2)*r12 187 188 if (d==3) { 189 ddSig.theta[1,3,d+1,d+1]<-ddSig.theta[3,1,d+1,d+1]<- -1/4*s1^(-3/2)*s3^(1/2)*r13 190 ddSig.theta[1,3,d+1,d+3]<-ddSig.theta[3,1,d+1,d+3]<- 1/4*s1^(-1/2)*s3^(-1/2)*r13 191 192 ddSig.theta[2,3,d+2,d+2]<-ddSig.theta[3,2,d+2,d+2]<- -1/4*s2^(-3/2)*s3^(1/2)*r23 193 ddSig.theta[2,3,d+2,d+3]<-ddSig.theta[3,2,d+2,d+3]<- 1/4*s2^(-1/2)*s3^(-1/2)*r23 194 195 ddSig.theta[1,3,d+3,d+3]<-ddSig.theta[3,1,d+3,d+3]<- -1/4*s1^(1/2)*s3^(-3/2)*r13 196 ddSig.theta[2,3,d+3,d+3]<-ddSig.theta[3,2,d+3,d+3]<- -1/4*s2^(1/2)*s3^(-3/2)*r23 197 198 } 199 200 if (!fix.rho) { 201 ddSig.theta[1,2,d+1,2*d+1]<-ddSig.theta[2,1,d+1,2*d+1]<- 1/2*s1^(-1/2)*s2^(1/2) 202 ddSig.theta[1,2,d+2,2*d+1]<-ddSig.theta[2,1,d+2,2*d+1]<- 1/2*s1^(1/2)*s2^(-1/2) 203 204 if (d==3) { 205 ddSig.theta[1,3,d+1,2*d+2]<-ddSig.theta[3,1,d+1,2*d+2]<- 1/2*s1^(-1/2)*s3^(1/2) 206 ddSig.theta[2,3,d+2,2*d+3]<-ddSig.theta[3,2,d+2,2*d+3]<- 1/2*s2^(-1/2)*s3^(1/2) 207 ddSig.theta[1,3,d+3,2*d+2]<-ddSig.theta[3,1,d+3,2*d+2]<- 1/2*s1^(1/2)*s3^(-1/2) 208 ddSig.theta[2,3,d+3,2*d+3]<-ddSig.theta[3,2,d+3,2*d+3]<- 1/2*s2^(1/2)*s3^(-1/2) 209 } 210 } 211 if (fix.rho) { 212 if (d==3) { 213 214 ddSig.theta[1,2,d+1,2*d+1]<-ddSig.theta[2,1,d+1,2*d+1]<- 1/2*s1^(-1/2)*s3^(1/2) 215 ddSig.theta[2,3,d+2,2*d+2]<-ddSig.theta[3,2,d+2,2*d+2]<- 1/2*s2^(-1/2)*s3^(1/2) 216 ddSig.theta[1,3,d+3,2*d+1]<-ddSig.theta[3,1,d+3,2*d+1]<- 1/2*s1^(1/2)*s3^(-1/2) 217 ddSig.theta[2,3,d+3,2*d+2]<-ddSig.theta[3,2,d+3,2*d+2]<- 1/2*s2^(1/2)*s3^(-1/2) 218 } 219 } 220 221 for (i in 1:(p-1)) 222 for (j in (i+1):p) { 223 ddSig.theta[,,j,i]<-ddSig.theta[,,i,j] 224 ddu.theta[,j,i]<-ddu.theta[,i,j] 225 } 226 227 return(list(ddu.theta=ddu.theta, ddSig.theta=ddSig.theta)) 228} 229 230##assuming the order of sufficient statistics 231## 2d, mean(W1), mean(W2), mean(W1^2) mean(W2^2), mean(W1W2) 232## 3d, mean(X), mean(W1), mean(W2), mean(X^2),mean(W1^2) mean(W2^2), 233## mean(XW1), mean(XW2), mean(W1W2) 234 235suff<-function(mu, suff.stat,n) { 236 237 d<-length(mu) 238 p<-d+d+d*(d-1)/2 239 u1<-mu[1] 240 u2<-mu[2] 241 if (d==3) u3<-mu[3] 242 243 S1<-n*suff.stat[1] 244 S2<-n*suff.stat[2] 245 S11<-n*suff.stat[d+1] 246 S22<-n*suff.stat[d+2] 247 S12<-n*suff.stat[2*d+1] 248 if (d==3) { 249 S3<-n*suff.stat[d] 250 S33<-n*suff.stat[2*d] 251 S13<-n*suff.stat[2*d+2] 252 S23<-n*suff.stat[2*d+3] 253 } 254 255 Vv<-rep(0,d) 256 Vv[1]<-S1-n*u1 257 Vv[2]<-S2-n*u2 258 if (d==3) Vv[3]<-S3-n*u3 259 260 Ss<-matrix(0,d,d) 261 Ss[1,1]<-S11-2*S1*u1+n*u1^2 262 Ss[2,2]<-S22-2*S2*u2+n*u2^2 263 Ss[1,2]<-Ss[2,1]<-S12-S1*u2-S2*u1+n*u1*u2 264 if (d==3) { 265 Ss[3,3]<-S33-2*S3*u3+n*u3^2 266 Ss[1,3]<-Ss[3,1]<-S13-S1*u3-S3*u1+n*u1*u3 267 Ss[2,3]<-Ss[3,2]<-S23-S3*u2-S2*u3+n*u2*u3 268 } 269 return(list(Ss=Ss, Vv=Vv)) 270} 271 272 273 274#du.theta and dSig.theta are the second derivative of mu and Sigma 275#with respect to theta 276#ddu.theta[n.u, n.theta, n.theta] 277#ddSig.theta[n.u, n.u, n.theta, n.theta] 278 279##comput the gradient vector (expected first derivatives) for MVN 280##not actually used here. 281 282Dcom.mvn<-function(mu, Sigma, suff.stat,n, fix.rho=FALSE) { 283 d<-dim(Sigma)[1] 284 p<-d*2+0.5*d*(d-1) 285 286 if (fix.rho) { 287 p<-p-1 288 } 289 290 Dcom<-rep(0,p) 291 invSigma<-solve(Sigma) 292 293 temp<-suff(mu, suff.stat, n) 294 Ss<-temp$Ss 295 Vv<-temp$Vv 296 297 temp<-d1st.mvn(mu=mu, Sigma=Sigma, fix.rho=fix.rho) 298 du.theta<-temp$du.theta 299 dSig.theta<-temp$dSig.theta 300 301 for (i in 1:p) 302 Dcom[i]<- -n/2*t(vec(invSigma))%*%vec(dSig.theta[,,i])+ 0.5*tr(invSigma%*%dSig.theta[,,i]%*%invSigma%*%Ss)+ t(du.theta[,i])%*%invSigma%*%Vv 303 304 Dcom 305} 306 307 308#compute the information matrix of MVN 309# -1*second derivatives 310 311Icom.mvn<-function(mu, Sigma, suff.stat,n, fix.rho=FALSE) { 312 d<-dim(Sigma)[1] 313 p<-d*2+1/2*d*(d-1) 314 315 if (fix.rho) 316 { 317 p<-p-1 318 } 319 320 Icom<-matrix(0,p,p) 321 322 invSigma<-solve(Sigma) 323 324 temp<-suff(mu, suff.stat, n) 325 Ss<-temp$Ss 326 Vv<-temp$Vv 327 328 temp<-d1st.mvn(mu, Sigma, fix.rho) 329 du.theta<-temp$du.theta 330 dSig.theta<-temp$dSig.theta 331 332 temp<-d2nd.mvn(mu, Sigma, fix.rho) 333 ddu.theta<-temp$ddu.theta 334 ddSig.theta<-temp$ddSig.theta 335 336 for (i in 1:p) { 337 dinvSig.theta.i<- -invSigma%*%dSig.theta[,,i]%*%invSigma 338 for (j in 1:i) { 339 dinvSig.theta.j<- -invSigma%*%dSig.theta[,,j]%*%invSigma 340 ddinvSig.theta.ij<- -dinvSig.theta.j%*%dSig.theta[,,i]%*%invSigma -invSigma%*%ddSig.theta[,,i,j]%*%invSigma-invSigma%*%dSig.theta[,,i]%*%dinvSig.theta.j 341 342 a1<- -n/2*(t(vec(dinvSig.theta.j))%*%vec(dSig.theta[,,i]) + t(vec(invSigma))%*%vec(ddSig.theta[,,i,j])) 343 344 a2<- t(du.theta[,j])%*%dinvSig.theta.i%*%Vv - 0.5*tr(ddinvSig.theta.ij%*%Ss) 345 346 a3<- t(ddu.theta[,i,j])%*%invSigma%*%Vv + t(du.theta[,i])%*%dinvSig.theta.j%*%Vv - n*t(du.theta[,i])%*%invSigma%*%du.theta[,j] 347 348 Icom[i,j]<-a1+a2+a3 349 350 if (i!=j) Icom[j,i]<-Icom[i,j] 351 } 352 } 353 -Icom 354} 355 356 357###compute the information matrix for various parameter transformation 358### "Fisher" transformation (variance stablization?) 359### unit scale transformation: first order approximation of mean and var, rho 360 361##express T1 and T2 in more general form 362 363Icom.transform<-function(Icom, Dvec, theta, transformation="Fisher", context, fix.rho) 364 { 365 366 if (!context) { 367 368 mu<-theta[1:2] 369 sigma<-theta[3:4] 370 rho<-theta[5] 371 } 372 if (context) { 373 374 mu<-theta[1:3] # x,w1,w2 375 sigma<-theta[4:6] #x, w1, w2 376 rho<-theta[7:9] #r_xw1, r_xw2, r_w1w2 377 } 378 379 ##T1: d(theta)/d(f(theta)), theta is the MVN parameterization 380 ##T2, d2(theta)/d(f(theta))(d(f(theta))') 381 382 ### transformation=Fisher, Icom_normal==>Icom_fisher 383 384 Imat<- -Icom 385 n.par<-dim(Imat)[1] 386 387 if (transformation=="Fisher") { 388 if (!context) { 389 T1<-c(1,1,sigma[1], sigma[2]) 390 391 T2<-matrix(0, n.par^2, n.par) 392 T2[(2*n.par+3), 3]<-sigma[1] 393 T2[(3*n.par+4), 4]<-sigma[2] 394 395 if (!fix.rho) { 396 T1<-c(T1, (1-(rho[1]^2))) 397 T2[(4*n.par+5),5]<- -2*rho[1]*(1-rho[1]^2) 398 } 399 400 T1<-diag(T1) 401 } 402 403 if (context) { 404 T1<-c(1,1,1,sigma[1:3],(1-(rho[1:2]^2))) 405 406 T2<-matrix(0, n.par^2, n.par) 407 T2[(3*n.par+4), 4]<-sigma[1] 408 T2[(4*n.par+5), 5]<-sigma[2] 409 T2[(5*n.par+6), 6]<-sigma[3] 410 T2[(6*n.par+7),7]<- -2*rho[1]*(1-rho[1]^2) 411 T2[(7*n.par+8),8]<- -2*rho[2]*(1-rho[2]^2) 412 413 if (!fix.rho) { 414 T1<-c(T1, (1-(rho[3]^2))) 415 T2[(8*n.par+9),9]<- -2*rho[3]*(1-rho[3]^2) 416 } 417 418 T1<-diag(T1) 419 } 420} 421 ### transformation=unitscale, Icom_normal==>Icom_unitscale 422 if (transformation=="unitscale") { 423 424 T1<-matrix(0,n.par,n.par) 425 T1[1,1]<-exp(-mu[1])*(1+exp(mu[1]))^2 426 T1[1,3]<-1/(sigma[1]*2*exp(2*mu[1])*(1+exp(mu[1]))^(-4)*(1-2*(1+exp(mu[1]))^(-1))) 427 428 T1[2,2]<-exp(-mu[2])*(1+exp(mu[2]))^2 429 T1[2,4]<-1/(sigma[2]*2*exp(2*mu[2])*(1+exp(mu[2]))^(-4)*(1-2*(1+exp(mu[2]))^(-1))) 430 431 T1[3,3]<-2*sigma[1]^0.5*(1+exp(mu[1]))^4*exp(-2*mu[1]) 432 T1[4,4]<-2*sigma[2]^0.5*(1+exp(mu[2]))^4*exp(-2*mu[2]) 433 434 435 # T2<-matrix(0, n.par^2, n.par) 436 # T2[1,1]<- 437 # T2[(1*n.par+2), (1*n.par+2)]<- 438 439 440 ##compute T1 and T2 441 442 } 443 444 Icom.tran<-matrix(NA, n.par, n.par) 445 Icom.tran<-T1%*%Imat%*%t(T1) 446 447 temp1<-matrix(0,n.par,n.par) 448 for (i in 1:n.par) 449 for (j in 1:n.par) 450 temp1[i,j]<- Dvec%*%T2[((i-1)*n.par+(1:n.par)),j] 451 452 Icom.tran<-Icom.tran+temp1 453 return(-Icom.tran) 454} 455 456 457ecoINFO<-function(theta.em, suff.stat, DM, context=TRUE, fix.rho=FALSE, sem=TRUE, r12=0, n) 458 { 459 460 if (context) fix.rho<-FALSE 461 ndim<-2 462 if (context) ndim<-3 463 464 n.var<-2*ndim+ ndim*(ndim-1)/2 465 466 n.par<-n.var 467 if (context) { 468 n.par<-n.var-2 469 } 470 471 if (!context & fix.rho) n.par<-n.par-1 472 473 mu<-param.pack(theta.em, fix.rho=fix.rho, r12=r12, dim=ndim)$mu 474 Sigma<-param.pack(theta.em, fix.rho=fix.rho, r12=r12, dim=ndim)$Sigma 475 476 theta.fisher<-param.trans(theta.em) 477 478 Icom<-Icom.mvn(mu=mu, Sigma=Sigma, fix.rho=fix.rho, suff.stat=suff.stat, n=n) 479 Dvec<-Dcom.mvn(mu=mu, Sigma=Sigma, fix.rho=fix.rho, suff.stat=suff.stat, n=n) 480 481 482 theta.icom<-theta.em 483 if (fix.rho) theta.icom<-c(theta.em[-n.var], r12) 484 485 Icom.fisher<-Icom.transform(Icom=Icom, Dvec=Dvec, theta=theta.icom, transformation="Fisher", context=context, fix.rho=fix.rho) 486 487 488 Vcom.fisher <- solve(Icom.fisher) 489 490 if (!context) { 491 dV <- Vcom.fisher%*%DM%*%solve(diag(1,n.par)-DM) 492 Vobs.fisher <- Vcom.fisher+dV } 493 494 ###verify with the parameters. 495 ###repartition Icom 496 if (context & !fix.rho) { 497 index<-c(1,4,2,3,5,6,7,8,9) 498 Itemp<-Icom.fisher[index,index] 499 invItemp<-solve(Itemp) 500 A1<-invItemp[1:2,1:2] 501 A2<-invItemp[1:2,3:9] 502 A3<-invItemp[3:9, 1:2] 503 A4<-invItemp[3:9, 3:9] 504 dV1<-(A4-t(A2)%*%solve(A1)%*%A2)%*%DM%*%solve(diag(rep(1,7))-DM) 505 dV<-matrix(0,9,9) 506 dV[3:9,3:9]<-dV1 507 Vobs.fisher<-invItemp+dV 508 509 index2<-c(1,3,4,2,5,6,7,8,9) 510 Vobs.fisher<-Vobs.fisher[index2,index2] 511 } 512 513 514 515 Iobs.fisher <- solve(Vobs.fisher) 516 517 518 ##transform Iobs.fisher to Iobs via delta method 519 ##V(theta)=d(fisher^(-1))V(bvn.trans(theta))d(fisher^(-1))' 520 521 if (!context) { 522 grad.invfisher <- c(1,1, exp(theta.fisher[3:4])) 523 if (! fix.rho) 524 grad.invfisher <- c(grad.invfisher,4*exp(2*theta.fisher[5])/(exp(2*theta.fisher[5])+1)^2) 525 } 526 527 if (context) { 528 grad.invfisher <- c(1,1, 1, exp(theta.fisher[4:6])) 529 grad.invfisher <- c(grad.invfisher,4*exp(2*theta.fisher[7:8])/(exp(2*theta.fisher[7:8])+1)^2) 530 if (!fix.rho) 531 grad.invfisher <- c(grad.invfisher,4*exp(2*theta.fisher[9])/(exp(2*theta.fisher[9])+1)^2) 532 } 533 534 535 Vobs<-diag(grad.invfisher)%*%Vobs.fisher%*%diag(grad.invfisher) 536 Iobs<-solve(Vobs) 537 ## obtain a symmetric Cov matrix 538 Vobs.sym <- 0.5*(Vobs+t(Vobs)) 539 540###unitscale transformation 541 542#theta.unit<-param.trans(theta.em, transformation="unitscale") 543#Icom.unit<-Icom.transform(Icom, Dvec,theta.em, transformation="unitscale") 544#Vobs.unit<-delta method 545 546if (!context) { 547 names(mu)<-c("W1","W2") 548 colnames(Sigma)<-rownames(Sigma)<-c("W1","W2") 549 names(suff.stat)<-c("S1","S2","S11","S22","S12") 550 if (!fix.rho) colnames(DM)<-rownames(DM)<-c("u1","u2","s1","s2","r12") 551 if (fix.rho) colnames(DM)<-rownames(DM)<-c("u1","u2","s1","s2") 552} 553if (context) { 554 names(mu)<-c("X","W1","W2") 555 colnames(Sigma)<-rownames(Sigma)<-c("X","W1","W2") 556 names(suff.stat)<-c("Sx","S1","S2","Sxx","S11","S22","Sx1","Sx2","S12") 557 if (!fix.rho) { 558 colnames(DM)<-rownames(DM)<-c("u1","u2","s1","s2","r1x","r2x","r12") 559 colnames(Icom)<-rownames(Icom)<-c("ux","u1","u2","sx","s1","s2","r1x","r2x","r12") } 560 if (fix.rho) { 561 colnames(DM)<-rownames(DM)<-c("u1","u2","s1","s2","r1x","r2x") 562colnames(Icom)<-rownames(Icom)<-c("ux","u1","u2","sx","s1","s2","r1x","r2x") } 563} 564 565colnames(Iobs)<-colnames(Iobs.fisher)<-colnames(Icom.fisher)<-colnames(Vobs)<-colnames(Vobs.sym)<-colnames(Icom) 566rownames(Iobs)<-rownames(Iobs.fisher)<-rownames(Icom.fisher)<-rownames(Vobs)<-rownames(Vobs.sym)<-rownames(Icom) 567 568 res.out<-list(mu=mu, Sigma=Sigma, suff.stat=suff.stat, context=context, fix.rho=fix.rho) 569 res.out$DM<-DM 570 res.out$Icom<-Icom 571 res.out$Iobs<-Iobs 572 res.out$Fmis<-1-diag(Iobs)/diag(Icom) 573 res.out$Vcom<-Vcom<-solve(Icom) 574 res.out$Vobs.original<-Vobs 575 res.out$VFmis<-1-diag(Vcom)/diag(Vobs) 576 res.out$Vobs<-Vobs.sym 577 res.out$Icom.trans<-Icom.fisher 578 res.out$Iobs.trans<-Iobs.fisher 579 res.out$Fmis.trans<-1-diag(Iobs.fisher)/diag(Icom.fisher) 580 res.out$Imiss<-res.out$Icom-res.out$Iobs 581 res.out$Ieigen<-eigen(res.out$Imiss)[[1]][1] 582res.out 583} 584