#Convert standard format to the table format
logLinear <- function(marker,n,pedInfo,missing="Direct"){
# find child and his/her parents
child <- (pedInfo$patid>0 | pedInfo$matid>0)
patid <- pedInfo$patid[child]
matid <- pedInfo$matid[child]
chtid <- pedInfo$id[child]
subPedid <- pedInfo$pedid[child]

# sort the data, and write the sorted data into a new dataset called newdata
data <- data.frame(cbind(pedInfo,marker))
newdata <- data.frame(matrix(0,nrow=n*3,ncol=8))
names(newdata) <- names(data)

for (i in 1:n){
	temp <- data[(data$id==matid[i] & data$pedid==subPedid[i]),]
	if(dim(temp)[1] > 0) newdata[(3*(i-1)+1),] <- temp
	temp <- data[(data$id==patid[i] & data$pedid==subPedid[i]),]
	if(dim(temp)[1] > 0) newdata[(3*(i-1)+2),] <- temp
	temp <- data[(data$id==chtid[i] & data$pedid==subPedid[i]),]
	if(dim(temp)[1] > 0) newdata[(i*3),] <- temp
	}
pedInfo <- newdata

nCopiesmarker <- (pedInfo[,7]+pedInfo[,8])-2
nCopiesmarker[nCopiesmarker < 0] <- NA

# Use pedInfo to generate matrix of MFC assuming that data may be out of order
MFC <- matrix(nCopiesmarker,n,3,byrow=T)

# Error in genotyping code
charMFC <- paste(MFC[,1],MFC[,2],MFC[,3],sep="")
# Code that shows which pedigree has the error in genotyping
errorMFC <- c("000","010","011","021","0NA0","0NA1","100","101",
	"110","111","112","121","122","1NA0","1NA1","1NA2","201",
	"211","212","222","2NA1","2NA2","NA00","NA01","NA10","NA11",
	"NA12","NA21","NA22")
numError <- ifelse(is.element(charMFC,errorMFC),0,1)
if(sum(numError)>0){stop("Genotyping errors in the data")}

numberMFC <- MFC[,1]*100+MFC[,2]*10+MFC[,3]*1
countMFC <- c(numberMFC,c(222,212,211,201,122,121,112,111,110, 	101,100,021,011,010,000))
MFCdata <- summary(factor(countMFC))-1

# Dummy variables for Poisson model
nCopy1 <- c(0,0,1,1,0,1,0,1,0,1,0,1,1,0,0)
nCopy2 <- c(0,0,0,0,0,0,0,0,1,0,1,0,0,1,1)
hetOff <- c(0,0,0,0,0,0,0,1,0,0,0,0,0,0,0)
intStrat <- c(6,5,5,3,5,5,4,4,4,2,2,3,2,2,1)

#Poisson model fit
outNew <- glm(MFCdata[1:15]~C(factor(intStrat),treatment,base=6)+
	nCopy1+nCopy2,family=poisson,na.action=na.omit,
	offset=log(2)*hetOff)
outRed <- glm(MFCdata[1:15]~C(factor(intStrat),treatment,base=6),
	family=poisson,na.action=na.omit,offset=log(2)*hetOff)
beta <- as.numeric(c(outNew$coef[2:6]+outNew$coef[1],
	outNew$coef[1],outNew$coef[7:8]))
rej <- (1-pchisq(anova(outRed,outNew)$Deviance[2],2))

MFCfunctions <- MFCdata[c(1:3,5:9,4,12,10,11,13:15)]
paEst <- (sum(MFCfunctions[2:5])+2*sum(MFCfunctions[6:10])+3*sum(
	MFCfunctions[11:14])+4*MFCfunctions[15])/(4*sum(
	MFCfunctions[1:15]))

if(sum(is.na(marker))>0 | dim(marker)[1] < n*3){
#######  define observed incompelete families M(i,?,k) ########

Mdad=array(0,c(3,3))
for(i in 1:3){
  for (k in 1:3){
    Mdad[i,k]=sum((MFC[,1]==(i %% 3))&is.na(MFC[,2])&(MFC[,
    3]==(k %% 3)),na.rm=T)
}}
########define M(?,j,k)
Mmom=array(0,c(3,3))
for(j in 1:3){
  for (k in 1:3){
    Mmom[j,k]=sum(is.na(MFC[,1])&(MFC[,2]==(j %% 3))&(MFC[,
    3]==(k %% 3)),na.rm=T)
}}

M=Mmom+Mdad
###(sum(M)+sum(MFCdata[1:15]))==n
M22 = M[2,2]
M21 = M[2,1]
M12 = M[1,2]
M11 = M[1,1]
M01 = M[3,1]
M10 = M[1,3]
M00 = M[3,3]

#MLE estimate of allele frequency
paEst <- (sum(MFCfunctions[2:5])+2*sum(MFCfunctions[6:10])+3*sum(
	MFCfunctions[11:14])+4*MFCfunctions[15]+M01+M10+M11+2*(
	M12+M21)+3*M22)/(4*sum(MFCfunctions[1:15])+3*(
	M00+M01+M10+M12+M21+M22)+2*M11)
source("triadMultinomial.R")
# Complete data guess at multinomial under H_A
multi.opt <- optim(c(0,0,paEst),p_zip,dp_zip,MFCfunctions,M22,M21,
	M12,M11,M01,M10,M00,method="BFGS",control=list(fnscale=-1))
#Initializaiton 
u1= log(n)+4*log(multi.opt$par[3])
u2= log(n)+3*log(multi.opt$par[3])+log(1-multi.opt$par[3])
u3= log(n)+2*log(multi.opt$par[3])+2*log(1-multi.opt$par[3])
u4= log(n)+2*log(multi.opt$par[3])+2*log(1-multi.opt$par[3])
u5= log(n)+log(multi.opt$par[3])+3*log(1-multi.opt$par[3])
u6= log(n)+4*log(1-multi.opt$par[3])
#Initializaiton 
u1r= log(n)+4*log(paEst)
u2r= log(n)+3*log(paEst)+log(1-paEst)
u3r= log(n)+2*log(paEst)+2*log(1-paEst)
u4r= log(n)+2*log(paEst)+2*log(1-paEst)
u5r= log(n)+log(paEst)+3*log(1-paEst)
u6r= log(n)+4*log(1-paEst)

############################## ZIP ##########################
if(missing=="ZIP"){
source("fullZip.R")
source("reducedZip.R")

oldpar = c(u1,u2,u3,u4,u5,u6,multi.opt$par[1:3])
zip.opt <- try(optim(oldpar,f_zip,df_zip,c(MFCfunctions[1:15],
	M22,M21,M12,M11,M01,M10,M00),method="BFGS",
	control=list(fnscale=-1)))
if(class(zip.opt)=="try-error"){
	zip.na.full=TRUE
	zip.par.full=rep(NA,9)
	} else{
zip.par.full <- zip.opt$par
zip.l.full=zip.opt$value
zip.na.full=(zip.opt$converge>0)}

redpar = c(u1r,u2r,u3r,u4r,u5r,u6r,paEst)
zip.opt <- try(optim(redpar,f_red_zip,df_red_zip,c(MFCfunctions[
	1:15],M22,M21,M12,M11,M01,M10,M00),method="BFGS",
	control=list(fnscale=-1)))
if(class(zip.opt)=="try-error"){zip.na.red=TRUE} else{
zip.par.red <- zip.opt$par
zip.l.red=zip.opt$value
zip.na.red=(zip.opt$converge>0)}

zip.na = zip.na.full || zip.na.red
if (!(zip.na)){
zip.lrt=-2*(zip.l.red-zip.l.full)
rej <- (1-pchisq(zip.lrt,2))
beta <- zip.par.full
names(beta) <- c("mu1","mu2","mu3","mu4","mu5","mu6",
	"beta1","beta2","pa")
} else{
rej <- beta <- NULL
}
}
############################## Direct ##########################
if(missing=="Direct"){
source("full.R")
source("reduced.R")

oldpar = c(u1,u2,u3,u4,u5,u6,multi.opt$par[1:2])
direct.opt <- try(optim(oldpar,f_full,df_full,c(MFCfunctions[
	1:15],M22,M21,M12,M11,M01,M10,M00),method="BFGS",
	control=list(fnscale=-1)))
if(class(direct.opt)=="try-error"){
	direct.na.full=TRUE
	direct.par.full=rep(NA,8)
	} else{
direct.par.full <- direct.opt$par
direct.l.full=direct.opt$value
direct.na.full=(direct.opt$converge>0)}

redpar = c(u1r,u2r,u3r,u4r,u5r,u6r)
direct.opt <- try(optim(redpar,f_red,df_red,c(MFCfunctions[
	1:15],M22,M21,M12,M11,M01,M10,M00),method="BFGS",
	control=list(fnscale=-1)))
if(class(direct.opt)=="try-error"){direct.na.red=TRUE} else{
direct.par.red <- direct.opt$par
direct.l.red=direct.opt$value
direct.na.red=(direct.opt$converge>0)}

direct.na = direct.na.full || direct.na.red
if (!(direct.na)){
direct.lrt=-2*(direct.l.red-direct.l.full)
rej <- (1-pchisq(direct.lrt,2))
beta <- direct.par.full
names(beta) <- c("mu1","mu2","mu3","mu4","mu5","mu6",
	"beta1","beta2")
} else{
rej <- beta <- NULL}
	}
}
print(list(coef=beta,LRTp_value=rej,maf=paEst))
}