Note on the EM Algorithm in Linear Regression …

1 International Mathematical Forum, 4, 2009 , no. 38, 1883 - 1889 Note on the EM Algorithm in LinearRegression ModelJi-Xia Wang and Yu MiaoCollege of Mathematics and Information ScienceHenan Normal UniversityHenan Province, 453007, Regression model has been used extensively in the fields ofinformation processing and data analysis. In the present paper, we con-sider the Linear model with missing data. Using the EM (Expectationand Maximization) Algorithm , the asymptotic variances and the stan-dard errors for the MLE of the unknown parameters are Subject Classification:93C05; 93C41 Keywords:Conditional expectation; maximum likelihood estimator; EMalgorithm; Newton-Raphson iteration1 IntroductionAs a typical statistical model , Linear Regression model has been widely used inthe fields of information processing and data analysis. In fact, there have beenseveral statistical methods for its learning or modeling ( , the expectation-maximization (EM) Algorithm [2] for maximum likelihood and the self-organizingnetwork with hyper-ellipsoidal clustering [5]).

2 Generally, the parameters of lin-ear regressive model can be estimated via the EM Algorithm under the maxi-mum likelihood framework, since the EM Algorithm owns certain good conver-gence behaviors in certain situations. However, in some applications, there aremany data sets including missing observations [9], which cause many problemsif the missing data is related to the values of the missing item [8], for instance,in [4], Little and Rubin showed that this can cause bias and inefficiency forsome estimations. So, an new Algorithm for estimating unknown parametersis proposed based on the likelihood function. In [1], Baker and Laird used Wang and Y. MiaoEM Algorithm to obtain maximum likelihood estimates (MLE) of the unknownparameters in the model with the incomplete data. Ibrahim and Lipsitz [3]established Bayesian methods for estimation in generalized Linear the present paper, we discuss the Linear Regression model with miss-ing data and propose a method for estimating parameters by using Newton-Raphson iteration to solve the score equation.

3 Moreover, the standard errorsof these estimators are calculated by the observed Fisher information Linear Regression model with missing dataSuppose thaty1,y2,.. ,ynare independent identically distributed normal ran-dom variables with unit variances. LetXi=(X1i,X2i)Tis a 2 1 random vec-tor of covariation, whereX1iandX2iare independent observations and follownormal distributions with means 1, 2and variances 21, 22, respectively. Fornotation convenience, let Xi=(1,X1i,X2i)Tand assume that T=( 0, 1, 2)are Regression coefficients. It is also supposed thatp(yi|Xi, )=1 2 exp yi XiT 22 .(1)We assume thatX1iis completely observed, andX2iis partially missing foreveryiand our objective is to estimate , 1, 2, 21, 22and their standarderrors from the known data with missing value indicators are introduced in [6] asri= 0,ifyiis observed,1,ifyiis 0,ifx2iis observed,1,ifx2iis missing.(2)with probabilitiesp(ri)= i,p(si)= i. Following the reference [8], for anyi=1,2.

4 ,n, the missing-data mechanism is defined aslogit( i) log i1 i= 1X1i+ 2X2i+yi (3)andlogit( i) log i1 i= 1X1i+ 2X2i+yi ,(4)where =( 1, 2)T, =( 1, 2)T, and are parameters determining themissing mechanism. Then the conditional probability functions forriandsiare derived by Eqs. (2)-(4) asp(ri|Xi,yi, , )=exp{ri(XTi +yi )}1 + exp{XTi +yi },Note on the EM algorithm1885p(si|Xi,yi, , )=exp{si(XTi +yi )}1 + exp{XTi +yi }.Now we derive the joint probability function ofyi,x2i,ri,siasp(yi,x2i,ri,si|x1i)=p(ri |Xi,yi, , )p(si|Xi,yi, , )p(yi|Xi, )p(x2i|X1i) exp{ri(XTi +yi )}1 + exp{XTi +yi } exp{si(XTi +yi )}1 + exp{XTi +yi } (2 ) 12 exp (yi XiT )22 (2 22) 12 exp (xTherefore, we can write down the complete-data log-likelihoodl( )bylogL( |yi,Xi,ri,si)=n i=1log exp{ri(XTi +yi )}1 + exp{ri(XTi +yi )} +n i=1log exp{ri(XTi +yi )}1 + exp{si(XTi +yi )} +n2log(2 ) n i=1 yi XiT 22 n2log(2 22) n i=1(x2i 2)22 22,where =( , , , , , 2, 22) is the parameter related to developing EM al-gorithm.)

