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