#General dependence initialcov<-function(X,meanvec) { for (i in 1:ncol(X)) { X[is.na(X[,i]),i]=meanvec[i] } covmat=cov(X)*(nrow(X)-1)/nrow(X) return(covmat) } empredict<-function(X, emmean, emcov) { completeX=X rowmiss=sort(unique(row(X)[is.na(X)])) nrowmiss=length(rowmiss) #number of rows with missing T2=0 for (i in 1:nrowmiss) { wholerow=X[rowmiss[i],] x2=wholerow[!is.na(wholerow)] colmiss=col(t(as.matrix(wholerow)))[is.na(t(as.matrix(wholerow)))] mu1=emmean[colmiss] #missing mu2=emmean[-colmiss] #available sigma11=emcov[colmiss,colmiss] sigma12=emcov[colmiss,-colmiss] sigma21=emcov[-colmiss,colmiss] sigma22=emcov[-colmiss,-colmiss] x1=mu1+sigma12%*%solve(sigma22)%*%(x2-mu2) wholerow[is.na(wholerow)]=x1 completeX[rowmiss[i],]=wholerow matprod=completeX[rowmiss[i],]%*%t(completeX[rowmiss[i],]) x1x1=sigma11-sigma12%*%solve(sigma22)%*%sigma21+x1%*%t(x1) matprod[colmiss,colmiss]=x1x1 T2=T2+matprod } T1=apply(completeX,2,sum) rownomiss=seq(1,nrow(X))[-rowmiss] nrownomiss=length(rownomiss) for (i in 1:nrownomiss) { wholerow=X[rownomiss[i],] T2=T2+wholerow%*%t(wholerow) } return(list(T1=T1, T2=T2)) } emestimate<-function(T1, T2, n) { means=T1/n covs=T2/n-means%*%t(means) return(list(means=means, covs=covs)) } EM<-function(X, delta=0.0001) #X is a n*p matrix with missing { X=as.matrix(X) #initial estimate emmean=apply(X,2,mean,na.rm=TRUE) emcov=initialcov(X, emmean) if (min(eigen(emcov)$values)<0) { stop("the initial estimate of covariance matrix is negative definite") } n=nrow(X) #number of rows meandiff=9999999 covdiff=9999999 while(meandiff>delta & covdiff>delta) { #prediction T1=empredict(X, emmean, emcov)$T1 T2=empredict(X, emmean, emcov)$T2 #estimate tempmean=emestimate(T1, T2, n)$means tempcov=emestimate(T1, T2, n)$covs meandiff=sqrt(sum((tempmean-emmean)^2)) covdiff=sqrt(sum((tempcov-emcov)^2)) emmean=tempmean emcov=tempcov } vars=diag(emcov) emcorr=emcov/(sqrt(vars%*%t(vars))) return(list(emmean=emmean, emcov=emcov, emcorr=emcorr)) } mortality=read.table("E:/Study/RA/2009/20090210/mortalitydata.txt") mortality=as.matrix(mortality) mortmat=matrix(as.numeric(mortality[,3:13]),nrow=112,ncol=11,byrow=F) emout=EM(mortmat,delta=0.0000001) > emout $emmean [1] -3.8569634 -3.5538932 -3.5791981 -2.8521299 -2.0904094 -1.2307405 [7] -0.7073208 -0.3072631 -0.4205684 -0.3617684 -0.1566562 $emcov [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.68143978 0.5106313 0.580665459 0.62811454 0.4873835 0.3408798 [2,] 0.51063125 1.0896582 0.954125488 1.05307029 0.9274666 0.6462927 [3,] 0.58066546 0.9541255 1.672149763 1.66067714 1.3151113 0.8551698 [4,] 0.62811454 1.0530703 1.660677142 2.72075293 2.1762361 1.4492108 [5,] 0.48738346 0.9274666 1.315111275 2.17623612 2.6043326 1.9213868 [6,] 0.34087981 0.6462927 0.855169768 1.44921076 1.9213868 1.9210149 [7,] 0.19041837 0.3992355 0.438061048 0.70453013 1.1077855 1.1993219 [8,] 0.07970121 0.2489312 0.291751910 0.33725877 0.5972812 0.6647720 [9,] 0.01154820 0.1125049 0.013661102 0.03347331 0.1738759 0.1776150 [10,] 0.05180805 0.1736873 -0.005410955 0.12691112 0.1883594 0.1705523 [11,] 0.13418132 0.2151106 0.077037932 0.41686552 0.4890497 0.3880390 [,7] [,8] [,9] [,10] [,11] [1,] 0.1904184 0.07970121 0.01154820 0.051808050 0.13418132 [2,] 0.3992355 0.24893125 0.11250493 0.173687267 0.21511061 [3,] 0.4380610 0.29175191 0.01366110 -0.005410955 0.07703793 [4,] 0.7045301 0.33725877 0.03347331 0.126911122 0.41686552 [5,] 1.1077855 0.59728116 0.17387586 0.188359431 0.48904971 [6,] 1.1993219 0.66477197 0.17761504 0.170552263 0.38803901 [7,] 1.1302852 0.67118374 0.31611448 0.235565969 0.26095642 [8,] 0.6711837 0.67159869 0.36033514 0.209899701 0.19769466 [9,] 0.3161145 0.36033514 0.47303160 0.213011813 0.14005889 [10,] 0.2355660 0.20989970 0.21301181 0.357099935 0.17746338 [11,] 0.2609564 0.19769466 0.14005889 0.177463384 0.47939606 $emcorr [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.00000000 0.5925821 0.543969645 0.46129692 0.3658549 0.2979354 [2,] 0.59258206 1.0000000 0.706842668 0.61160025 0.5505603 0.4467029 [3,] 0.54396964 0.7068427 1.000000000 0.77857969 0.6301975 0.4771439 [4,] 0.46129692 0.6116002 0.778579688 1.00000000 0.8175480 0.6339013 [5,] 0.36585485 0.5505603 0.630197524 0.81754796 1.0000000 0.8590159 [6,] 0.29793545 0.4467029 0.477143912 0.63390128 0.8590159 1.0000000 [7,] 0.21697069 0.3597411 0.318642007 0.40175473 0.6456742 0.8139100 [8,] 0.11781382 0.2909909 0.275309502 0.24949611 0.4516226 0.5852648 [9,] 0.02034022 0.1567045 0.015360419 0.02950591 0.1566556 0.1863242 [10,] 0.10502402 0.2784376 -0.007002312 0.12875386 0.1953189 0.2059194 [11,] 0.23476377 0.2976250 0.086043843 0.36500952 0.4376809 0.4043550 [,7] [,8] [,9] [,10] [,11] [1,] 0.2169707 0.1178138 0.02034022 0.105024016 0.23476377 [2,] 0.3597411 0.2909909 0.15670445 0.278437639 0.29762502 [3,] 0.3186420 0.2753095 0.01536042 -0.007002312 0.08604384 [4,] 0.4017547 0.2494961 0.02950591 0.128753856 0.36500952 [5,] 0.6456742 0.4516226 0.15665562 0.195318896 0.43768085 [6,] 0.8139100 0.5852648 0.18632421 0.205919386 0.40435498 [7,] 1.0000000 0.7703578 0.43231996 0.370786268 0.35450856 [8,] 0.7703578 1.0000000 0.63930281 0.428609718 0.34841174 [9,] 0.4323200 0.6393028 1.00000000 0.518279138 0.29411580 [10,] 0.3707863 0.4286097 0.51827914 1.000000000 0.42891047 [11,] 0.3545086 0.3484117 0.29411580 0.428910466 1.00000000 Order1 antedependence atdp<-function(covmat) { n=nrow(covmat) newmat=matrix(0,nrow=n,ncol=n) diagonal=diag(covmat) subdiag1=rep(0,n-1) for (i in 1:(n-1)) { subdiag1[i]=covmat[i+1,i]/sqrt(diagonal[i]*diagonal[i+1]) } diag(newmat)=diagonal for (i in 2:n) { for (j in 1:(i-1)) { sindex=min(i,j) lindex=max(i,j) newmat[i,j]=sqrt(diagonal[i]*diagonal[j])*prod(subdiag1[sindex:(lindex-1)]) newmat[j,i]=newmat[i,j] } } return(newmat=newmat) } initialcov<-function(X, meanvec) { for (i in 1:ncol(X)) { X[is.na(X[,i]),i]=meanvec[i] } covmat=cov(X)*(nrow(X)-1)/nrow(X) covmat.atdp=atdp(covmat) return(covmat.atdp) } empredict<-function(X, emmean, emcov) { completeX=X rowmiss=sort(unique(row(X)[is.na(X)])) nrowmiss=length(rowmiss) #number of rows with missing T2=0 for (i in 1:nrowmiss) { wholerow=X[rowmiss[i],] x2=wholerow[!is.na(wholerow)] colmiss=col(t(as.matrix(wholerow)))[is.na(t(as.matrix(wholerow)))] mu1=emmean[colmiss] #missing mu2=emmean[-colmiss] #available sigma11=emcov[colmiss,colmiss] sigma12=emcov[colmiss,-colmiss] sigma21=emcov[-colmiss,colmiss] sigma22=emcov[-colmiss,-colmiss] x1=mu1+sigma12%*%solve(sigma22)%*%(x2-mu2) wholerow[is.na(wholerow)]=x1 completeX[rowmiss[i],]=wholerow matprod=completeX[rowmiss[i],]%*%t(completeX[rowmiss[i],]) x1x1=sigma11-sigma12%*%solve(sigma22)%*%sigma21+x1%*%t(x1) matprod[colmiss,colmiss]=x1x1 T2=T2+matprod } T1=apply(completeX,2,sum) rownomiss=seq(1,nrow(X))[-rowmiss] nrownomiss=length(rownomiss) for (i in 1:nrownomiss) { wholerow=X[rownomiss[i],] T2=T2+wholerow%*%t(wholerow) } return(list(T1=T1, T2=T2)) } emestimate<-function(T1, T2, n) { means=T1/n covs=T2/n-means%*%t(means) covs.atdp=atdp(covs) return(list(means=means, covs.atdp=covs.atdp)) } ATDPEM<-function(X, delta=0.0001) #X is a matrix with missing { X=as.matrix(X) #initial estimate emmean=apply(X,2,mean,na.rm=TRUE) emcov=initialcov(X, emmean) if (min(eigen(emcov)$values)<0) { stop("the initial estimate of covariance matrix is negative definite") } n=nrow(X) meandiff=9999999 covdiff=9999999 while(meandiff>delta & covdiff>delta) { #prediction T1=empredict(X, emmean, emcov)$T1 T2=empredict(X, emmean, emcov)$T2 #estimate tempmean=emestimate(T1, T2, n)$means tempcov=emestimate(T1, T2, n)$covs.atdp meandiff=sqrt(sum((tempmean-emmean)^2)) covdiff=sqrt(sum((tempcov-emcov)^2)) emmean=tempmean emcov=tempcov } vars=diag(emcov) emcorr=emcov/(sqrt(vars%*%t(vars))) return(list(emmean.ad1=emmean, emcov.ad1=emcov, emcorr.ad1=emcorr)) } emad1out=ATDPEM(mortmat,delta=0.0000001) > emad1out $emmean.ad1 [1] -3.8655356 -3.5459324 -3.5898244 -2.8568650 -2.0806447 -1.2163396 [7] -0.6772133 -0.2680650 -0.3612898 -0.4002523 -0.3451259 $emcov.ad1 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.68538788 0.51714140 0.46542510 0.46781019 0.36919819 0.27533996 [2,] 0.51714140 1.06544100 0.95889246 0.96380634 0.76064087 0.56726938 [3,] 0.46542510 0.95889246 1.69300398 1.70167985 1.34297441 1.00156104 [4,] 0.46781019 0.96380634 1.70167985 2.77770687 2.19218041 1.63488037 [5,] 0.36919819 0.76064087 1.34297441 2.19218041 2.56568578 1.91343253 [6,] 0.27533996 0.56726938 1.00156104 1.63488037 1.91343253 1.93272432 [7,] 0.17813619 0.36700524 0.64797812 1.05771558 1.23792996 1.25041113 [8,] 0.10823872 0.22299892 0.39372305 0.64268682 0.75218829 0.75977207 [9,] 0.05877436 0.12108993 0.21379430 0.34898332 0.40844337 0.41256141 [10,] 0.02479584 0.05108565 0.09019594 0.14722973 0.17231485 0.17405218 [11,] 0.01008339 0.02077430 0.03667875 0.05987191 0.07007293 0.07077943 [,7] [,8] [,9] [,10] [,11] [1,] 0.17813619 0.1082387 0.05877436 0.02479584 0.01008339 [2,] 0.36700524 0.2229989 0.12108993 0.05108565 0.02077430 [3,] 0.64797812 0.3937230 0.21379430 0.09019594 0.03667875 [4,] 1.05771558 0.6426868 0.34898332 0.14722973 0.05987191 [5,] 1.23792996 0.7521883 0.40844337 0.17231485 0.07007293 [6,] 1.25041113 0.7597721 0.41256141 0.17405218 0.07077943 [7,] 1.18975914 0.7229188 0.39254985 0.16560966 0.06734623 [8,] 0.72291884 0.7108371 0.38598936 0.16284191 0.06622070 [9,] 0.39254985 0.3859894 0.48429028 0.20431329 0.08308530 [10,] 0.16560966 0.1628419 0.20431329 0.34558732 0.14053529 [11,] 0.06734623 0.0662207 0.08308530 0.14053529 0.41070859 $emcorr.ad1 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.00000000 0.60516861 0.43206849 0.33904571 0.27841307 0.23923024 [2,] 0.60516861 1.00000000 0.71396383 0.56025001 0.46005867 0.39531172 [3,] 0.43206849 0.71396383 1.00000000 0.78470363 0.64437252 0.55368592 [4,] 0.33904571 0.56025001 0.78470363 1.00000000 0.82116674 0.70559877 [5,] 0.27841307 0.46005867 0.64437252 0.82116674 1.00000000 0.85926370 [6,] 0.23923024 0.39531172 0.55368592 0.70559877 0.85926370 1.00000000 [7,] 0.19726699 0.32597030 0.45656417 0.58183008 0.70854072 0.82459055 [8,] 0.15507058 0.25624360 0.35890278 0.45737367 0.55698026 0.64820644 [9,] 0.10201566 0.16857395 0.23610993 0.30089058 0.36641837 0.42643297 [10,] 0.05094853 0.08418899 0.11791772 0.15027039 0.18299620 0.21296861 [11,] 0.01900515 0.03140472 0.04398643 0.05605483 0.06826242 0.07944292 [,7] [,8] [,9] [,10] [,11] [1,] 0.19726699 0.1550706 0.1020157 0.05094853 0.01900515 [2,] 0.32597030 0.2562436 0.1685739 0.08418899 0.03140472 [3,] 0.45656417 0.3589028 0.2361099 0.11791772 0.04398643 [4,] 0.58183008 0.4573737 0.3008906 0.15027039 0.05605483 [5,] 0.70854072 0.5569803 0.3664184 0.18299620 0.06826242 [6,] 0.82459055 0.6482064 0.4264330 0.21296861 0.07944292 [7,] 1.00000000 0.7860949 0.5171451 0.25827195 0.09634227 [8,] 0.78609492 1.0000000 0.6578660 0.32855060 0.12255806 [9,] 0.51714511 0.6578660 1.0000000 0.49941873 0.18629639 [10,] 0.25827195 0.3285506 0.4994187 1.00000000 0.37302643 [11,] 0.09634227 0.1225581 0.1862964 0.37302643 1.00000000 > emad1out$emcorr.ad1-emout$emcorr [,1] [,2] [,3] [,4] [,5] [1,] 0.00000000 0.01258655 -0.111901150 -0.122251203 -0.0874417869 [2,] 0.01258655 0.00000000 0.007121160 -0.051350244 -0.0905016560 [3,] -0.11190115 0.00712116 0.000000000 0.006123939 0.0141749986 [4,] -0.12225120 -0.05135024 0.006123939 0.000000000 0.0036187843 [5,] -0.08744179 -0.09050166 0.014174999 0.003618784 0.0000000000 [6,] -0.05870521 -0.05139114 0.076542005 0.071697488 0.0002477733 [7,] -0.01970370 -0.03377075 0.137922166 0.180075351 0.0628665417 [8,] 0.03725676 -0.03474726 0.083593274 0.207877556 0.1053576589 [9,] 0.08167544 0.01186950 0.220749508 0.271384669 0.2097627422 [10,] -0.05407548 -0.19424865 0.124920032 0.021516535 -0.0123226997 [11,] -0.21575862 -0.26622031 -0.042057416 -0.308954695 -0.3694184349 [,6] [,7] [,8] [,9] [,10] [1,] -0.0587052068 -0.01970370 0.03725676 0.08167544 -0.054075484 [2,] -0.0513911351 -0.03377075 -0.03474726 0.01186950 -0.194248651 [3,] 0.0765420046 0.13792217 0.08359327 0.22074951 0.124920032 [4,] 0.0716974884 0.18007535 0.20787756 0.27138467 0.021516535 [5,] 0.0002477733 0.06286654 0.10535766 0.20976274 -0.012322700 [6,] 0.0000000000 0.01068054 0.06294163 0.24010875 0.007049226 [7,] 0.0106805417 0.00000000 0.01573712 0.08482514 -0.112514314 [8,] 0.0629416288 0.01573712 0.00000000 0.01856317 -0.100059123 [9,] 0.2401087539 0.08482514 0.01856317 0.00000000 -0.018860404 [10,] 0.0070492260 -0.11251431 -0.10005912 -0.01886040 0.000000000 [11,] -0.3249120561 -0.25816630 -0.22585368 -0.10781941 -0.055884034 [,11] [1,] -0.21575862 [2,] -0.26622031 [3,] -0.04205742 [4,] -0.30895469 [5,] -0.36941843 [6,] -0.32491206 [7,] -0.25816630 [8,] -0.22585368 [9,] -0.10781941 [10,] -0.05588403 [11,] 0.00000000