[R] EM Algorithm for missing data

Sudhi.Upadhyaya sudhi.upadhyaya at ndsu.edu
Wed Dec 7 07:44:15 CET 2011


Dear all,

I need help with EM algorithm. 

I am modeling this alogirthm based on " Incomplete Data in Generalized
Linear Models" by Joseph G Ibrahim. 

i have half way through developing the R code based on this this paper, I
need little help in tweaking my code furthure.

----------------------------------------------------------------------------
a=read.table("C:/sudhi/ex3.csv",header=T,sep=",")

#Set design matrix
n=1558
X1=as.matrix(rep(1,n))
X2=as.matrix((a[,1]))
X3=as.matrix((a[,2]))
X4=as.matrix((a[,3]))
X5=as.matrix((a[,4]))
X6=as.matrix((a[,5]))
X7=as.matrix((a[,6]))
X8=as.matrix((a[,7]))
X=cbind(X1,X2,X3,X4,X5,X6,X7,X8)


#Start Newton-Raphson Method
#Set Initial parameter
B0=0.001
B1=0.001
B2=0.001
B3=0.001
B4=0.001
B5=0.001
B6=0.001
B7=0.001

#Parameters
B=matrix(c(B0,B1,B2,B3,B4,B4,B6,B7),nrow=8,ncol=1,byrow=F)

#Start Iteration
Diff=1
AA = 0
while(Diff>.000001){

AA= AA +1
plot_colors <- c("blue","red","forestgreen")
plot(B, type="l", lwd=2, col=plot_colors[1], add=TRUE) 
split.screen(c(1,1))
#Probability for each comb
P=(1+exp(-(X%*%B)))^-1

#Set V matrix
ni=1
V=diag(n)
for(i in 1:n)
{
V[i,i]=diag(ni%*%P[i]%*%(1-P[i]))
}

#Set S matrix
#Response
R=a[,8]
S=matrix(data=1,nrow=n,ncol=1)
for (i in 1:n)
{S[i]=R[i]-ni*P[i]}

new.B=B+solve(t(X)%*%V%*%X)%*%(t(X)%*%S)
Diff=max(abs(new.B-B))


B=new.B}


#Check the second derivative
-diag(t(X)%*%V%*%X)

#Check the first derivative
t(X)%*%S

#Maximum Likelihood Estimates

---------------------------------------------------

Please see the pdf file attached for Joseph G Ibrahim's logic
Sincerely
Sudhi
http://r.789695.n4.nabble.com/file/n4167861/EM_on_Logit_Model_-Ibrahim.PDF
EM_on_Logit_Model_-Ibrahim.PDF 

--
View this message in context: http://r.789695.n4.nabble.com/EM-Algorithm-for-missing-data-tp4167861p4167861.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list