Nothing
R/cma.uni.mix.hl.R
In macc: Mediation Analysis of Causality under Confounding
cma.uni.mix.hl <-
function(dat,delta,tol=1e-4,max.itr=50,alpha=0.05,random.indep=TRUE,optimizer=c("bobyqa","Nelder_Mead","optimx"),
mix.pkg=c("lme4","nlme"),random.var.equal=TRUE,u.int=FALSE,Sigma.update=FALSE,var.constraint=FALSE,
random.var.update=TRUE)
{
N<-length(dat)
K<-length(dat[[1]])
n<-matrix(NA,N,K)
sigma1=sigma2<-matrix(NA,N,K)
At=Ct=Bt=C2t<-matrix(NA,N,K)
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
n[i,k]<-nrow(dd)
re2<-cma.uni.delta(dd,delta=delta)
sigma1[i,k]<-re2$Sigma[1,1]
sigma2[i,k]<-re2$Sigma[2,2]
At[i,k]<-re2$Coefficients[1,1]
Bt[i,k]<-re2$Coefficients[3,1]
Ct[i,k]<-re2$Coefficients[2,1]
C2t[i,k]<-re2$Coefficients[4,1]
}
}
coe<-rep(c("A","C","B"),each=N*K)
Sub<-factor(rep(1:N,K))
Sub.all<-rep(Sub,3)
coet<-c(c(At),c(Ct),c(Bt))
if(random.var.equal==FALSE)
{
if(random.indep==TRUE)
{
vec.At<-c(At)
vec.Bt<-c(Bt)
vec.Ct<-c(Ct)
if(mix.pkg[1]=="nlme")
{
if(optimizer[1]=="optimx")
{
fit.A<-lme(vec.At~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.B<-lme(vec.Bt~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.C<-lme(vec.Ct~1,random=~1|Sub,control=lmeControl(opt="optim"))
}else
{
fit.A<-lme(vec.At~1,random=~1|Sub)
fit.B<-lme(vec.Bt~1,random=~1|Sub)
fit.C<-lme(vec.Ct~1,random=~1|Sub)
}
Afix<-as.numeric(summary(fit.A)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit.B)$coefficients$fixed[1])
Cfix<-as.numeric(summary(fit.C)$coefficients$fixed[1])
Phi<-diag(c(as.numeric(VarCorr(fit.A)[1,1]),as.numeric(VarCorr(fit.B)[1,1]),as.numeric(VarCorr(fit.C)[1,1])))
Lambda<-diag(c(as.numeric(VarCorr(fit.A)[2,1]),as.numeric(VarCorr(fit.B)[2,1]),as.numeric(VarCorr(fit.C)[2,1])))
u<-cbind(as.matrix(ranef(fit.A)),as.matrix(ranef(fit.B)),as.matrix(ranef(fit.C)))
colnames(u)<-c("A","B","C")
b0<-c(Afix,Bfix,Cfix)
cor.AB=cor.AC=cor.BC<-0
LL<-as.numeric(logLik(fit.A,REML=FALSE)+logLik(fit.B,REML=FALSE)+logLik(fit.C,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-intervals(fit.A)
intB<-intervals(fit.B)
intC<-intervals(fit.C)
Phi.confint[1,]<-as.matrix(intA[[2]]$Sub[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(intB[[2]]$Sub[1,c(1,3)])^2
Phi.confint[3,]<-as.matrix(intC[[2]]$Sub[1,c(1,3)])^2
Lambda.confint[1,]<-as.matrix(intA[[3]][c(1,3)])^2
Lambda.confint[2,]<-as.matrix(intB[[3]][c(1,3)])^2
Lambda.confint[3,]<-as.matrix(intC[[3]][c(1,3)])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-intervals(fit.A)
intB<-intervals(fit.B)
intC<-intervals(fit.C)
Phi.confint[1,]<-as.matrix(intA[[2]]$Sub[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(intB[[2]]$Sub[1,c(1,3)])^2
Phi.confint[3,]<-as.matrix(intC[[2]]$Sub[1,c(1,3)])^2
Lambda.confint[1,]<-as.matrix(intA[[3]][c(1,3)])^2
Lambda.confint[2,]<-as.matrix(intB[[3]][c(1,3)])^2
Lambda.confint[3,]<-as.matrix(intC[[3]][c(1,3)])^2
}
}else
{
if(optimizer[1]=="optimx")
{
fit.A<-lmer(vec.At~1+(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
fit.B<-lmer(vec.Bt~1+(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
fit.C<-lmer(vec.Ct~1+(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
}else
{
fit.A<-lmer(vec.At~1+(1|Sub),control= lmerControl(optimizer=optimizer[1]))
fit.B<-lmer(vec.Bt~1+(1|Sub),control= lmerControl(optimizer=optimizer[1]))
fit.C<-lmer(vec.Ct~1+(1|Sub),control= lmerControl(optimizer=optimizer[1]))
}
Afix<-as.numeric(summary(fit.A)$coefficients[1])
Bfix<-as.numeric(summary(fit.B)$coefficients[1])
Cfix<-as.numeric(summary(fit.C)$coefficients[1])
Phi<-diag(c(as.numeric((attr(VarCorr(fit.A)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit.B)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit.C)[[1]],"stddev")^2)[1])))
Lambda<-diag(c(as.numeric((attr(VarCorr(fit.A),"sc")^2)[1]),
as.numeric((attr(VarCorr(fit.B),"sc")^2)[1]),
as.numeric((attr(VarCorr(fit.C),"sc")^2)[1])))
u<-cbind(as.matrix(ranef(fit.A)$Sub),as.matrix(ranef(fit.B)$Sub),as.matrix(ranef(fit.C)$Sub))
colnames(u)<-c("A","B","C")
b0<-c(Afix,Bfix,Cfix)
cor.AB=cor.AC=cor.BC<-0
LL<-as.numeric(logLik(fit.A,REML=FALSE)+logLik(fit.B,REML=FALSE)+logLik(fit.C,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-confint(fit.A)
intB<-confint(fit.B)
intC<-confint(fit.C)
Phi.confint[1,]<-as.matrix(intA[1,])^2
Phi.confint[2,]<-as.matrix(intB[1,])^2
Phi.confint[3,]<-as.matrix(intC[1,])^2
Lambda.confint[1,]<-as.matrix(intA[2,])^2
Lambda.confint[2,]<-as.matrix(intB[2,])^2
Lambda.confint[3,]<-as.matrix(intC[2,])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-confint(fit.A)
intB<-confint(fit.B)
intC<-confint(fit.C)
Phi.confint[1,]<-as.matrix(intA[1,])^2
Phi.confint[2,]<-as.matrix(intB[1,])^2
Phi.confint[3,]<-as.matrix(intC[1,])^2
Lambda.confint[1,]<-as.matrix(intA[2,])^2
Lambda.confint[2,]<-as.matrix(intB[2,])^2
Lambda.confint[3,]<-as.matrix(intC[2,])^2
}
}
}else
{
if(mix.pkg[1]=="nlme")
{
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
}
if(optimizer[1]=="optimx")
{
fit0<-lme(c(coet)~0+coe,random=~0+coe|Sub.all,control=lmeControl(opt="optim"))
}else
{
fit0<-lme(c(coet)~0+coe,random=~0+coe|Sub.all)
}
fit<-NULL
try(fit<-update(fit0,weights=varIdent(form=~1|factor(coe))))
if(is.null(fit)==TRUE)
{
warning("Equal-varaince assumption in random error is applied instead.")
fit<-fit0
Lambda<-diag(rep(summary(fit)$sigma^2,3))
if(!is.matrix(var.constraint))
{
if(var.constraint)
{
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-(intervals(fit)[[3]][c(1,3)])^2
}
}
}else
{
lambda.tmp<-summary(fit)$sigma*coef(fit$modelStruct$varStruct, unconstrained=FALSE, allCoef=TRUE)
Lambda<-diag(c((lambda.tmp[1])^2,(lambda.tmp[3])^2,(lambda.tmp[2])^2))
if(is.matrix(var.constraint)==FALSE)
{
if(var.constraint==TRUE)
{
int.fit<-intervals(fit)
Lambda.confint[1,]<-(int.fit[[4]][c(1,3)])^2
Lambda.confint[2,]<-(int.fit[[3]][2,c(1,3)])^2*Lambda[1,1]
Lambda.confint[3,]<-(int.fit[[3]][1,c(1,3)])^2*Lambda[1,1]
}
}
}
Afix<-as.numeric(summary(fit)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit)$coefficients$fixed[2])
Cfix<-as.numeric(summary(fit)$coefficients$fixed[3])
b0<-c(Afix,Bfix,Cfix)
u<-as.matrix(ranef(fit))
Phi<-diag(c(as.numeric(VarCorr(fit)[1,1]),as.numeric(VarCorr(fit)[2,1]),as.numeric(VarCorr(fit)[3,1])))
cor.AB<-as.numeric(VarCorr(fit)[2,3])
cor.AC<-as.numeric(VarCorr(fit)[3,3])
cor.BC<-as.numeric(VarCorr(fit)[3,4])
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
LL<-as.numeric(logLik(fit,REML=FALSE))
if(!is.matrix(var.constraint))
{
if(var.constraint)
{
int.fit<-intervals(fit)
Phi.confint[1,]<-as.matrix(int.fit[[2]]$Sub.all[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(int.fit[[2]]$Sub.all[2,c(1,3)])^2
Phi.confint[3,]<-as.matrix(int.fit[[2]]$Sub.all[3,c(1,3)])^2
}
}
}else
{
warning("Equal-varaince assumption in random error is applied instead.")
if(optimizer[1]=="optimx")
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
}else
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer=optimizer[1]))
}
u<-as.matrix(ranef(fit)$Sub.all)
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
s<-0
while(length(which(abs(cor.t)<1e-06))>0&s<20)
{
s<-s+1
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
}
b0<-c(summary(fit)$coefficients[1,1],summary(fit)$coefficients[2,1],summary(fit)$coefficients[3,1])
Phi<-diag(c(as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[2]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[3])))
cor.AB<-attr(VarCorr(fit)[[1]],"correlation")[1,2]
cor.AC<-attr(VarCorr(fit)[[1]],"correlation")[1,3]
cor.BC<-attr(VarCorr(fit)[[1]],"correlation")[2,3]
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
Lambda<-diag(rep(attr(VarCorr(fit),"sc")^2,3))
LL<-as.numeric(logLik(fit,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-as.matrix(re.confint[1,])^2
Phi.confint[2,]<-as.matrix(re.confint[4,])^2
Phi.confint[3,]<-as.matrix(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(re.confint[7,])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-as.matrix(re.confint[1,])^2
Phi.confint[2,]<-as.matrix(re.confint[4,])^2
Phi.confint[3,]<-as.matrix(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(re.confint[7,])^2
}
}
}
}else
{
if(mix.pkg[1]=="nlme")
{
if(optimizer[1]=="optimx")
{
fit<-lme(c(coet)~0+coe,random=~0+coe|Sub.all,control=lmeControl(opt="optim"))
}else
{
fit<-lme(c(coet)~0+coe,random=~0+coe|Sub.all)
}
Afix<-as.numeric(summary(fit)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit)$coefficients$fixed[2])
Cfix<-as.numeric(summary(fit)$coefficients$fixed[3])
b0<-c(Afix,Bfix,Cfix)
u<-as.matrix(ranef(fit))
Lambda<-diag(rep(summary(fit)$sigma^2,3))
Phi<-diag(c(as.numeric(VarCorr(fit)[1,1]),as.numeric(VarCorr(fit)[2,1]),as.numeric(VarCorr(fit)[3,1])))
if(random.indep)
{
cor.AB=cor.AC=cor.BC<-0
}else
{
cor.AB<-as.numeric(VarCorr(fit)[2,3])
cor.AC<-as.numeric(VarCorr(fit)[3,3])
cor.BC<-as.numeric(VarCorr(fit)[3,4])
}
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
LL<-as.numeric(logLik(fit,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)>=4)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-matrix(rep(var.constraint[4,],3),nrow=3,ncol=2,byrow=TRUE)
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
int.fit<-intervals(fit)
Phi.confint[1,]<-as.matrix(int.fit[[2]]$Sub.all[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(int.fit[[2]]$Sub.all[2,c(1,3)])^2
Phi.confint[3,]<-as.matrix(int.fit[[2]]$Sub.all[3,c(1,3)])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(int.fit[[3]][c(1,3)])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
int.fit<-intervals(fit)
Phi.confint[1,]<-as.matrix(int.fit[[2]]$Sub.all[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(int.fit[[2]]$Sub.all[2,c(1,3)])^2
Phi.confint[3,]<-as.matrix(int.fit[[2]]$Sub.all[3,c(1,3)])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(intervals(fit)[[3]][c(1,3)])^2
}
}else
{
if(optimizer[1]=="optimx")
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
u<-as.matrix(ranef(fit)$Sub.all)
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
s<-0
while(length(which(abs(cor.t)<1e-06))>0&s<20)
{
s<-s+1
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
}
}else
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control=lmerControl(optimizer=optimizer[1]))
u<-as.matrix(ranef(fit)$Sub.all)
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
s<-0
while(length(which(abs(cor.t)<1e-06))>0&s<20)
{
s<-s+1
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
}
}
b0<-c(summary(fit)$coefficients[1,1],summary(fit)$coefficients[2,1],summary(fit)$coefficients[3,1])
Phi<-diag(c(as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[2]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[3])))
if(random.indep==TRUE)
{
cor.AB=cor.AC=cor.BC<-0
}else
{
cor.AB<-attr(VarCorr(fit)[[1]],"correlation")[1,2]
cor.AC<-attr(VarCorr(fit)[[1]],"correlation")[1,3]
cor.BC<-attr(VarCorr(fit)[[1]],"correlation")[2,3]
}
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
Lambda<-diag(rep(attr(VarCorr(fit),"sc")^2,3))
LL<-as.numeric(logLik(fit,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)>=4)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-matrix(rep(var.constraint[4,],3),nrow=3,ncol=2,byrow=TRUE)
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-(re.confint[1,])^2
Phi.confint[2,]<-(re.confint[4,])^2
Phi.confint[3,]<-(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-(re.confint[7,])^2
}
}else
if(var.constraint==TRUE)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-(re.confint[1,])^2
Phi.confint[2,]<-(re.confint[4,])^2
Phi.confint[3,]<-(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-(re.confint[7,])^2
}
}
}
diff<-100
s<-0
while(s<max.itr&diff>=tol)
{
if(random.indep&!random.var.equal&u.int)
{
V.alpha<-matrix(Phi[1,1],K,K)+diag(rep(Lambda[1,1],K))
V.beta<-matrix(Phi[2,2],K,K)+diag(rep(Lambda[2,2],K))
V.gamma<-matrix(Phi[3,3],K,K)+diag(rep(Lambda[3,3],K))
for(i in 1:N)
{
# A
z1z=z1m=z1z.d=z1m.d<-NULL
for(k in 1:K)
{
dd<-dat[[i]][[k]]
inv.O1<-diag(rep(1/sigma1[i,k],nrow(dd)))
inv.O2<-diag(rep(1/sigma2[i,k],nrow(dd)))
O12<-diag(rep(sqrt(sigma2[i,k]/sigma1[i,k]),nrow(dd)))
z1z<-c(z1z,(t(dd$Z)%*%inv.O1%*%(dd$Z))[1,1])
z1m<-c(z1m,(t(dd$Z)%*%inv.O1%*%(dd$M))[1,1])
z1z.d<-c(z1z.d,(delta^2*t(dd$Z)%*%O12%*%inv.O2%*%O12%*%(dd$Z)/(1-delta^2))[1,1])
z1m.d<-c(z1m.d,(delta*t(dd$Z)%*%O12%*%inv.O2%*%(dd$R-dd$Z*Ct[i,k]-dd$M*Bt[i,k]-delta*O12%*%dd$M)/(1-delta^2))[1,1])
}
At[i,]<-solve(diag(z1z+z1z.d)+solve(V.alpha))%*%(z1m-z1m.d+solve(V.alpha)%*%rep(b0[1],K))
# B
m2m=m2r.d<-NULL
for(k in 1:K)
{
dd<-dat[[i]][[k]]
inv.O1<-diag(rep(1/sigma1[i,k],nrow(dd)))
inv.O2<-diag(rep(1/sigma2[i,k],nrow(dd)))
O12<-diag(rep(sqrt(sigma2[i,k]/sigma1[i,k]),nrow(dd)))
m2m<-c(m2m,(t(dd$M)%*%inv.O2%*%(dd$M)/(1-delta^2))[1,1])
m2r.d<-c(m2r.d,(t(dd$M)%*%inv.O2%*%(dd$R-dd$Z*Ct[i,k]-delta*O12%*%(dd$M-dd$Z*At[i,k]))/(1-delta^2))[1,1])
}
Bt[i,]<-solve(diag(m2m)+solve(V.beta))%*%(m2r.d+solve(V.beta)%*%rep(b0[2],K))
# C
z2z=z2r.d<-NULL
for(k in 1:K)
{
dd<-dat[[i]][[k]]
inv.O1<-diag(rep(1/sigma1[i,k],nrow(dd)))
inv.O2<-diag(rep(1/sigma2[i,k],nrow(dd)))
O12<-diag(rep(sqrt(sigma2[i,k]/sigma1[i,k]),nrow(dd)))
z2z<-c(z2z,(t(dd$Z)%*%inv.O2%*%(dd$Z)/(1-delta^2))[1,1])
z2r.d<-c(z2r.d,(t(dd$Z)%*%inv.O2%*%(dd$R-dd$M*Bt[i,k]-delta*O12%*%(dd$M-dd$Z*At[i,k]))/(1-delta^2))[1,1])
}
Ct[i,]<-solve(diag(z2z)+solve(V.gamma))%*%(z2r.d+solve(V.gamma)%*%rep(b0[3],K))
}
# update sigma_1ik and sigma_2ik
if(Sigma.update)
{
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
Y<-cbind(dd$M,dd$R)
X<-cbind(dd$Z,dd$M)
Theta<-matrix(c(At[i,k],0,Ct[i,k],Bt[i,k]),2,2)
S<-t(Y-X%*%Theta)%*%(Y-X%*%Theta)
sigma1[i,k]<-(S[1,1]-delta*S[1,2]*sqrt(S[1,1]/S[2,2]))/(nrow(dd)*(1-delta^2))
sigma2[i,k]<-(S[2,2]-delta*S[1,2]*sqrt(S[2,2]/S[1,1]))/(nrow(dd)*(1-delta^2))
}
}
}
vec.At<-c(At)
vec.Bt<-c(Bt)
vec.Ct<-c(Ct)
fit.A<-lme(vec.At~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.B<-lme(vec.Bt~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.C<-lme(vec.Ct~1,random=~1|Sub,control=lmeControl(opt="optim"))
Afix<-as.numeric(summary(fit.A)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit.B)$coefficients$fixed[1])
Cfix<-as.numeric(summary(fit.C)$coefficients$fixed[1])
if(random.var.update)
{
Phi<-diag(c(as.numeric(VarCorr(fit.A)[1,1]),as.numeric(VarCorr(fit.B)[1,1]),as.numeric(VarCorr(fit.C)[1,1])))
Lambda<-diag(c(as.numeric(VarCorr(fit.A)[2,1]),as.numeric(VarCorr(fit.B)[2,1]),as.numeric(VarCorr(fit.C)[2,1])))
if(sum(var.constraint==FALSE)==0)
{
# Phi
if(Phi[1,1]<Phi.confint[1,1])
{
Phi[1,1]<-Phi.confint[1,1]
}else
if(Phi[1,1]>Phi.confint[1,2])
{
Phi[1,1]<-Phi.confint[1,2]
}
if(Phi[2,2]<Phi.confint[2,1])
{
Phi[2,2]<-Phi.confint[2,1]
}else
if(Phi[2,2]>Phi.confint[2,2])
{
Phi[2,2]<-Phi.confint[2,2]
}
if(Phi[3,3]<Phi.confint[3,1])
{
Phi[3,3]<-Phi.confint[3,1]
}else
if(Phi[3,3]>Phi.confint[3,2])
{
Phi[3,3]<-Phi.confint[3,2]
}
# Lambda
if(Lambda[1,1]<Lambda.confint[1,1])
{
Lambda[1,1]<-Lambda.confint[1,1]
}else
if(Lambda[1,1]>Lambda.confint[1,2])
{
Lambda[1,1]<-Lambda.confint[1,2]
}
if(Lambda[2,2]<Lambda.confint[2,1])
{
Lambda[2,2]<-Lambda.confint[2,1]
}else
if(Lambda[2,2]>Lambda.confint[2,2])
{
Lambda[2,2]<-Lambda.confint[2,2]
}
if(Lambda[3,3]<Lambda.confint[3,1])
{
Lambda[3,3]<-Lambda.confint[3,1]
}else
if(Lambda[3,3]>Lambda.confint[3,2])
{
Lambda[3,3]<-Lambda.confint[3,2]
}
}
}
u<-cbind(as.matrix(ranef(fit.A)),as.matrix(ranef(fit.B)),as.matrix(ranef(fit.C)))
colnames(u)<-c("A","B","C")
bnew<-c(Afix,Bfix,Cfix)
diff<-max(abs(bnew-b0))
b0<-bnew
s<-s+1
}else
{
if((random.var.equal&u.int)|(!random.indep&u.int))
{
warning("The u.int=TRUE argument will be ignored and the full h-likelihood method will be applied.")
}
bik.mat<-array(NA,c(N,K,3))
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
colnames(dd)<-c("Z","M","R")
n[i,k]<-nrow(dd)
Z<-matrix(dd[,1])
M<-matrix(dd[,2])
R<-matrix(dd[,3])
r<-delta*sqrt(sigma2[i,k]/sigma1[i,k])
X<-cbind(Z,M)
Pik<-matrix(c(-r,0,0,1,1,0),2,3)
Qik<-c(0,r)
Vik<-matrix(c(1,0,0),1,3)
bik<-solve(sigma1[i,k]*t(Pik)%*%t(X)%*%X%*%Pik+sigma2[i,k]*(1-delta^2)*t(Vik)%*%t(Z)%*%Z%*%Vik+
sigma1[i,k]*sigma2[i,k]*(1-delta^2)*solve(Lambda))%*%
(sigma1[i,k]*t(Pik)%*%t(X)%*%(R-X%*%Qik)+sigma2[i,k]*(1-delta^2)*t(Vik)%*%t(Z)%*%M+
sigma1[i,k]*sigma2[i,k]*(1-delta^2)*solve(Lambda)%*%(b0+u[i,]))
At[i,k]<-bik[1,1]
Bt[i,k]<-bik[2,1]
Ct[i,k]<-bik[3,1]
bik.mat[i,k,]<-bik
}
}
e<-c(sum(At-b0[1])/N,sum(Bt-b0[2])/N,sum(Ct-b0[3])/N)
for(i in 1:N)
{
u[i,]<-solve(solve(Phi)+K*solve(Lambda))%*%solve(Lambda)%*%(apply(cbind(At[i,],Bt[i,],Ct[i,]),2,sum)-K*b0-e)
}
bnew<-c(mean(At),mean(Bt),mean(Ct))-apply(u,2,mean)
if(random.var.update==TRUE)
{
Phi<-t(u)%*%u/N
if(det(Phi)<10e-5|random.indep==TRUE)
{
Phi<-diag(as.vector(diag(Phi)))
}
if(random.var.equal==FALSE)
{
Lambda.new<-matrix(0,3,3)
for(i in 1:N)
{
for(k in 1:K)
{
Lambda.new<-Lambda.new+(bik.mat[i,k,]-bnew-u[i,])%*%t(bik.mat[i,k,]-bnew-u[i,])/(N*K)
}
}
Lambda<-diag(as.vector(diag(Lambda.new)))
}else
{
theta2<-0
for(i in 1:N)
{
for(k in 1:K)
{
theta2<-theta2+(t(bik.mat[i,k,]-bnew-u[i,])%*%(bik.mat[i,k,]-bnew-u[i,])/(3*N*K))[1,1]
}
}
Lambda<-diag(rep(theta2,3))
}
Phi<-t(u)%*%u/N
if(random.indep)
{
Phi<-diag(as.vector(diag(Phi)))
}
# add variance constraint
if(sum(var.constraint==FALSE)==0)
{
# Phi
if(Phi[1,1]<Phi.confint[1,1])
{
Phi[1,1]<-Phi.confint[1,1]
}else
if(Phi[1,1]>Phi.confint[1,2])
{
Phi[1,1]<-Phi.confint[1,2]
}
if(Phi[2,2]<Phi.confint[2,1])
{
Phi[2,2]<-Phi.confint[2,1]
}else
if(Phi[2,2]>Phi.confint[2,2])
{
Phi[2,2]<-Phi.confint[2,2]
}
if(Phi[3,3]<Phi.confint[3,1])
{
Phi[3,3]<-Phi.confint[3,1]
}else
if(Phi[3,3]>Phi.confint[3,2])
{
Phi[3,3]<-Phi.confint[3,2]
}
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
# Lambda
if(Lambda[1,1]<Lambda.confint[1,1])
{
Lambda[1,1]<-Lambda.confint[1,1]
}else
if(Lambda[1,1]>Lambda.confint[1,2])
{
Lambda[1,1]<-Lambda.confint[1,2]
}
if(Lambda[2,2]<Lambda.confint[2,1])
{
Lambda[2,2]<-Lambda.confint[2,1]
}else
if(Lambda[2,2]>Lambda.confint[2,2])
{
Lambda[2,2]<-Lambda.confint[2,2]
}
if(Lambda[3,3]<Lambda.confint[3,1])
{
Lambda[3,3]<-Lambda.confint[3,1]
}else
if(Lambda[3,3]>Lambda.confint[3,2])
{
Lambda[3,3]<-Lambda.confint[3,2]
}
}
}
# update sigma_1ik and sigma_2ik
if(Sigma.update)
{
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
Y<-cbind(dd$M,dd$R)
X<-cbind(dd$Z,dd$M)
Theta<-matrix(c(At[i,k],0,Ct[i,k],Bt[i,k]),2,2)
S<-t(Y-X%*%Theta)%*%(Y-X%*%Theta)
sigma1[i,k]<-(S[1,1]-delta*S[1,2]*sqrt(S[1,1]/S[2,2]))/(nrow(dd)*(1-delta^2))
sigma2[i,k]<-(S[2,2]-delta*S[1,2]*sqrt(S[2,2]/S[1,1]))/(nrow(dd)*(1-delta^2))
}
}
}
diff<-max(abs(b0-bnew))
b0<-bnew
s<-s+1
}
}
HL<-cma.h(dat,delta=delta,A.ik=At,B.ik=Bt,C.ik=Ct,b=b0,u=u,Phi=Phi,Lambda=Lambda,random.indep=random.indep,
u.int=u.int,Sigma.update=Sigma.update)$h
zc<-qnorm(1-alpha/2)
coe.re<-matrix(NA,6,4)
colnames(coe.re)<-c("Estimate","LB","UB","SE")
rownames(coe.re)<-c("A","C","B","C2","AB.prod","AB.diff")
# A, C, B, AB.prod, AB.diff
coe.re[1,4]<-sqrt(Phi[1,1]/N+Lambda[1,1]/(N*K))
coe.re[2,4]<-sqrt(Phi[3,3]/N+Lambda[3,3]/(N*K))
coe.re[3,4]<-sqrt(Phi[2,2]/N+Lambda[2,2]/(N*K))
coe.re[1,1:3]<-c(b0[1],b0[1]-zc*coe.re[1,4],b0[1]+zc*coe.re[1,4])
coe.re[2,1:3]<-c(b0[3],b0[3]-zc*coe.re[2,4],b0[3]+zc*coe.re[2,4])
coe.re[3,1:3]<-c(b0[2],b0[2]-zc*coe.re[3,4],b0[2]+zc*coe.re[3,4])
coe.re[5,1]<-b0[1]*b0[2]
coe.re[5,4]<-sqrt((b0[1]*coe.re[3,4])^2+(b0[2]*coe.re[1,4])^2)
coe.re[5,2:3]<-c(coe.re[5,1]-zc*coe.re[5,4],coe.re[5,1]+zc*coe.re[5,4])
# C2, AB.diff
if(optimizer[1]=="optimx")
{
fit.C2<-lmer(c(C2t)~(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
}else
{
fit.C2<-lmer(c(C2t)~(1|Sub),control=lmerControl(optimizer=optimizer[1]))
}
C2fix<-summary(fit.C2)$coefficients[1]
s2.C2<-c(as.numeric(VarCorr(fit.C2)),attr(VarCorr(fit.C2),"sc")^2)
e.C2<-C2fix+mean(ranef(fit.C2)$Sub[,"(Intercept)"])
coe.re[4,4]<-sqrt(s2.C2[1]/N+s2.C2[2]/(N*K))
coe.re[4,1:3]<-c(C2fix,C2fix-zc*coe.re[4,4],C2fix+zc*coe.re[4,4])
coe.re[6,1]<-coe.re[4,1]-coe.re[2,1]
coe.re[6,4]<-sqrt(coe.re[4,4]^2+coe.re[2,4]^2)
coe.re[6,2:3]<-c(coe.re[6,1]-zc*coe.re[6,4],coe.re[6,1]+zc*coe.re[6,4])
s2.C2<-data.frame(C2=s2.C2[1],Error=s2.C2[2],Total=sum(s2.C2))
colnames(s2.C2)<-c("Random(C')","Var(Error)","Var(C')")
var.comp<-data.frame(delta=delta,A=Phi[1,1],C=Phi[3,3],B=Phi[2,2],Lambda[1,1],Lambda[3,3],Lambda[2,2])
colnames(var.comp)<-c("delta","Random(A)","Random(C)","Random(B)","Var(Error A)","Var(Error C)","Var(Error B)")
if(random.indep)
{
Phi<-diag(as.vector(diag(Phi)))
}
cor.comp.temp<-sqrt(diag(1/diag(Phi)))%*%Phi%*%sqrt(diag(1/diag(Phi)))
cor.comp<-cor.comp.temp
cor.comp[1,2]=cor.comp[2,1]<-cor.comp.temp[1,3]
cor.comp[1,3]=cor.comp[3,1]<-cor.comp.temp[1,2]
colnames(cor.comp)=rownames(cor.comp)<-c("A","C","B")
re<-list(delta=delta,Coefficients=coe.re,Cor.comp=cor.comp,Var.comp=var.comp,Var.C2=s2.C2,logLik=LL,HL=HL)
if(sum(var.constraint==FALSE)==0)
{
re.var.constraint<-rbind(Phi.confint,Lambda.confint)
rownames(re.var.constraint)<-c("Random(A)","Random(B)","Random(C)","Var(Error A)","Var(Error B)","Var(Error C)")
re<-list(delta=delta,Coefficients=coe.re,Cor.comp=cor.comp,Var.comp=var.comp,Var.constraint=re.var.constraint,
Var.C2=s2.C2,logLik=LL,HL=HL,convergence=(s<max.itr|max.itr==0))
}
return(re)
}
Try the macc package in your browser
Any scripts or data that you put into this service are public.
macc documentation built on May 2, 2019, 12:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
cma.uni.mix.hl <-
function(dat,delta,tol=1e-4,max.itr=50,alpha=0.05,random.indep=TRUE,optimizer=c("bobyqa","Nelder_Mead","optimx"),
mix.pkg=c("lme4","nlme"),random.var.equal=TRUE,u.int=FALSE,Sigma.update=FALSE,var.constraint=FALSE,
random.var.update=TRUE)
{
N<-length(dat)
K<-length(dat[[1]])
n<-matrix(NA,N,K)
sigma1=sigma2<-matrix(NA,N,K)
At=Ct=Bt=C2t<-matrix(NA,N,K)
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
n[i,k]<-nrow(dd)
re2<-cma.uni.delta(dd,delta=delta)
sigma1[i,k]<-re2$Sigma[1,1]
sigma2[i,k]<-re2$Sigma[2,2]
At[i,k]<-re2$Coefficients[1,1]
Bt[i,k]<-re2$Coefficients[3,1]
Ct[i,k]<-re2$Coefficients[2,1]
C2t[i,k]<-re2$Coefficients[4,1]
}
}
coe<-rep(c("A","C","B"),each=N*K)
Sub<-factor(rep(1:N,K))
Sub.all<-rep(Sub,3)
coet<-c(c(At),c(Ct),c(Bt))
if(random.var.equal==FALSE)
{
if(random.indep==TRUE)
{
vec.At<-c(At)
vec.Bt<-c(Bt)
vec.Ct<-c(Ct)
if(mix.pkg[1]=="nlme")
{
if(optimizer[1]=="optimx")
{
fit.A<-lme(vec.At~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.B<-lme(vec.Bt~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.C<-lme(vec.Ct~1,random=~1|Sub,control=lmeControl(opt="optim"))
}else
{
fit.A<-lme(vec.At~1,random=~1|Sub)
fit.B<-lme(vec.Bt~1,random=~1|Sub)
fit.C<-lme(vec.Ct~1,random=~1|Sub)
}
Afix<-as.numeric(summary(fit.A)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit.B)$coefficients$fixed[1])
Cfix<-as.numeric(summary(fit.C)$coefficients$fixed[1])
Phi<-diag(c(as.numeric(VarCorr(fit.A)[1,1]),as.numeric(VarCorr(fit.B)[1,1]),as.numeric(VarCorr(fit.C)[1,1])))
Lambda<-diag(c(as.numeric(VarCorr(fit.A)[2,1]),as.numeric(VarCorr(fit.B)[2,1]),as.numeric(VarCorr(fit.C)[2,1])))
u<-cbind(as.matrix(ranef(fit.A)),as.matrix(ranef(fit.B)),as.matrix(ranef(fit.C)))
colnames(u)<-c("A","B","C")
b0<-c(Afix,Bfix,Cfix)
cor.AB=cor.AC=cor.BC<-0
LL<-as.numeric(logLik(fit.A,REML=FALSE)+logLik(fit.B,REML=FALSE)+logLik(fit.C,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-intervals(fit.A)
intB<-intervals(fit.B)
intC<-intervals(fit.C)
Phi.confint[1,]<-as.matrix(intA[[2]]$Sub[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(intB[[2]]$Sub[1,c(1,3)])^2
Phi.confint[3,]<-as.matrix(intC[[2]]$Sub[1,c(1,3)])^2
Lambda.confint[1,]<-as.matrix(intA[[3]][c(1,3)])^2
Lambda.confint[2,]<-as.matrix(intB[[3]][c(1,3)])^2
Lambda.confint[3,]<-as.matrix(intC[[3]][c(1,3)])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-intervals(fit.A)
intB<-intervals(fit.B)
intC<-intervals(fit.C)
Phi.confint[1,]<-as.matrix(intA[[2]]$Sub[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(intB[[2]]$Sub[1,c(1,3)])^2
Phi.confint[3,]<-as.matrix(intC[[2]]$Sub[1,c(1,3)])^2
Lambda.confint[1,]<-as.matrix(intA[[3]][c(1,3)])^2
Lambda.confint[2,]<-as.matrix(intB[[3]][c(1,3)])^2
Lambda.confint[3,]<-as.matrix(intC[[3]][c(1,3)])^2
}
}else
{
if(optimizer[1]=="optimx")
{
fit.A<-lmer(vec.At~1+(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
fit.B<-lmer(vec.Bt~1+(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
fit.C<-lmer(vec.Ct~1+(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
}else
{
fit.A<-lmer(vec.At~1+(1|Sub),control= lmerControl(optimizer=optimizer[1]))
fit.B<-lmer(vec.Bt~1+(1|Sub),control= lmerControl(optimizer=optimizer[1]))
fit.C<-lmer(vec.Ct~1+(1|Sub),control= lmerControl(optimizer=optimizer[1]))
}
Afix<-as.numeric(summary(fit.A)$coefficients[1])
Bfix<-as.numeric(summary(fit.B)$coefficients[1])
Cfix<-as.numeric(summary(fit.C)$coefficients[1])
Phi<-diag(c(as.numeric((attr(VarCorr(fit.A)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit.B)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit.C)[[1]],"stddev")^2)[1])))
Lambda<-diag(c(as.numeric((attr(VarCorr(fit.A),"sc")^2)[1]),
as.numeric((attr(VarCorr(fit.B),"sc")^2)[1]),
as.numeric((attr(VarCorr(fit.C),"sc")^2)[1])))
u<-cbind(as.matrix(ranef(fit.A)$Sub),as.matrix(ranef(fit.B)$Sub),as.matrix(ranef(fit.C)$Sub))
colnames(u)<-c("A","B","C")
b0<-c(Afix,Bfix,Cfix)
cor.AB=cor.AC=cor.BC<-0
LL<-as.numeric(logLik(fit.A,REML=FALSE)+logLik(fit.B,REML=FALSE)+logLik(fit.C,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-confint(fit.A)
intB<-confint(fit.B)
intC<-confint(fit.C)
Phi.confint[1,]<-as.matrix(intA[1,])^2
Phi.confint[2,]<-as.matrix(intB[1,])^2
Phi.confint[3,]<-as.matrix(intC[1,])^2
Lambda.confint[1,]<-as.matrix(intA[2,])^2
Lambda.confint[2,]<-as.matrix(intB[2,])^2
Lambda.confint[3,]<-as.matrix(intC[2,])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
intA<-confint(fit.A)
intB<-confint(fit.B)
intC<-confint(fit.C)
Phi.confint[1,]<-as.matrix(intA[1,])^2
Phi.confint[2,]<-as.matrix(intB[1,])^2
Phi.confint[3,]<-as.matrix(intC[1,])^2
Lambda.confint[1,]<-as.matrix(intA[2,])^2
Lambda.confint[2,]<-as.matrix(intB[2,])^2
Lambda.confint[3,]<-as.matrix(intC[2,])^2
}
}
}else
{
if(mix.pkg[1]=="nlme")
{
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
}
if(optimizer[1]=="optimx")
{
fit0<-lme(c(coet)~0+coe,random=~0+coe|Sub.all,control=lmeControl(opt="optim"))
}else
{
fit0<-lme(c(coet)~0+coe,random=~0+coe|Sub.all)
}
fit<-NULL
try(fit<-update(fit0,weights=varIdent(form=~1|factor(coe))))
if(is.null(fit)==TRUE)
{
warning("Equal-varaince assumption in random error is applied instead.")
fit<-fit0
Lambda<-diag(rep(summary(fit)$sigma^2,3))
if(!is.matrix(var.constraint))
{
if(var.constraint)
{
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-(intervals(fit)[[3]][c(1,3)])^2
}
}
}else
{
lambda.tmp<-summary(fit)$sigma*coef(fit$modelStruct$varStruct, unconstrained=FALSE, allCoef=TRUE)
Lambda<-diag(c((lambda.tmp[1])^2,(lambda.tmp[3])^2,(lambda.tmp[2])^2))
if(is.matrix(var.constraint)==FALSE)
{
if(var.constraint==TRUE)
{
int.fit<-intervals(fit)
Lambda.confint[1,]<-(int.fit[[4]][c(1,3)])^2
Lambda.confint[2,]<-(int.fit[[3]][2,c(1,3)])^2*Lambda[1,1]
Lambda.confint[3,]<-(int.fit[[3]][1,c(1,3)])^2*Lambda[1,1]
}
}
}
Afix<-as.numeric(summary(fit)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit)$coefficients$fixed[2])
Cfix<-as.numeric(summary(fit)$coefficients$fixed[3])
b0<-c(Afix,Bfix,Cfix)
u<-as.matrix(ranef(fit))
Phi<-diag(c(as.numeric(VarCorr(fit)[1,1]),as.numeric(VarCorr(fit)[2,1]),as.numeric(VarCorr(fit)[3,1])))
cor.AB<-as.numeric(VarCorr(fit)[2,3])
cor.AC<-as.numeric(VarCorr(fit)[3,3])
cor.BC<-as.numeric(VarCorr(fit)[3,4])
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
LL<-as.numeric(logLik(fit,REML=FALSE))
if(!is.matrix(var.constraint))
{
if(var.constraint)
{
int.fit<-intervals(fit)
Phi.confint[1,]<-as.matrix(int.fit[[2]]$Sub.all[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(int.fit[[2]]$Sub.all[2,c(1,3)])^2
Phi.confint[3,]<-as.matrix(int.fit[[2]]$Sub.all[3,c(1,3)])^2
}
}
}else
{
warning("Equal-varaince assumption in random error is applied instead.")
if(optimizer[1]=="optimx")
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
}else
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer=optimizer[1]))
}
u<-as.matrix(ranef(fit)$Sub.all)
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
s<-0
while(length(which(abs(cor.t)<1e-06))>0&s<20)
{
s<-s+1
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
}
b0<-c(summary(fit)$coefficients[1,1],summary(fit)$coefficients[2,1],summary(fit)$coefficients[3,1])
Phi<-diag(c(as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[2]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[3])))
cor.AB<-attr(VarCorr(fit)[[1]],"correlation")[1,2]
cor.AC<-attr(VarCorr(fit)[[1]],"correlation")[1,3]
cor.BC<-attr(VarCorr(fit)[[1]],"correlation")[2,3]
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
Lambda<-diag(rep(attr(VarCorr(fit),"sc")^2,3))
LL<-as.numeric(logLik(fit,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)==6)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-var.constraint[4:6,]
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-as.matrix(re.confint[1,])^2
Phi.confint[2,]<-as.matrix(re.confint[4,])^2
Phi.confint[3,]<-as.matrix(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(re.confint[7,])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-as.matrix(re.confint[1,])^2
Phi.confint[2,]<-as.matrix(re.confint[4,])^2
Phi.confint[3,]<-as.matrix(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(re.confint[7,])^2
}
}
}
}else
{
if(mix.pkg[1]=="nlme")
{
if(optimizer[1]=="optimx")
{
fit<-lme(c(coet)~0+coe,random=~0+coe|Sub.all,control=lmeControl(opt="optim"))
}else
{
fit<-lme(c(coet)~0+coe,random=~0+coe|Sub.all)
}
Afix<-as.numeric(summary(fit)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit)$coefficients$fixed[2])
Cfix<-as.numeric(summary(fit)$coefficients$fixed[3])
b0<-c(Afix,Bfix,Cfix)
u<-as.matrix(ranef(fit))
Lambda<-diag(rep(summary(fit)$sigma^2,3))
Phi<-diag(c(as.numeric(VarCorr(fit)[1,1]),as.numeric(VarCorr(fit)[2,1]),as.numeric(VarCorr(fit)[3,1])))
if(random.indep)
{
cor.AB=cor.AC=cor.BC<-0
}else
{
cor.AB<-as.numeric(VarCorr(fit)[2,3])
cor.AC<-as.numeric(VarCorr(fit)[3,3])
cor.BC<-as.numeric(VarCorr(fit)[3,4])
}
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
LL<-as.numeric(logLik(fit,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)>=4)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-matrix(rep(var.constraint[4,],3),nrow=3,ncol=2,byrow=TRUE)
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
int.fit<-intervals(fit)
Phi.confint[1,]<-as.matrix(int.fit[[2]]$Sub.all[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(int.fit[[2]]$Sub.all[2,c(1,3)])^2
Phi.confint[3,]<-as.matrix(int.fit[[2]]$Sub.all[3,c(1,3)])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(int.fit[[3]][c(1,3)])^2
}
}else
if(var.constraint)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
int.fit<-intervals(fit)
Phi.confint[1,]<-as.matrix(int.fit[[2]]$Sub.all[1,c(1,3)])^2
Phi.confint[2,]<-as.matrix(int.fit[[2]]$Sub.all[2,c(1,3)])^2
Phi.confint[3,]<-as.matrix(int.fit[[2]]$Sub.all[3,c(1,3)])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-as.matrix(intervals(fit)[[3]][c(1,3)])^2
}
}else
{
if(optimizer[1]=="optimx")
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
u<-as.matrix(ranef(fit)$Sub.all)
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
s<-0
while(length(which(abs(cor.t)<1e-06))>0&s<20)
{
s<-s+1
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
}
}else
{
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control=lmerControl(optimizer=optimizer[1]))
u<-as.matrix(ranef(fit)$Sub.all)
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
s<-0
while(length(which(abs(cor.t)<1e-06))>0&s<20)
{
s<-s+1
fit<-lmer(c(coet)~0+coe+(0+coe|Sub.all),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
cor.t<-1-c(attr(VarCorr(fit)[[1]],"correlation")-diag(diag(attr(VarCorr(fit)[[1]],"correlation"))))
}
}
b0<-c(summary(fit)$coefficients[1,1],summary(fit)$coefficients[2,1],summary(fit)$coefficients[3,1])
Phi<-diag(c(as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[1]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[2]),
as.numeric((attr(VarCorr(fit)[[1]],"stddev")^2)[3])))
if(random.indep==TRUE)
{
cor.AB=cor.AC=cor.BC<-0
}else
{
cor.AB<-attr(VarCorr(fit)[[1]],"correlation")[1,2]
cor.AC<-attr(VarCorr(fit)[[1]],"correlation")[1,3]
cor.BC<-attr(VarCorr(fit)[[1]],"correlation")[2,3]
}
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
Lambda<-diag(rep(attr(VarCorr(fit),"sc")^2,3))
LL<-as.numeric(logLik(fit,REML=FALSE))
if(is.matrix(var.constraint))
{
if(nrow(var.constraint)>=4)
{
Phi.confint<-var.constraint[1:3,]
Lambda.confint<-matrix(rep(var.constraint[4,],3),nrow=3,ncol=2,byrow=TRUE)
colnames(Phi.confint)=colnames(Lambda.confint)<-c("LB","UB")
rownames(Phi.confint)=rownames(Lambda.confint)<-c("A","B","C")
}else
{
warning("The number of intervals is not correct. The constraint intervals will be estimated instead.")
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-(re.confint[1,])^2
Phi.confint[2,]<-(re.confint[4,])^2
Phi.confint[3,]<-(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-(re.confint[7,])^2
}
}else
if(var.constraint==TRUE)
{
Lambda.confint<-matrix(NA,3,2)
colnames(Lambda.confint)<-c("LB","UB")
rownames(Lambda.confint)<-c("A","B","C")
Phi.confint<-Lambda.confint
re.confint<-confint(fit)
Phi.confint[1,]<-(re.confint[1,])^2
Phi.confint[2,]<-(re.confint[4,])^2
Phi.confint[3,]<-(re.confint[6,])^2
Lambda.confint[1,]=Lambda.confint[2,]=Lambda.confint[3,]<-(re.confint[7,])^2
}
}
}
diff<-100
s<-0
while(s<max.itr&diff>=tol)
{
if(random.indep&!random.var.equal&u.int)
{
V.alpha<-matrix(Phi[1,1],K,K)+diag(rep(Lambda[1,1],K))
V.beta<-matrix(Phi[2,2],K,K)+diag(rep(Lambda[2,2],K))
V.gamma<-matrix(Phi[3,3],K,K)+diag(rep(Lambda[3,3],K))
for(i in 1:N)
{
# A
z1z=z1m=z1z.d=z1m.d<-NULL
for(k in 1:K)
{
dd<-dat[[i]][[k]]
inv.O1<-diag(rep(1/sigma1[i,k],nrow(dd)))
inv.O2<-diag(rep(1/sigma2[i,k],nrow(dd)))
O12<-diag(rep(sqrt(sigma2[i,k]/sigma1[i,k]),nrow(dd)))
z1z<-c(z1z,(t(dd$Z)%*%inv.O1%*%(dd$Z))[1,1])
z1m<-c(z1m,(t(dd$Z)%*%inv.O1%*%(dd$M))[1,1])
z1z.d<-c(z1z.d,(delta^2*t(dd$Z)%*%O12%*%inv.O2%*%O12%*%(dd$Z)/(1-delta^2))[1,1])
z1m.d<-c(z1m.d,(delta*t(dd$Z)%*%O12%*%inv.O2%*%(dd$R-dd$Z*Ct[i,k]-dd$M*Bt[i,k]-delta*O12%*%dd$M)/(1-delta^2))[1,1])
}
At[i,]<-solve(diag(z1z+z1z.d)+solve(V.alpha))%*%(z1m-z1m.d+solve(V.alpha)%*%rep(b0[1],K))
# B
m2m=m2r.d<-NULL
for(k in 1:K)
{
dd<-dat[[i]][[k]]
inv.O1<-diag(rep(1/sigma1[i,k],nrow(dd)))
inv.O2<-diag(rep(1/sigma2[i,k],nrow(dd)))
O12<-diag(rep(sqrt(sigma2[i,k]/sigma1[i,k]),nrow(dd)))
m2m<-c(m2m,(t(dd$M)%*%inv.O2%*%(dd$M)/(1-delta^2))[1,1])
m2r.d<-c(m2r.d,(t(dd$M)%*%inv.O2%*%(dd$R-dd$Z*Ct[i,k]-delta*O12%*%(dd$M-dd$Z*At[i,k]))/(1-delta^2))[1,1])
}
Bt[i,]<-solve(diag(m2m)+solve(V.beta))%*%(m2r.d+solve(V.beta)%*%rep(b0[2],K))
# C
z2z=z2r.d<-NULL
for(k in 1:K)
{
dd<-dat[[i]][[k]]
inv.O1<-diag(rep(1/sigma1[i,k],nrow(dd)))
inv.O2<-diag(rep(1/sigma2[i,k],nrow(dd)))
O12<-diag(rep(sqrt(sigma2[i,k]/sigma1[i,k]),nrow(dd)))
z2z<-c(z2z,(t(dd$Z)%*%inv.O2%*%(dd$Z)/(1-delta^2))[1,1])
z2r.d<-c(z2r.d,(t(dd$Z)%*%inv.O2%*%(dd$R-dd$M*Bt[i,k]-delta*O12%*%(dd$M-dd$Z*At[i,k]))/(1-delta^2))[1,1])
}
Ct[i,]<-solve(diag(z2z)+solve(V.gamma))%*%(z2r.d+solve(V.gamma)%*%rep(b0[3],K))
}
# update sigma_1ik and sigma_2ik
if(Sigma.update)
{
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
Y<-cbind(dd$M,dd$R)
X<-cbind(dd$Z,dd$M)
Theta<-matrix(c(At[i,k],0,Ct[i,k],Bt[i,k]),2,2)
S<-t(Y-X%*%Theta)%*%(Y-X%*%Theta)
sigma1[i,k]<-(S[1,1]-delta*S[1,2]*sqrt(S[1,1]/S[2,2]))/(nrow(dd)*(1-delta^2))
sigma2[i,k]<-(S[2,2]-delta*S[1,2]*sqrt(S[2,2]/S[1,1]))/(nrow(dd)*(1-delta^2))
}
}
}
vec.At<-c(At)
vec.Bt<-c(Bt)
vec.Ct<-c(Ct)
fit.A<-lme(vec.At~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.B<-lme(vec.Bt~1,random=~1|Sub,control=lmeControl(opt="optim"))
fit.C<-lme(vec.Ct~1,random=~1|Sub,control=lmeControl(opt="optim"))
Afix<-as.numeric(summary(fit.A)$coefficients$fixed[1])
Bfix<-as.numeric(summary(fit.B)$coefficients$fixed[1])
Cfix<-as.numeric(summary(fit.C)$coefficients$fixed[1])
if(random.var.update)
{
Phi<-diag(c(as.numeric(VarCorr(fit.A)[1,1]),as.numeric(VarCorr(fit.B)[1,1]),as.numeric(VarCorr(fit.C)[1,1])))
Lambda<-diag(c(as.numeric(VarCorr(fit.A)[2,1]),as.numeric(VarCorr(fit.B)[2,1]),as.numeric(VarCorr(fit.C)[2,1])))
if(sum(var.constraint==FALSE)==0)
{
# Phi
if(Phi[1,1]<Phi.confint[1,1])
{
Phi[1,1]<-Phi.confint[1,1]
}else
if(Phi[1,1]>Phi.confint[1,2])
{
Phi[1,1]<-Phi.confint[1,2]
}
if(Phi[2,2]<Phi.confint[2,1])
{
Phi[2,2]<-Phi.confint[2,1]
}else
if(Phi[2,2]>Phi.confint[2,2])
{
Phi[2,2]<-Phi.confint[2,2]
}
if(Phi[3,3]<Phi.confint[3,1])
{
Phi[3,3]<-Phi.confint[3,1]
}else
if(Phi[3,3]>Phi.confint[3,2])
{
Phi[3,3]<-Phi.confint[3,2]
}
# Lambda
if(Lambda[1,1]<Lambda.confint[1,1])
{
Lambda[1,1]<-Lambda.confint[1,1]
}else
if(Lambda[1,1]>Lambda.confint[1,2])
{
Lambda[1,1]<-Lambda.confint[1,2]
}
if(Lambda[2,2]<Lambda.confint[2,1])
{
Lambda[2,2]<-Lambda.confint[2,1]
}else
if(Lambda[2,2]>Lambda.confint[2,2])
{
Lambda[2,2]<-Lambda.confint[2,2]
}
if(Lambda[3,3]<Lambda.confint[3,1])
{
Lambda[3,3]<-Lambda.confint[3,1]
}else
if(Lambda[3,3]>Lambda.confint[3,2])
{
Lambda[3,3]<-Lambda.confint[3,2]
}
}
}
u<-cbind(as.matrix(ranef(fit.A)),as.matrix(ranef(fit.B)),as.matrix(ranef(fit.C)))
colnames(u)<-c("A","B","C")
bnew<-c(Afix,Bfix,Cfix)
diff<-max(abs(bnew-b0))
b0<-bnew
s<-s+1
}else
{
if((random.var.equal&u.int)|(!random.indep&u.int))
{
warning("The u.int=TRUE argument will be ignored and the full h-likelihood method will be applied.")
}
bik.mat<-array(NA,c(N,K,3))
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
colnames(dd)<-c("Z","M","R")
n[i,k]<-nrow(dd)
Z<-matrix(dd[,1])
M<-matrix(dd[,2])
R<-matrix(dd[,3])
r<-delta*sqrt(sigma2[i,k]/sigma1[i,k])
X<-cbind(Z,M)
Pik<-matrix(c(-r,0,0,1,1,0),2,3)
Qik<-c(0,r)
Vik<-matrix(c(1,0,0),1,3)
bik<-solve(sigma1[i,k]*t(Pik)%*%t(X)%*%X%*%Pik+sigma2[i,k]*(1-delta^2)*t(Vik)%*%t(Z)%*%Z%*%Vik+
sigma1[i,k]*sigma2[i,k]*(1-delta^2)*solve(Lambda))%*%
(sigma1[i,k]*t(Pik)%*%t(X)%*%(R-X%*%Qik)+sigma2[i,k]*(1-delta^2)*t(Vik)%*%t(Z)%*%M+
sigma1[i,k]*sigma2[i,k]*(1-delta^2)*solve(Lambda)%*%(b0+u[i,]))
At[i,k]<-bik[1,1]
Bt[i,k]<-bik[2,1]
Ct[i,k]<-bik[3,1]
bik.mat[i,k,]<-bik
}
}
e<-c(sum(At-b0[1])/N,sum(Bt-b0[2])/N,sum(Ct-b0[3])/N)
for(i in 1:N)
{
u[i,]<-solve(solve(Phi)+K*solve(Lambda))%*%solve(Lambda)%*%(apply(cbind(At[i,],Bt[i,],Ct[i,]),2,sum)-K*b0-e)
}
bnew<-c(mean(At),mean(Bt),mean(Ct))-apply(u,2,mean)
if(random.var.update==TRUE)
{
Phi<-t(u)%*%u/N
if(det(Phi)<10e-5|random.indep==TRUE)
{
Phi<-diag(as.vector(diag(Phi)))
}
if(random.var.equal==FALSE)
{
Lambda.new<-matrix(0,3,3)
for(i in 1:N)
{
for(k in 1:K)
{
Lambda.new<-Lambda.new+(bik.mat[i,k,]-bnew-u[i,])%*%t(bik.mat[i,k,]-bnew-u[i,])/(N*K)
}
}
Lambda<-diag(as.vector(diag(Lambda.new)))
}else
{
theta2<-0
for(i in 1:N)
{
for(k in 1:K)
{
theta2<-theta2+(t(bik.mat[i,k,]-bnew-u[i,])%*%(bik.mat[i,k,]-bnew-u[i,])/(3*N*K))[1,1]
}
}
Lambda<-diag(rep(theta2,3))
}
Phi<-t(u)%*%u/N
if(random.indep)
{
Phi<-diag(as.vector(diag(Phi)))
}
# add variance constraint
if(sum(var.constraint==FALSE)==0)
{
# Phi
if(Phi[1,1]<Phi.confint[1,1])
{
Phi[1,1]<-Phi.confint[1,1]
}else
if(Phi[1,1]>Phi.confint[1,2])
{
Phi[1,1]<-Phi.confint[1,2]
}
if(Phi[2,2]<Phi.confint[2,1])
{
Phi[2,2]<-Phi.confint[2,1]
}else
if(Phi[2,2]>Phi.confint[2,2])
{
Phi[2,2]<-Phi.confint[2,2]
}
if(Phi[3,3]<Phi.confint[3,1])
{
Phi[3,3]<-Phi.confint[3,1]
}else
if(Phi[3,3]>Phi.confint[3,2])
{
Phi[3,3]<-Phi.confint[3,2]
}
Phi[1,2]=Phi[2,1]<-cor.AB*sqrt(Phi[1,1]*Phi[2,2])
Phi[1,3]=Phi[3,1]<-cor.AC*sqrt(Phi[1,1]*Phi[3,3])
Phi[2,3]=Phi[3,2]<-cor.BC*sqrt(Phi[2,2]*Phi[3,3])
# Lambda
if(Lambda[1,1]<Lambda.confint[1,1])
{
Lambda[1,1]<-Lambda.confint[1,1]
}else
if(Lambda[1,1]>Lambda.confint[1,2])
{
Lambda[1,1]<-Lambda.confint[1,2]
}
if(Lambda[2,2]<Lambda.confint[2,1])
{
Lambda[2,2]<-Lambda.confint[2,1]
}else
if(Lambda[2,2]>Lambda.confint[2,2])
{
Lambda[2,2]<-Lambda.confint[2,2]
}
if(Lambda[3,3]<Lambda.confint[3,1])
{
Lambda[3,3]<-Lambda.confint[3,1]
}else
if(Lambda[3,3]>Lambda.confint[3,2])
{
Lambda[3,3]<-Lambda.confint[3,2]
}
}
}
# update sigma_1ik and sigma_2ik
if(Sigma.update)
{
for(i in 1:N)
{
for(k in 1:K)
{
dd<-dat[[i]][[k]]
Y<-cbind(dd$M,dd$R)
X<-cbind(dd$Z,dd$M)
Theta<-matrix(c(At[i,k],0,Ct[i,k],Bt[i,k]),2,2)
S<-t(Y-X%*%Theta)%*%(Y-X%*%Theta)
sigma1[i,k]<-(S[1,1]-delta*S[1,2]*sqrt(S[1,1]/S[2,2]))/(nrow(dd)*(1-delta^2))
sigma2[i,k]<-(S[2,2]-delta*S[1,2]*sqrt(S[2,2]/S[1,1]))/(nrow(dd)*(1-delta^2))
}
}
}
diff<-max(abs(b0-bnew))
b0<-bnew
s<-s+1
}
}
HL<-cma.h(dat,delta=delta,A.ik=At,B.ik=Bt,C.ik=Ct,b=b0,u=u,Phi=Phi,Lambda=Lambda,random.indep=random.indep,
u.int=u.int,Sigma.update=Sigma.update)$h
zc<-qnorm(1-alpha/2)
coe.re<-matrix(NA,6,4)
colnames(coe.re)<-c("Estimate","LB","UB","SE")
rownames(coe.re)<-c("A","C","B","C2","AB.prod","AB.diff")
# A, C, B, AB.prod, AB.diff
coe.re[1,4]<-sqrt(Phi[1,1]/N+Lambda[1,1]/(N*K))
coe.re[2,4]<-sqrt(Phi[3,3]/N+Lambda[3,3]/(N*K))
coe.re[3,4]<-sqrt(Phi[2,2]/N+Lambda[2,2]/(N*K))
coe.re[1,1:3]<-c(b0[1],b0[1]-zc*coe.re[1,4],b0[1]+zc*coe.re[1,4])
coe.re[2,1:3]<-c(b0[3],b0[3]-zc*coe.re[2,4],b0[3]+zc*coe.re[2,4])
coe.re[3,1:3]<-c(b0[2],b0[2]-zc*coe.re[3,4],b0[2]+zc*coe.re[3,4])
coe.re[5,1]<-b0[1]*b0[2]
coe.re[5,4]<-sqrt((b0[1]*coe.re[3,4])^2+(b0[2]*coe.re[1,4])^2)
coe.re[5,2:3]<-c(coe.re[5,1]-zc*coe.re[5,4],coe.re[5,1]+zc*coe.re[5,4])
# C2, AB.diff
if(optimizer[1]=="optimx")
{
fit.C2<-lmer(c(C2t)~(1|Sub),control= lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B")))
}else
{
fit.C2<-lmer(c(C2t)~(1|Sub),control=lmerControl(optimizer=optimizer[1]))
}
C2fix<-summary(fit.C2)$coefficients[1]
s2.C2<-c(as.numeric(VarCorr(fit.C2)),attr(VarCorr(fit.C2),"sc")^2)
e.C2<-C2fix+mean(ranef(fit.C2)$Sub[,"(Intercept)"])
coe.re[4,4]<-sqrt(s2.C2[1]/N+s2.C2[2]/(N*K))
coe.re[4,1:3]<-c(C2fix,C2fix-zc*coe.re[4,4],C2fix+zc*coe.re[4,4])
coe.re[6,1]<-coe.re[4,1]-coe.re[2,1]
coe.re[6,4]<-sqrt(coe.re[4,4]^2+coe.re[2,4]^2)
coe.re[6,2:3]<-c(coe.re[6,1]-zc*coe.re[6,4],coe.re[6,1]+zc*coe.re[6,4])
s2.C2<-data.frame(C2=s2.C2[1],Error=s2.C2[2],Total=sum(s2.C2))
colnames(s2.C2)<-c("Random(C')","Var(Error)","Var(C')")
var.comp<-data.frame(delta=delta,A=Phi[1,1],C=Phi[3,3],B=Phi[2,2],Lambda[1,1],Lambda[3,3],Lambda[2,2])
colnames(var.comp)<-c("delta","Random(A)","Random(C)","Random(B)","Var(Error A)","Var(Error C)","Var(Error B)")
if(random.indep)
{
Phi<-diag(as.vector(diag(Phi)))
}
cor.comp.temp<-sqrt(diag(1/diag(Phi)))%*%Phi%*%sqrt(diag(1/diag(Phi)))
cor.comp<-cor.comp.temp
cor.comp[1,2]=cor.comp[2,1]<-cor.comp.temp[1,3]
cor.comp[1,3]=cor.comp[3,1]<-cor.comp.temp[1,2]
colnames(cor.comp)=rownames(cor.comp)<-c("A","C","B")
re<-list(delta=delta,Coefficients=coe.re,Cor.comp=cor.comp,Var.comp=var.comp,Var.C2=s2.C2,logLik=LL,HL=HL)
if(sum(var.constraint==FALSE)==0)
{
re.var.constraint<-rbind(Phi.confint,Lambda.confint)
rownames(re.var.constraint)<-c("Random(A)","Random(B)","Random(C)","Var(Error A)","Var(Error B)","Var(Error C)")
re<-list(delta=delta,Coefficients=coe.re,Cor.comp=cor.comp,Var.comp=var.comp,Var.constraint=re.var.constraint,
Var.C2=s2.C2,logLik=LL,HL=HL,convergence=(s<max.itr|max.itr==0))
}
return(re)
}
Try the macc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.