5 The complete-data log-likelihood specifies a model for the joint char-acterization of the observed data and the associated missing-data E-step of EM algorithmThe MLE of is a point which maximizes the observed-data likelihood functionL( |(y, X)obs,ri,si), where (y, X)obsis the observed components of (y, X). Let (r)be ther-st iteration estimate of and define the conditional expectation ofl( )-with respect to the conditional distribution of the missing data (y, X)misgiven the observed datayi,Xi,ri,siand the value (r)as the following:Q( | (r))=E[l( )|(y, X)obs,r,s, (r)].(5)The EM Algorithm is composed of E-step and M-step iterations. Nowfor the expectation of the complete-data log-likelihood in the E-step of EMalgorithm, we consider four possible-cases: response variableyiis missing, acovariancex2iis missing, both of them are missing, and no missing the expected log-likelihood function can be written by(6)= Wang and Y. Miaowherex2i,misdenotes the missing components ofx2i.

6 Eqs.( ) and ( ) leadto the conditional expectation ofl( ), which is our target quantity asQ( | (r))=n1 i=1l( )+n2 i=n1+1 l( )p yi,mis|Xi,ri,si, (r) dyi,mis+n3 i=n2+1 l( )p x2i,mis|Xi,obs,yi,ri,si, (r) dx2i,mis+n i=n3+1 yi=1 l( )p yi,mis,x2i,mis|Xi,obs,ri,si, (r) dyi,misdx2i,miswheren1,n2,n3are corresponding sample sizes,yi,misis the missing compo-nents ofyi,Xi,obsis the observed component ofXi, andp(yi,mis,x2i,mis|Xi,obs,ri,si),p(yi,m is|Xi,ri,si) andp(yi,mis,x2i,mis|Xi,obs,ri,si) are the conditional probabil-ities of the missing data given the observed data. These conditional probabil-ities are regarded as the weights inQ( | (r)). The weights have the followingform:p yi,mis,x2i,mis|Xi,obs,ri,si, (r) =p yi|Xi, (r) p(x2i|x1i)p ri|yi,Xi, (r) p si|yi,Xi, (r) y1=1 p(yi|Xi, (r))p(x2i|x1i)p(ri|yi,Xi, (r))p(si|yi,Xi, (r)) p yi,x2i,ri,si|x1i, (r) ,p x2i,mis|Xi,obs,yi,ri,si, (r) =p x2i|x1i, (r) p si|yi,Xi, (r) p(x2i|x1i, (r))p(si|yi,Xi, (r)) exp{ri(XTi +yi )}1 + exp{XTi +yi } (2 22) 12 exp (x2i 2)22 22 ,andp yi,mis|Xi,ri,si, (r) =p yi|Xi, (r) p ri|yi,Xi, (r) yi=1p(yi|Xi, (r))p(ri|yi,Xi, (r)) p yi|Xi, (r) p ri|yi,Xi, (r).

7 Then the conditional expectationQ( | (r)) is to be calculated by a Metropolis-Hastings(MH) Algorithm [7].Note on the EM algorithm18874 M-step of EM Algorithm and convergenceNow we need to find a value of , saying (r), at whichQ( | (r)) will attainthe maximum. The Newton-Raphson method will be used to solve the scoreequation. The parameters (r+1)in the M-step at the (r+1)stEM iterationand the (r+1)stNewton-Raphson iteration take the following form (for forexample): (r+1)= (r)+ 2Q( | (r)) T 1 = (r) Q( | (r)) = (r).The derivatives of the parameter used in the iteration are given as follows: Q( | (r)) =n1 i=1 Xi yi XiT +n2 i=n1+1E Xi yi XiT |Xi, (r) +n3 i=n2+1E Xi yi XiT |Xobs,yi, (r) +n i=n3+1E Xi yi XiT |Xobs, (r) ,and 2Q( | (r)) T=n1 i=1 XiT Xi+n2 i=n1+1E XiT Xi|Xi, (r) +n3 i=n2+1E XiT Xi|Xobs,yi, (r) +n i=n3+1E XiT Xi|Xobs, (r) .The derivatives of other components of used in the iteration are given in thereference [6].

8 The (r+1)stestimates of 2, 22are obtained by solving the score equations: Q( | (r)) 2=n i=1E(x2i|x1i,yi,ri,si) n 2=0, Q( | (r)) 22=n i=1E (x2i 2)2|x1i,yi,ri,si n 22= , we can take (r+1)2, 2(r+1)2by (r+1)2=1nE(x2i|x1i,yi,ri,si), 2(r+1)2=1nE (x2i 2)2|x1i,yi,ri,si , Wang and Y. Miaowhich are approximated by the sample averages of simulated and given sequence{Q( | (r))}often exhibits an increasing trend, and then fluc-tuate around the value ofQ( | (r))ifrbecomes large enough. The sequence{ (r)}would also fluctuate the MLE (r)whenris sufficiently large. To monitorthe convergence of the EM Algorithm we can plot{Q( | (r))}as well as{ (r)}against iteration number. We terminate the Algorithm when the sequenceof{Q( | (r))}become stationary. Otherwise, we continue by increasing theMonte Carlo precision in the E-step provided calculation is Standard errors of estimatesIt is well know that the distribution of maximum likelihood estimates asymp-totically tends to a normal distributionMV N( , V( )) under some regularityconditions.

9 The expected Fisher information matrixI( ) which gives the in-verse of variance matrix of is approximated by the observed informationmatrixJ (Y):V( ) 1=nE 2logL( ) 2 = n 2logL( ) 2 dx n i=1 2logL( ) 2 = nJ( ).By using the following relation which is obtained in [9]:observed informa-tion=complete information-missing information,we haveI( ) J (Y)= 2logL( ) 2= 2Q( | (r)) 2 Var n i=1 logL( ) = ,whereVar( ) is the conditional variance given (y, X)obs,r,s, and (r). Thedetails are to be provided in the reference [6]. authors acknowledge the financial support of the Foundation for Dis-tinguished Young Scholars of Henan Province (084100510013).References[1] S. G. Baker and N. M. Laird, Regression analysis for categorical variableswith outcome subject to nonignorable nonresponse, J. Am. Stat. Assoc,1988,83 on the EM algorithm1889[2] A. P. Dempster, N. M. Laird and D. B. Rubin,Maximum likelihood fromincomplete data via the EM Stat.

10 Soc. B, 1977,39:1-38.[3] J. G. Ibrahim, S. R. Lipsitz,Missing covariates in generalized Linear mod-els when the missing data mechanism is non-ignorable, J. Royal Stat. , 1999,61: 173-190.[4] R. J. A. Little and D. B. Rubin,Statistical Analysis with Missing Data,New York, Wiley, 2002.[5] J. Mao and A. K. Jain,A self-organizing network for hyperellipsoidalclustering, IEEE Trans. Neural Networks, 1996,7(1): 16-29.[6] J. S. Park, G. Q. Qian and Y. Jun,Monte Carlo EM Algorithm in logisticlinear models involving non-ignorable missing data, Appl. Math. Comput.,2008,197: 440-450.[7] C. P. Robert and G. Casella,Monte Carlo Statistical Methods, Berlin:Springer, 1999.[8] M. M. Rueda,S. Gonzalez and A. Arcos,Indirect methods of imputation ofmissing data based on available units, Appl. Math. Comput., 2005,164:249-261.[9] Y. G. Smirlis and E. K. Despotis,Data envelopment analysis with missingvalues: An interval DEA approach, Appl.