Nothing
R/R_steps.R
In pfda: Paired Functional Data Analysis
{ # global scope variables
# these variables are defined for the sole purpose of confusing the checker
# Each of these should be masked or ignored by a different frame path.
# This is necessary due to the eval/structure programming method employed here.
subject<-
y<-
Z<-
B<-
Bt<-
Bx<-
fname<-
s.xi<-
s.eps<-
ka<-
kb<-
kz<-
nobs<-
M<-
N<-
p<-
K<-
Kt<-
Kx<-
ln<-
lf<-
lg<-
min.v<-
max.I<-
tol<-
name.t<-
name.x<-
tbase<-
xbase<-
NULL
}
{ # general steps
.gen.tm<-function(y,B,subject,tf,alpha,sigma,l,K){
left<-sigma*l*K
right<-matrix(0,NCOL(B),1)
for(i in seq(nlevels(subject))){
ix <- i == as.integer(subject)
Bi<-B[ix,]
left = left + crossprod(Bi)
right = right + crossprod(Bi,y[ix]-Bi%*%tf%*%alpha[i,])
}
solve(left,right)
}
.gen.tf<-function(y,B,subject,tm,tf,alpha,sigma,aa,lf,K){
n<-nlevels(subject)
k<-NCOL(tf)
for(est in seq(k)){
left <- sigma*lf*K
right <- matrix(0,NCOL(B),1)
for(i in seq(n)){
ix <- i == as.integer(subject)
Bi <- B[ix,]
left = left + crossprod(Bi)*(aa[est,est,i])
sne <- matrix(0,NROW(Bi),1)
for(l in setdiff(seq(k),est))sne = sne + Bi%*%tf[,l]*aa[est,l,i]
right = right + crossprod(Bi, (y[ix]-Bi%*%tm)*alpha[i,est] - sne )
}
tf[,est]<-solve(left,right)
}
tf*rep(sign(tf[1,]),each=NROW(tf))
}
.gen.dual.lambda<-function(sum.aa,sum.ab){
t(sum.ab)%*%solve(sum.aa)
}
.gen.orthog<-function(tf,alpha,sum.aa){
n<-NROW(alpha)
k<-NCOL(tf)
Sa<-sum.aa/n
e<-eigen(tf%*%tcrossprod(Sa,tf))
qf = e$vectors[,seq_len(k),drop=FALSE]
qf = qf*rep(sign(qf[1,]),each=NROW(qf))
Da = x=e$values[seq_len(k)]
R = crossprod(qf,tf)
list(tf=qf,D=Da,d=tcrossprod(alpha,R),transformation=R)
}
.gen.dual.sigmas<-function(By,Bz,tf,tg,lambda,Da,Db,s.eps,s.xi){
Seta <- diag(Db,length(Db))-lambda%*%tcrossprod(diag(Da,length(Da)),lambda)
Setai<- solve(Seta)
phi <- By %*% tf
psi <- Bz %*% tg
Sa <- diag(1/Da,length(Da)) + crossprod(lambda,Setai)%*%lambda + crossprod(phi)/s.eps
Sab1 <- -crossprod(lambda,Setai)
Sb <- Setai+crossprod(psi)/s.xi
Saa <- solve(Sa - Sab1%*%solve(Sb,t(Sab1)))
Sab <- -Saa%*%Sab1%*%solve(Sb)
Sbb <- solve(Sb - crossprod(Sab1,solve(Sa,Sab1)))
list(Saa=Saa,Sab=Sab,Sbb=Sbb)
}
.gen.symblock.solve<-function(A,B,C){
#! Solves matrixof form
#! [[ A B]
#! [t(B) C]]
#! Returning the block inverses
#!
stopifnot(NROW(A)==NROW(B))
stopifnot(NCOL(C)==NCOL(B))
svAB = solve(A,B)
Ci = solve(C-crossprod(B,svAB))
Bi = -svAB%*%Ci
Ai = solve(A)+ svAB%*%tcrossprod(Ci,svAB)
list(Ai,Bi,Ci)
}
.X.handle.z<-expression({ # handle Z
stopifnot(exists('Z',inherits=FALSE))
if(missing(Z)||is.null(Z))
Z = matrix(nrow=length(y),ncol=0)
else {
if(is(Z,"formula"))
Z = model.matrix(Z)
else Z = as.matrix(Z)
}
})
.X.subset<-expression({
stopifnot(exists('subset',inherits=FALSE))
if(!missing(subset) && !is.null(subset)){
if(exists('Z',inherits=FALSE))Z<-subset(Z,subset)
y<-subset(y,subset)
if(exists('z',inherits=FALSE))z<-subset(z,subset)
t<-subset(t,subset)
if(exists('x',inherits=FALSE))x<-subset(x,subset)
subject<-subject[subject,drop=T]
}
})
.X.single.knots<-expression({ # knots identification
if(is.list(knots)){
Bt = knots[[1]]
Kt = knots[[2]]
if(!all(dim(Kt)==dim(Bt)[2]))stop("If knots are provided as a list the first component must be the evaluated coefficient matrix for the basis functions, and the second must be the penalty matrix to be used.")
}else {
if(!exists("tbase",inherits=FALSE) || is.null(tbase)){
if(is.null(knots)){
kt<-expand.knots(unique(quantile(t,seq(0,1,length.out=control$nknots))))
tbase = OBasis(kt)
} else if(is(knots,'numeric')) {
kt<-knots
tbase = OBasis(kt)
} else {
tbase = knots
}
Bt = evaluate(tbase,t)
Kt = penaltyMatrix(tbase)
}
}
})
.X.single.penalties<-expression({ # penalties
if(is.null(penalties)){
penalties<- if(is.null(df)) rep(NA,2) else c(l.from.df(df[1],Bt,Kt),l.from.df(df[2],Bt,Kt))
}
})
.X.single.optimize.npc<-expression({
# funcall$k=1
# if(exists('y'))funcall$y=y
# if(exists('Z'))funcall$y=Z
# if(exists('t'))funcall$y=t
# if(exists('subject'))funcall$y=subject
# model.k0<-eval(funcall,env=attr(funcall,'envir'))
message("optimizing number of principal components using AIC")
k=1
model.k0 <- RecallWith(k=k)
while(TRUE){
# funcall$k<-funcall$k+1
# model.k1<-try(eval(funcall,env=attr(funcall,'envir')),silent=TRUE)
k=k+1
model.k1 <- RecallWith(k=k)
if(class(model.k1)[1]=="try-error")return(model.k0)
else if(AIC(model.k0)<AIC(model.k1))return(model.k0)
else model.k0<-model.k1
}
})
.F.single.optimize.npc<-function(){
message("optimizing number of principal components using AIC")
# localfuncs("RecallWith")
k=1
model.k0 <- tryCatch(RecallWith(k=k,fname=fname),warning=function(w)stop(w$message))
tryCatch(while(TRUE){
k=k+1
model.k1 <- RecallWith(k=k,fname=fname)
if(AIC(model.k0)<AIC(model.k1))return(model.k0)
else model.k0<-model.k1
}, warning=function(w){
if(w$message=="EM-algorithm did not converge")model.k0 else stop(w$message)
}, error=function(e)model.k0)
}
.X.dual.k<-expression({ # resloves k input for dual pc (cc/bc/add)
if(is.null(k)) k<-rep(NA,2)
if(length(k)!=2) k<-rep(as.vector(k),length.out=2)
ka = as.integer(k[1])
kb = as.integer(k[2])
# if(names(k)==NULL) names(k)<-c(name.t,name.x)
})
.X.dual.penalties <-expression({ # resolve penalties input
if(is.null(penalties)){
penalties<- if(exists('Bx',inherits=FALSE))
if(is.null(df)) matrix(NA,2,2) else c(l.from.df(df[1],Bt,Kt),l.from.df(df[2],Bx,Kx),l.from.df(df[3],Bt,Kt),l.from.df(df[4],Bx,Kx))
else
if(is.null(df)) matrix(NA,2,2) else c(l.from.df(df[1],Bt,Kt),l.from.df(df[2],Bt,Kt),l.from.df(df[3],Bt,Kt),l.from.df(df[4],Bt,Kt))
}
if(!class(penalties)=='matrix') penalties<-matrix(penalties,2,2)
# if(is.null(colnames(penalties))) colnames(penalties)<-c(name.t,name.x)
})
.X.read.parameters <- expression({
if(exists('Z',inherits=FALSE)){ kz = if(is.null(Z))0L else NCOL(Z)}
k = as.integer(k)
N = nlevels(subject)
nobs= table(subject)
M = length(y)
p = ncol(Bt)
})
.X.binary.y<-expression({ # enforces the conditions on a binary y
y<-if(is(y,'factor')){
stopifnot(nlevels(y)==2)
as.integer(y)-1
} else if(is(y,'logical')){
as.integer(y)
} else {
stopifnot(identical(as.integer(sort(unique(y))), as.integer(0:1)))
as.integer(y)
}
})
.X.funcall.replacevars<-expression({
funcall$y=y
if(exists("z",inherits = FALSE))funcall$z=z
if(exists("Z",inherits = FALSE))funcall$Z=Z
funcall$t = t
if(exists("x",inherits = FALSE))funcall$x = x
funcall$subject = subject
funcall$k=k
funcall$penalties = penalties
funcall$df = df
if(exists("bases",inherits=FALSE))funcall$bases=bases
funcall$knots = knots
funcall$control = control
})
.X.optimize.penalties<-expression({
eval(.X.funcall.replacevars)
pix<-which(is.na(penalties))
if(control$penalty.method=='CV'){
message("optimizing penalties using cross-validation")
if(is.null(control$folds))control$folds<-cv.folds(nlevels(subject),10)
ews<-function(s,p){
ix = !(subject %in% s)
fc<-funcall
fc$penalties=p
fc$subset=ix
model<-eval(fc,env=attr(fc,'envir'))
logLik(model,newdata=list(y=y[!ix],Z=Z[!ix,,drop=FALSE],Bt=Bt[!ix,,drop=FALSE],Bx=Bx[!ix,,drop=FALSE],subject=subject[!ix,drop=TRUE]),n2L=TRUE)
}
cvf<-function(pen,...){
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
cvl<-try(lapply(control$folds,ews,p=p),TRUE)
if(class(cvl)=="try-error")return(NA)
(sum(unlist(cvl)))
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(1,length(pix))
optimpar<-optim(log(control$optim.start),cvf,method=control$optim.method)
penalties[pix]<-exp(optimpar$par)
funcall$penalties=penalties
eval(funcall,env=attr(funcall,'envir'))
} else
if(control$penalty.method=='AIC') {
message("optimizing penalties using AIC")
aicf<-function(pen){
fc<-funcall
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
fc$penalties<-p
m<-try(eval(fc,env=attr(fc,'envir')),silent=TRUE)
if(class(m)[1]=="try-error") NA else AIC(m)
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(l.from.df(2.1,Bt,Kt),length(pix))
optimpar<-optim(log(control$optim.start),aicf,method=control$optim.method)
penalties[pix] <- exp(optimpar$par)
funcall$penalties=penalties
eval(funcall,env=attr(funcall,'envir'))
}
})
.F.optimize.penalties<-function(){
# eval(.X.funcall.replacevars)
pix<-which(is.na(penalties))
if(control$penalty.method=='CV'){
message("optimizing penalties using cross-validation")
if(is.null(control$folds))control$folds<-cv.folds(nlevels(subject),10)
ews<-function(s,p){
ix = !(subject %in% s)
fc<-funcall
fc$penalties=p
fc$subset=ix
model<-eval(fc,env=attr(fc,'envir'))
logLik(model,newdata=list(y=y[!ix],Z=Z[!ix,,drop=FALSE],Bt=Bt[!ix,,drop=FALSE],Bx=Bx[!ix,,drop=FALSE],subject=subject[!ix,drop=TRUE]),n2L=TRUE)
}
cvf<-function(pen,...){
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
cvl<-try(lapply(control$folds,ews,p=p),TRUE)
if(class(cvl)=="try-error")return(NA)
(sum(unlist(cvl)))
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(1,length(pix))
optimpar<-optim(log(control$optim.start),cvf,method=control$optim.method)
penalties[pix]<-exp(optimpar$par)
funcall$penalties=penalties
eval(funcall,env=attr(funcall,'envir'))
} else
if(control$penalty.method=='AIC') {
message("optimizing penalties using AIC")
aicf<-function(pen){
# fc<-funcall
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
# fc$penalties<-p
tryCatch({
m<-RecallWith(penalties=p,fname=fname)
if(class(m)[1]=="try-error") NA else AIC(m)
}
,error=function(e)NA
,warning=function(w){
if(w$message=="EM-algorithm did not converge") return(NA)
else stop(w$message)
}
)
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(l.from.df(2.1,Bt,Kt),length(pix))
optimpar<-optim(log(control$optim.start),aicf,method=control$optim.method)
penalties[pix] <- exp(optimpar$par)
# funcall$penalties=penalties
# eval(funcall,env=attr(funcall,'envir'))
RecallWith(penalties=penalties,fname=fname)
}
}
}
{ # utilities and convenience functions
.u.single.resid<-function(y,B,subject,tm,tf,alpha){
for(i in seq_len(nlevels(subject))){
ix = i == as.integer(subject)
y[ix]<-y[ix]-B[ix,,drop=FALSE]%*%(tm+tf%*%alpha[i,])
}
y
}
.u.orthogonalize<-function(U,v){
# U is assumed to be orthogonal already
u = as.vector(if(ncol(U)>0) v-U%*%as.vector(crossprod(v,U)) else v)
u/sqrt(crossprod(u))
}
cv.folds<-function(n, folds = 10){
split(sample(1:n), rep(1:folds, length = n))
}
positive.first.row<-function(X){
stopifnot(is.matrix(X))
X*rep(sign(X[1,]),each=nrow(X))
}
.pfda.df<-Vectorize(function(B,l,K,sigma=1){
tr<-function(x)sum(diag(x))
while(TRUE){
L <- crossprod(B)+sigma*l*K
if(kappa(L)<1e14)break
else if(l>1) l <- l*.90 else NA #l <- l*1.1
}
tr(solve(L,crossprod(B)))
},'l')
l.from.df<-function(df,B,K,sigma=1)if(is.na(df)) NA else uniroot(function(l).pfda.df(B,l,K,sigma)-df,c(1e-4,1e10))$root
RecallWith<-function(...,fname){
calls <- unlist(lapply(sys.calls(),function(x)deparse(x[[1]])[1]))
if(missing(fname)){
iRW<-max(which(calls=="RecallWith"))
if(iRW==1) stop("RecallWith must be called within another function")
unacceptable.calls<-c(which(calls=="{"),unlist(sapply(c("^eval","^function","^try","TryCatch","do.in.envir","^debug:::","^get","assign","mlocal"),grep,calls)))
acceptable.calls<- setdiff(seq_len(iRW-1),unacceptable.calls)
w = max(acceptable.calls)
fname<-calls[w]# deparse(match.call(call = sys.call(w))[[1]])
} else {
w = max(which(calls==fname))
}
oldArgs<-formals(sys.function(w))
if(is.null(oldArgs))oldArgs<-list()
for(a in names(oldArgs)){
try(oldArgs[[a]]<-get(a,sys.frame(w)),silent=TRUE)
}
do.call(fname,as.list(modifyList(oldArgs,list(...))))
}
}
{ # single version
.single.c.i<-function(y,B,subject,k,min.v){
tm = solve(crossprod(B)+diag(min.v,NCOL(B)),crossprod(B,y))
Ry = y-B%*%tm
N = nlevels(subject)
tfa<-matrix(0,NCOL(B),N)
T=matrix(0,NCOL(B),NCOL(B))
for(i in seq_len(nlevels(subject))){
ix <- i==as.integer(subject)
tfa[,i] = solve(crossprod(B[ix,])+min.v*diag(NCOL(B)),crossprod(B[ix,],Ry[ix]))
Ry[ix]<-Ry[ix]-B[ix,]%*%tfa[,i]
T<-T+tcrossprod(tfa[,i])
}
e<-eigen(T)
tf<-e$vectors[,seq_len(k),drop=FALSE]
tf<-tf*rep(sign(tf[1,]),each=NROW(tf))
alpha<-matrix(0,N,k)
aa<-list()
for(i in seq_len(nlevels(subject))){
alpha[i,]<-solve(crossprod(tf),crossprod(tf,tfa[,i]))
aa[[i]] <- tcrossprod(alpha[i,])
}
sigma=crossprod(.u.single.resid(y,B,subject,tm,tf,alpha))/length(y)
return(list(tm=tm,tf=tf,alpha=alpha,aa=aa,Da=e$values[seq_len(k)],sigma=sigma))
}
.single.c.1.1<-function(yi,Bi,tm,tf,Da,sigma){
Sa = solve(crossprod(Bi%*%tf)+diag(1/Da,length(Da)))
mu = Sa%*%crossprod(Bi%*%tf,yi-Bi%*%tm)/sigma
aa = tcrossprod(mu)+Sa
list(alpha=mu,Sa=Sa,aa=aa)
}
.single.c.1<-function(y,B,subject,tm,tf,Da,sigma){
alpha<-matrix(nrow=nlevels(subject),ncol=NCOL(tf))
aa<-array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
Sa<-array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
r<-.single.c.1.1(y[ix],B[ix,],tm,tf,Da,sigma)
alpha[i,]<-r$alpha
Sa[,,i]<-r$Sa
aa[,,i]<-r$aa
}
list(alpha=alpha,Saa=Sa,aa=aa)
}
.single.c.2<-function(y,B,subject,tm,tf,alpha,Sa){
Ry <- y-B%*%tm
s<-0
for(i in seq_len(nlevels(subject))){
ix<- i==as.integer(subject)
phi=B[ix,,drop=FALSE]%*%tf
s<-s+crossprod(Ry[ix]-phi%*%alpha[i,])+sum(diag(phi%*%Sa[,,i]%*%t(phi)))
}
list(sigma=as.vector(s/length(subject)))
}
.single.c.3<-function(y,B,subject,tf,alpha,sigma,lm,K){
list(tm=.gen.tm(y,B,subject,tf,alpha,sigma,lm,K))
}
.single.c.4<-function(y,B,subject,tm,tf,alpha,sigma,aa,lf,K){
list(tf=.gen.tf(y,B,subject,tf,tm,alpha,sigma,aa,lf,K))
}
.single.c.5<-function(tf,alpha,aa){
saa<-matrix(0,NCOL(tf),NCOL(tf))
for(sub in seq(NROW(alpha))) saa = saa + aa[[sub]]
a<-.gen.orthog(tf,alpha,saa)
list(tf=a$tf,Da=a$D, alpha = a$d)
}
single.c.core<-function(y,B,subject,k,lm,lf,K,min.v,max.I,tol){
{ #Initial values
r<-.single.c.i(y,B,subject,k,min.v)
tm=r$tm
tf=r$tf
alpha=r$alpha
aa=r$aa
Da=r$Da
sigma=r$sigma
I=0
.cc<-numeric(0)
}
while(I<-I+1){
s1<-.single.c.1(y,B,subject,tm,tf,Da,sigma)
s2<-.single.c.2(y,B,subject,tm,tf,s1$alpha,s1$Saa)
s3<-.single.c.3(y,B,subject,tf,s1$alpha,s2$sigma,lm,K)
s4<-.single.c.4(y,B,subject,s3$tm,tf,s1$alpha,s2$sigma,s1$aa,lf,K)
s5<-.single.c.5(s4$tf,s1$alpha,s1$aa)
{ # convergence checks
ccl<-numeric(0)
ccl['tm'] <- sum(abs((tm-s3$tm)/s3$tm))
ccl['tf'] <- sum(abs((tf-s5$tf)/s5$tf))
ccl['Da'] <- sum(abs((Da-s5$Da)/s5$Da))
cc <- sum(ccl)
.cc<-c(.cc,cc)
if(cc<tol)break
if(I>=max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
{ # reassign values to
tm<-s3$tm
tf<-s5$tf
alpha<-s5$alpha
Da<-s5$Da
sigma<-s2$sigma
aa<-s5$aa
Saa<-s1$Saa
}
}
}
list(y=y,Bt=B,subject=subject,tm=tm,tf=tf,alpha=alpha, Da=Da,sigma=sigma,aa=aa,Saa=Saa,lm=lm,lf=lf,K=K,k=k,I=I,cc=cc)
}
.single.c.n2L<-function(y, subject, B, tm, tf, Da, sigma){
phi<-B%*%tf
L <- 0.0
for(i in seq_len(nlevels(subject))){
ix<-i==as(subject,"integer")
Sy<-phi[ix,,drop=FALSE]%*%tcrossprod(diag(x=Da,ncol=length(Da)),phi[ix,,drop=FALSE])+diag(sigma,nrow=sum(ix))
Ry<-y[ix]-B[ix,]%*%tm
D<-determinant(Sy,logarithm=TRUE)
if(D$sign<0)stop("Negative determinant encountered for variance matrix.")
L<-L+as.numeric(D$modulus+crossprod(Ry,solve(Sy,Ry)))
}
L
}
AIC.pfda.single.c<-function(object,...){
with(object,{
as.vector(.single.c.n2L(y, subject, Bt, tm, tf, Da, sigma))+
2*(.pfda.df(Bt,lm,K,sigma)+k*.pfda.df(Bt,lf,K,sigma))
})
}
# AIC.pfda.single.c.R<-function(object,...)
logLik.pfda.single.c.R<-logLik.pfda.single.c.rawC<-function(object,...,newdata=NULL, n2L=TRUE){
r<-with(object,with(newdata,.single.c.n2L(y,subject,B,tm,tf,Da,sigma)))
if(n2L) r else exp(-r/2)
}
single.c<-function(y,Z,t,subject,knots=NULL,penalties=NULL,df=NULL,k=NULL,control=pfdaControl(),subset=NULL,...){
{ # setup
fname = deparse(match.call()[[1L]])
localfuncs(c('.F.single.optimize.npc','.F.optimize.penalties'))
name.t = deparse(substitute(t))
eval(.X.handle.z)
eval(.X.subset)
eval(.X.single.knots)
eval(.X.single.penalties)
}
if(is.null(k)||any(is.na(k))){
# funcall <- match.call()
#stop("identification of number of principle components is not done yet.")
# eval(.X.single.optimize.npc)
.F.single.optimize.npc()
} else
if (any(is.na(penalties))) {
funcall <- match.call()
.F.optimize.penalties()
}
else {
rtn<-if(control$useC){
{ # setup for passing to Compiled code
eval(.X.read.parameters)
{ #Compute Memory Requirements
pfda_computeResid = M*k
pfda_s_i = p*p + M + N*p + p*p + max(
pfda_eigens =8*p,
pfda_computeResid )
single_c_resid = pfda_computeResid
pfda_m1 = M + M*k + 2*k^2
pfda_m2 = M + M*k + p^2
pfda_m3_for_subjectnum = p
pfda_m3_for_estcol = p^2 + pfda_m3_for_subjectnum
pfda_m3_core = M + pfda_m3_for_estcol
pfda_m3 = k^2 + pfda_m3_core
pfda_m5_0 = p^2 + p^2 + p + k + max(
pfda_eigens = 8*p,
pfda_matrix_outer_quadratic_form =k*p)
pfda_m5_1 = k^2 + pfda_m5_0
pfdaSingle_m5 = k^2 + pfda_m5_1
pfdaSingle_e_1 = k + 2*k^2
pfdaSingle_e = M + M*k + k*N + max(
pfda_computeResid ,
pfdaSingle_e_1 )
single_c_E = M + pfdaSingle_e
single_c_unpenalized = M + kz*kz + max(
single_c_resid ,
dsysv_=10*kz)
single_c_penalized = M + pfda_m2
single_c_princcomp = M + p*p*N + pfda_m3_core
single_c_variances = M + max(
pfda_m1 ,
pfdaSingle_m5 )
single_c_core = p*p*N + pfda_s_i + p + p*k + k + max(
single_c_E ,
single_c_unpenalized ,
single_c_penalized ,
single_c_princcomp ,
single_c_variances ,
single_c_resid )
ipl<-8*p
dpl<- single_c_core
}
}
structure(.C('single_c_core', residuals=y, Z=Z, Bt=Bt, tz=double(kz), tm=double(p), tf=matrix(0,p,k), alpha=matrix(0,N,k), Da=double(k), sigma=0, aa=array(0,dim=c(k,k,N)), Saa=array(0,dim=c(k,k,N)), nobs=nobs, N=N, M=M, kz=kz, k=k, p=p, lm=penalties[1], lf=penalties[2], K=Kt, minV=control$minimum.variance, Iterations=control$max.iterations, tol=control$convergence.tolerance, dl=control$C.debug, dp=double(dpl), ip=integer(max(6*p,kz)))
,class=c('pfda.single.c.rawC','pfda.single.c','list'))
} else {
structure(single.c.core(y,Bt,subject,k,lm=penalties[1],lf=penalties[2],Kt,control$minimum.variance,control$max.iterations,control$convergence.tolerance)
,class=c('pfda.single.c.R','pfda.single.c','list'))
}
rtn$tbase<-tbase
rtn$y<-y
rtn$subject<-subject
return(rtn)
}
}
print.pfda.single.c<-function(x,...){
cat('Univariate Functional Principal Component Model\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(NCOL(x$tf),' principal components\n')
cat('penalties are \n');print( penalty.pfda.single.c(x))
}
penalty.pfda.single.c<-function(x){
with(x,c(mean=lm,pc=lf))
}
plot.pfda.single.c<-function(x,...){
with(x,{
layout(matrix(1:2,nrow=2,ncol=1,byrow=T))
plot(tbase,tm, main=paste("Plot of mean curve for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab=attr(x,'name.y'))
plot(tbase,tf, main=paste("Principle components for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab='')
})
}
}
{ # single binaries
.roberts1<-function(c){
a<-(c+sqrt(c^2+4))/2
while(TRUE){
u=runif(1)
z=c-log(u)/a
p=exp(-(z-a)^2/2)
v=runif(1)
if(v<=p)return(z)
}
}
.reject1<-function(c){
while(TRUE){
z=rnorm(1)
if(z>c)return(z)
}
}
.rtruncnormlower<-function(m,mean,sd,lower){
c = (lower-mean)/sd
if(c>0){
replicate(m,.roberts1(c))
} else {
replicate(m,.reject1(c))
}
}
.single.b.i<-function(y,B,subject,k,min.v){
.single.c.i(y,B,subject,k,min.v)[1:5]
}
.single.b.1a.Sa.i<-function(B,tf,Da)solve(crossprod(B%*%tf)+solve(Da))
.single.b.1a.mu.i<-function(B,w, tm, tf, Sa)Sa%*%crossprod(B%*%tf, w-B%*%tm)
.single.b.1b.aa.i<-function(B,w,ww,tm,tf,Sa){
phi=B%*%tf
s1=-Sa%*%crossprod(phi,B%*%tm)
s2=Sa%*%crossprod(phi,tcrossprod(as.matrix(w),s1))
Sa%*%t(phi)%*%ww%*%phi%*%t(Sa) +
s2 + t(s2) +
tcrossprod(s1) +
Sa
}
.single.b.1.i<-function(B,w,ww,tm,tf,Da){
Sa = .single.b.1a.Sa.i(B,tf,Da)
alpha = .single.b.1a.mu.i(B, w, tm, tf, Sa)
aa = .single.b.1b.aa.i(B,w,ww,tm,tf,Sa)
list(alpha=alpha,aa=aa,Sa=Sa)
}
.single.b.1<-function(B,subject,w,ww,tm,tf,Da){
subject<-as.integer(as.factor(subject))
n<-max(subject)
k<-NCOL(tf)
aa<-vector('list',n)
alpha<-matrix(nrow=n,ncol=ncol(tf))
Sa<-array(0,dim=c(k,k,n))
for(i in seq(n)){
ix<-i == as.integer(subject)
r<-.single.b.1.i(B[ix,],w[ix],ww[[i]],tm,tf,Da)
Sa[,,i]<-r$Sa
alpha[i,]<-r$alpha
aa[[i]]<-r$aa
}
list(alpha=alpha,aa=structure(unlist(aa),dim=c(1,1,10)),Sa=Sa)
}
.single.b.2<-function(B,subject,w,tf,alpha,lm,K){
list(tm=.gen.tm(w,B,subject,tf,alpha,1,lm,K))
}
.single.b.3<-function(B,subject,w,tm,tf,alpha,aa,lf,K){
list(tf=.gen.tf(w,B,subject,tm,tf,alpha,1,aa,lf,K))
}
.single.b.4<-function(tf,alpha,aa){
saa<-matrix(0,NCOL(tf),NCOL(tf))
for(sub in seq(NROW(alpha))) saa = saa + aa[[sub]]
a<-.gen.orthog(tf,alpha,saa)
list(tf=a$tf,Da=a$D, alpha = a$d)
}
.single.b.w.genw<-function(yi,wi,Bi,tm,tf,Da,kr,j){
# R version of pfda_bin_s_gen_w in C code
Bij <- Bi[-j,,drop=F]
wij <- wi[-j]
# pfda_bin_single_generate_w_parms1
Ss <- solve(crossprod(Bij%*%tf)+solve(diag(Da,length(Da))))
mu <- Ss%*%t(Bij%*%tf)%*%(wij-Bij%*%tm)
# pfda_bin_single_generate_w_parms2
a <- t(Bi[j,])%*%tm + t(Bi[j,])%*%tf%*%mu
s <- 1+t(Bi[j,])%*%tf%*%Ss%*%t(tf)%*%Bi[j,]
a+if(yi[j]) s*.rtruncnormlower(kr,0,1,-a/s) else -s*.rtruncnormlower(kr,0,1,a/s)
}
.single.b.w.1<-function(yi,Bi,wi,tm,tf,Da,kr){
# R version of pfda_bin_single_approximate_moments_forobs
n<-NROW(Bi)
w_sim<-matrix(sapply(seq_len(n),.single.b.w.genw,kr=kr,yi=yi,Bi=Bi,wi=wi,tm=tm,tf=tf,Da=Da),nrow=kr)
wi<-apply(w_sim,2,mean)
# wwa<-array(apply(w_sim,2,tcrossprod),dim=c(n,n,kr))
# wws<-matrix(0,n,n)
# for(l in seq_len(kr))wws<-wws+wwa[,,l]
wwi<-crossprod(w_sim)/kr
return(list(wi=wi,wwi=wwi))
}
.single.b.w<-function(y,B,subject,w,ww,tm,tf,Da,weight,kr){
for(i in seq_len(nlevels(subject))){
ix <- i==as.integer(subject)
rtn <- .single.b.w.1(y[ix],B[ix,,drop=F],w[ix],tm,tf,Da,kr)
w[ix] <- (1-weight)*w[ix]+weight*rtn[[1]]
ww[[i]] <- (1-weight)*ww[[i]]+weight*rtn[[2]]
}
return(list(w=w,ww=ww))
}
.single.b.core<-function(y,B,subject,k,lm,lf,K,minimum.variance, max.iterations,convergence.tolerance){
{ # setup variables & initial values
M<-NROW(y)
N<-nlevels(subject)
nobs<-table(subject)
w <- y
ww<-vector('list',N)
for(i in seq_len(N)) ww[[i]]<-matrix(0,nobs[i],nobs[i])
Saa<-array(0,dim=c(k,k,N))
rtn<-.single.b.i(y,B,subject,k,minimum.variance)
tm<-rtn$tm
tf<-rtn$tf
Da<-rtn$Da
alpha<-rtn$alpha
old<-vector('list',0)
aa <- rtn$aa #vector('list',N); for(i in seq_len(N)) aa<-matrix(0,k,k)
I=0
.cc<-numeric(0)
}
while(I<-I+1){
{ #step w
r0=k0=100
kr=10
if(I<r0){
weight=1
rtn<-.single.b.w(y,B,subject,w,ww,tm,tf,Da,k0,weight)
} else {
weight=10/(10+I)
rtn<-.single.b.w(y,B,subject,w,ww,tm,tf,Da,kr,weight)
}
w<-rtn$w
ww<-rtn$ww
}
s1 <- .single.b.1(B,subject,w,ww,tm,tf,Da)
s2 <- .single.b.2(B,subject,w,tf,s1$alpha,lm,K)
s3 <- .single.b.3(B,subject,w,s2$tm,tf,s1$alpha,s1$aa,lf,K)
s4 <- .single.b.4(tf,alpha,aa)
{ #check convergence
ccl<-numeric(0)
ccl['tm'] <- sum(abs((tm-s2$tm)/s2$tm))
ccl['tf'] <- sum(abs((tf-s4$tf)/s4$tf))
ccl['Da'] <- sum(abs((Da-s4$Da)/s4$Da))
cc <- sum(ccl)
if(cc<convergence.tolerance)break
if(I>=max.iterations){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
.cc<-c(.cc,cc)
{ #reassign variable
tm <- s2$tm
tf <- s4$tf
Da <- s4$Da
alpha <- s4$alpha
aa <- s1$aa
Saa <- s1$Sa
}
}
}
list(tm=tm,tf=tf,Da=Da,aa=aa,Saa=Saa,I=I,cc=cc)
}
.single.b.updatePCS<-function(Yi, Da, phii, Y.rhoi, alpha){
ni = length(Yi)
muy = Y.rhoi+phii%*%alpha
mua = pnorm(muy)
Y.tilda = muy + (Yi-mua)/dnorm(mua)
A = tcrossprod(diag(Da,length(Da)), phii)
Q = solve(diag(1, ni) + phii%*%A)
list(alpha= as.vector((A %*% Q) %*% Y.tilda))
}
.single.b.estimatePCS<-function(Y, B, subject, tm, tf, Da){
Y.rho = B%*%tm
phi = B%*%tf
ofun<-function(i){
ix = as.integer(subject)==i
alpha = rep(0,NCOL(tf))
dif = 1
while(dif>1e-3){
ab<-.single.b.updatePCS(Y[ix], Da, phi[ix,,drop=F], Y.rho[ix,,drop=F], alpha)
dif<-mean(abs(alpha-ab$alpha))
alpha = ab$alpha
}
alpha
}
alpha<-matrix(unlist(lapply(seq_len(nlevels(subject)), ofun)),nlevels(subject),NCOL(tf),byrow=T)
alpha
}
.single.b.n2L.1<-function(yi,Bi,tm,tf,alpha,Da){
sum(.single.b.n2L.parts(yi,Bi,tm,tf,alpha,Da))
}
.single.b.n2L.parts<-function(yi,Bi,tm,tf,alpha,Da){
phi = Bi%*%tf
mui = pnorm(Bi%*%tm+phi%*%alpha)
nui = mui*(1-mui)
Wi = as.vector(dnorm(mui)^2/nui)
c(sum(log(Da)),
determinant(diag(NROW(Da))+crossprod(phi,diag(Wi)%*%phi%*%diag(Da,length(Da))))$modulus,
crossprod(alpha,solve(diag(Da,length(Da)))%*%alpha),
-2*sum(yi*log(mui)+ (1-yi)*log(1-mui)-Wi-mui)
)
}
.single.b.n2L.bySubject<-function(y, subject, B, tm, tf, Da, alpha){
sapply(seq_len(nlevels(subject)),function(si){
ix=as.integer(subject)==si
.single.b.n2L.1(y[ix], B[ix,], tm, tf, alpha[si,], Da)
})
}
.single.b.n2L<-function(y, subject, B, tm, tf, Da){
alpha<-.single.b.estimatePCS(y, B, subject, tm, tf, Da)
ll<-sum(sapply(seq_len(nlevels(subject)),function(si){
ix=as.integer(subject)==si
.single.b.n2L.1(y[ix], B[ix,], tm, tf, alpha[si,], Da)
}))
}
loglik.pfda.single.b<-function(object,...,newdata=NULL,n2L=TRUE){
y<-object@FittingData[[1]][[1]]
t<-object@FittingData[[2]]
subject<-object@FittingData[[3]]
B<-evaluate(object@Basis,t)
}
AIC.pfda.single.b<-function(object,...){
n2L<-with(object,.single.b.n2L(y, subject, Bt, tm, tf, Da))
n2L + with(object,2*(.pfda.df(Bt,lm,K)+k*.pfda.df(Bt,lf,K)))
}
single.b<-function(y,t,subject, knots=NULL, penalties=NULL, df=NULL, k=NULL, control=pfdaControl(),subset=NULL){
{ # setup
fname = deparse(match.call()[[1L]])
localfuncs(c('.F.single.optimize.npc','.F.optimize.penalties'))
eval(.X.subset)
eval(.X.single.knots)
eval(.X.single.penalties)
}
if(is.null(k)||any(is.na(k))){
# stop('number of principal components optimization not finished yet')
.F.single.optimize.npc()
} else
if(any(is.na(penalties))) {
funcall <- match.call()
.F.optimize.penalties()
}
else {
eval(.X.binary.y)
rtn<-if(control$useC){
{ # setup for passing to Compiled code
k = as.integer(k[1])
N = nlevels(subject)
nobs= table(subject)
M = length(y)
kr = with(control,max(binary.k0,binary.kr))
p = ncol(Bt)
{ #Compute Memory Requirements
ni = max(nobs)
.inits = 2*p^2 + M + p*N + 8*p
.step.W = ni*kr + kr + (ni * k + ni + k^2 + k + p + p^2 + k*p + 3*k)
.step.1.E = M + M*k + N*k + 2* k^2 + 2 * k + k *ni
.step.2 = p^2+ M+ M*k
.step.3 = M + p^2 + p + k^2
.step.4 = k + 2*k^2 + 2*p^2 + p + p*max(k,8)
dpl <- M + sum(nobs^2) + N*k^2 + N*p^2 + p + p*k + k + N + max(.inits, .step.W, .step.1.E, .step.2, .step.3, .step.4)
ipl <- 8*p
}
}
structure(.C('pfda_bin_single', y=y, nobs=nobs, M=M, N=N, k=k, Bt=Bt, p=p, lm=penalties[1], lf=penalties[2], K=Kt, tm=double(p), tf=matrix(0,p,k), Da=double(k), alpha=matrix(0,N,k), Saa=array(0,dim=c(k,k,N)), minV=control$minimum.variance, tol=control$convergence.tolerance, Iterations=control$max.iterations, burninlength=as.integer(control$binary.burnin), burningenerate=as.integer(control$binary.k0), weightedgenerate=as.integer(control$binary.kr), dl=control$C.debug, dp=double(dpl) , p=integer(ipl))
,class=c('pfda.single.b.rawC','pfda.single.b','list'))
} else {
structure(.single.b.core(y,Bt,subject,k,penalties[1],penalties[2],K=Kt,control$minimum.variance, control$max.iterations,control$convergence.tolerance)
,class=c('pfda.single.b.R','pfda.single.b','list'))
}
rtn$tbase<-tbase
rtn$y<-y
rtn$subject<-subject
return(rtn)
}
}
print.pfda.single.b<-function(x,...){
cat('Univariate Binary Functional Principal Component Model\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(NCOL(x$tf),' principal components\n')
cat('penalties are \n');print( penalty.pfda.single.b(x))
}
penalty.pfda.single.b<-function(x){
with(x,c(mean=lm,pc=lf))
}
}
{ # Dual (Continuous/Continuous) case
.dual.cc.i<-function(y,z,B,subject,ka,kb,min.v){
a<-.single.c.i(y,B,subject,ka,min.v)
b<-.single.c.i(z,B,subject,kb,min.v)
lambda=crossprod(b$alpha,a$alpha)%*%solve(crossprod(a$alpha))
ab<-array(0,dim=c(ka,kb,nlevels(subject)))
for(i in seq_len(nlevels(subject)))ab[,,i]<-tcrossprod(a$alpha[i,],b$alpha[i,])
{ list(tm=a$tm,tn=b$tm,
tf=a$tf,tg=b$tf,
alpha=a$alpha,beta=b$alpha,
Da=a$Da,Db=b$Da,
lambda=lambda,
s.eps=a$sigma,s.xi=b$sigma,
aa=a$aa,ab=ab,bb=b$aa
)}
}
.dual.cc.1.1<-function(yi,zi,Bi,tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi,Setai){
phi = Bi%*%tf
psi = Bi%*%tg
Saa = diag(1/Da,length(Da)) + t(lambda)%*%Setai%*%lambda+crossprod(phi)/s.eps
Sab = crossprod(-lambda,Setai)
Sbb = Setai+crossprod(psi)/s.xi
Sbbi = solve(Sbb)
Zaa = solve(Saa-Sab%*%Sbbi%*%t(Sab))
Zab = -Zaa%*%Sab%*%Sbbi
Zbb = solve(Sbb-t(Sab)%*%solve(Saa)%*%Sab)
mu.a = Zaa%*%crossprod(phi,yi-Bi%*%tm)/s.eps + Zab%*%crossprod(psi,zi-Bi%*%tn)
mu.b = t(Zab)%*%crossprod(phi,yi-Bi%*%tm)/s.eps + Zbb%*%crossprod(psi,zi-Bi%*%tn)
aa = tcrossprod(mu.a)+Zaa
ab = tcrossprod(mu.a,mu.b)+Zab
bb = tcrossprod(mu.b)+Zbb
list(alpha=mu.a, beta=mu.b, Saa=Zaa, Sab=Zab, Sbb = Zbb,aa=aa,ab=ab,bb=bb)
}
.dual.cc.1<-function(y,z,B,subject,tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi){
Setai = solve(Db-lambda%*%Da%*%t(lambda))
#TODO
{
Saa = array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
Sab = array(0,dim=c(NCOL(tf),NCOL(tg),nlevels(subject)))
Sbb = array(0,dim=c(NCOL(tg),NCOL(tg),nlevels(subject)))
alpha = matrix(nrow=nlevels(subject),ncol=NCOL(tf))
beta = matrix(nrow=nlevels(subject),ncol=NCOL(tg))
aa = array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
ab = array(0,dim=c(NCOL(tf),NCOL(tg),nlevels(subject)))
bb = array(0,dim=c(NCOL(tg),NCOL(tg),nlevels(subject)))
}
for(i in seq_len(nlevels(subject))){
ix = i == as.integer(subject)
r<-.dual.cc.1.1(y[ix],z[ix],B[ix,],tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi,Setai)
{ # assign values
Saa[,,i] <-r$Saa
Sab[,,i] <-r$Sab
Sbb[,,i] <-r$Sbb
alpha[i,]<-r$alpha
beta[i,] <-r$beta
aa[,,i] <-r$aa
ab[,,i] <-r$ab
bb[,,i] <-r$bb
}
}
list(alpha=alpha,beta=beta,Saa=Saa,Sab=Sab,Sbb=Sbb,aa=aa,ab=ab,bb=bb)
}
.dual.cc.2<-function(y,z,B,subject,tm,tn,tf,tg,alpha,beta,Saa,Sbb){
list(
s.eps = .single.c.2(y,B,subject,tm,tf,alpha,Saa)$sigma,
s.xi = .single.c.2(z,B,subject,tn,tg,beta ,Sbb)$sigma)
}
.dual.cc.3<-function(y,z,B,subject,tf,tg,alpha,beta,s.eps,s.xi,lm,ln,K){
list(
tm = .gen.tm(y,B,subject,tf,alpha,s.eps,lm,K),
tn = .gen.tm(y,B,subject,tg,beta ,s.xi ,ln,K))
}
.dual.cc.4<-function(y,z,B,subject,tm,tn,tf,tg,alpha,beta,s.eps,s.xi,aa,bb,lf,lg,K){
list(
tf = .gen.tf(y,B,subject,tm,tf,alpha,s.eps,aa,lf,K),
tg = .gen.tf(z,B,subject,tn,tg,beta ,s.xi ,bb,lg,K))
}
.dual.cc.5<-function(aa,ab){
sum.aa<-matrix(0,dim(aa)[1],dim(aa)[2])
sum.ab<-matrix(0,dim(ab)[1],dim(ab)[2])
for(i in seq_len(dim(aa)[3])){
sum.aa = sum.aa + aa[,,i]
sum.ab = sum.ab + ab[,,i]
}
list(lambda = .gen.dual.lambda(sum.aa,sum.ab))
}
.dual.cc.6<-function(tf,tg,alpha,beta,lambda,aa,bb){
sum.aa<-matrix(0,dim(aa)[1],dim(aa)[2])
sum.bb<-matrix(0,dim(bb)[1],dim(bb)[2])
for(sub in seq(NROW(alpha))){
sum.aa = sum.aa + aa[,,sub]
sum.bb = sum.bb + bb[,,sub]
}
a<-.gen.orthog(tf,alpha,sum.aa)
b<-.gen.orthog(tg,beta ,sum.bb)
list(tf = a$tf, tg = b$tf, Da=a$D, Db=b$D, alpha = a$d, beta = b$d,
lambda = b$transformation%*%lambda%*%solve(a$transformation))
}
.dual.cc.n2L<-function(t, y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, seps, sxi){
ka<-length(Da)
kb<-length(Db)
phi<-B%*%tf
psi<-B%*%tg
ll<-0
for(i in seq_len(nlevels(subject))){
ind<-i==as(subject,"integer")
ni<-sum(ind)
Sy<-phi[ind,,drop=FALSE]%*%diag(x=Da,ka,ka)%*%t(phi[ind,,drop=FALSE])+diag(seps,sum(ind))
Sz<-psi[ind,,drop=FALSE]%*%diag(x=Db,kb,kb)%*%t(psi[ind,,drop=FALSE])+diag(sxi,sum(ind))
Syz<-phi[ind,,drop=FALSE]%*%diag(Da,ka,ka)%*%t(lambda)%*%t(psi[ind,,drop=FALSE])
Ry<-y[ind]-B[ind,]%*%tm
Rz<-z[ind]-B[ind,]%*%tn
Cmat<-rbind(cbind(Sy,Syz),cbind(t(Syz),Sz))
Cimat<-solve(Cmat)
Sy.inv <-Cimat[seq_len(ni),seq_len(ni)] #solve(Sy-crossprod(Syz, solve(Sz,t(Syz))))
Sz.inv <- Cimat[seq_len(ni)+ni,seq_len(ni)+ni] #Szi+Szi%*%t(Syz)%*%Sy.inv%*%Syx%*%Szi #solve(Sz-tcrossprod(Syz,solve(Sy,t(Syz))))
Syz.inv<- Cimat[seq_len(ni),seq_len(ni)+ni] #-solve(Sy, Syz)%*%Sz.inv
Dety<-determinant(Sy,logarithm=TRUE)
Detzinv<-determinant(Sz.inv,logarithm=TRUE)
if((Dety$sign<0)|(Detzinv$sign<0))stop("Negative determinant encountered for variance matrix.")
ll<- ll + Dety$modulus-Detzinv$modulus
ll<- ll + crossprod(Ry,Sy.inv%*%Ry)
ll<- ll + crossprod(Rz,Sz.inv%*%Rz)
ll<- ll + 2*crossprod(Ry,Syz.inv%*%Rz)
}
ll
}
AIC.pfda.dual.cc<-function(object,...){
with(object,.dual.ge.AIC(
.dual.cc.n2L(t, y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, seps, sxi),0,
B , B , B , B ,
1 , ka , 1 , kb,
penalties[1] , penalties[2] , penalties[3] , penalties[4],
K , K , K , K ,
seps, seps, sxi, sxi)
)
# .dual.cc.AIC.rawC(object, t, y, z, B, subject, lm, ln, lf, lg, K))
}
logLik.pfda.dual.cc<-function(object,...,newdata=NULL,n2L=TRUE){
with(object,with(newdata,.dual.cc.n2L(t,y,z,B,subject,tm,tn,tf,tg,lambda,Da,Db,seps,sxi)))
}
dual.cc.core<-function(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol){
{ #initial values
init<-.dual.cc.i(y,z,B,subject,ka,kb,min.v)
tm<-init$tm
tn<-init$tn
tf<-init$tf
tg<-init$tg
alpha<-init$alpha
beta<-init$beta
aa<-init$aa
ab<-init$ab
bb<-init$bb
Da<-init$Da
Db<-init$Db
lambda<-init$lambda
n<-nlevels(subject)
Saa<-vector('list',n); for(i in seq_len(n))Saa[[i]]<-matrix(0,ka,ka)
Sab<-vector('list',n); for(i in seq_len(n))Sab[[i]]<-matrix(0,ka,kb)
Sbb<-vector('list',n); for(i in seq_len(n))Sbb[[i]]<-matrix(0,kb,kb)
s.eps<-as.vector(init$s.eps)
s.xi<-as.vector(init$s.xi)
I<-0
nobs<-table(subject)
.cc<-numeric(0)
}
while(I<-I+1){
s1<-.dual.cc.1(y,z,B,subject,tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi)
s2<-.dual.cc.2(y,z,B,subject,tm,tn,tf,tg,s1$alpha,s1$beta,s1$Saa,s1$Sbb)
s3<-.dual.cc.3(y,z,B,subject,tf,tg,alpha,beta,s2$s.eps,s2$s.xi,lm,ln,K)
s4<-.dual.cc.4(y,z,B,subject,s3$tm,s3$tn,tf,tg,s1$alpha,s1$beta,s2$s.eps,s2$s.xi,s1$aa,s1$bb,lf,lg,K)
s5<-.dual.cc.5(s1$aa,s1$ab)
s6<-.dual.cc.6(s4$tf,s4$tg,s1$alpha,s1$beta,s5$lambda,s1$aa,s1$bb)
{ # convergence
ccl<-numeric(0)
ccl['tm']<-sum(abs(tm-s3$tm))
ccl['tn']<-sum(abs(tn-s3$tn))
ccl['tf']<-sum(abs(tf-s6$tf))
ccl['tg']<-sum(abs(tf-s6$tg))
ccl['Da']<-sum(abs(Da-s6$Da))
ccl['Db']<-sum(abs(Da-s6$Db))
ccl['lambda']<-sum(abs(lambda-s6$lambda))
(cc<-sum(ccl))
{ # reassign values
aa<-s1$aa
ab<-s1$ab
bb<-s1$bb
Saa<-s1$Saa
Sab<-s1$Sab
Sbb<-s1$Sab
s.eps<-s2$s.eps
s.xi<-s2$s.xi
tm<-s3$tm
tn<-s3$tn
tf<-s6$tf
tg<-s6$tg
Da<-s6$Da
Db<-s6$Db
alpha<-s6$alpha
beta<-s6$beta
lambda<-s6$lambda
}
if(cc<tol)break
if(I>max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
.cc<-c(.cc,cc)
}
}
list(tm=tm,tn=tn,tf=tf,tg=tg,alpha=alpha,beta=beta,lambda=lambda,Da=Da,Db=Db,s.eps=s.eps,s.xi=s.xi)
}
dual.cc<-function(y,z,t,subject, knots=NULL, penalties=NULL,df=NULL, k=NULL, control=pfdaControl(),tbase=NULL, subset=NULL){
fname = deparse(match.call()[[1L]])
eval(.X.subset)
eval(.X.dual.k)
eval(.X.single.knots)
eval(.X.dual.penalties)
if(any(is.na(k))){
# stop("not finished with number of principal component optimization.")
if(is.na(ka)){
model.y.single <- single.c(y,NULL,t,subject, knots=knots, penalties=penalties[,1],k=NULL, control=control,subset=subset)
ka <- model.y.single$k
}
if(is.na(kb)){
model.z.single <- single.c(z,NULL,t,subject, knots=knots, penalties=penalties[,2],k=NULL, control=control,subset=subset)
kb<-model.z.single$k
}
Recall(y,z,t,subject, knots=knots,penalties=penalties,k=c(ka,kb), control=control)
}
else if(any(is.na(penalties))){
# funcall <- match.call()
# eval(.X.optimize.penalties)}
localfuncs('.F.optimize.penalties')
.F.optimize.penalties()}
else { # Pass to core function
rtn<-{ modifyList(
if(control$useC){
eval(.X.read.parameters)
{ #Compute Memory Requirements
k<-max(ka,kb)
dpl_1 <- M + M*ka + 2*k^2
dpl_2 <- p^2 + M + M*k
dpl_3 <- M + p^2 + p + k^2
dpl_4 <- k^2 + ka*kb
dpl_5_1 <- k + 2*p^2 + p +
max( outer_qf = k*p , eigens = 8*p)
dpl_5_2 <- ka^2 + kb*ka
dpl_5 <- ka^2 + kb^2 +
max(dpl_5_1 ,dpl_5_2 )
dpl_E_1 <- kb^2 + ka^2 +
max( outer_qf = ka*kb,
sym_inv = 2* kb^2,
inner_qf = ka*kb )
dpl_E_2_1 <- 2*k^2
dpl_E_2_2 <- 3*k^2
dpl_E_2 <- kb^2 + max(dpl_E_2_1,dpl_E_2_2,ka*kb)
dpl_E_3_1 <- 2*k*max(nobs)
dpl_E_3 <- dpl_E_3_1
dpl_E <- 2*M + M*ka + M*kb + ka^2 + kb^2 + ka*kb +
max(dpl_E_1, dpl_E_2, dpl_E_3)
dpl<- N*(p**2) + p*2 + p*ka+p*kb + ka*kb + ka + kb +
max(dpl_1, dpl_2, dpl_3, dpl_4, dpl_5, dpl_E )
ipl<-8*p
}
{ structure(.C("pfdaDual",
residuals.y=as.double(y), residuals.z=as.double(z),
B=(Bt),
tm=double(p), tn=double(p),
tf=matrix(0,p,ka), tg=matrix(0,p,kb),
alpha =matrix(0,N,ka),
beta =matrix(0,N,kb),
lambda =matrix(0,kb,ka),
Da =double(ka),
Db =double(kb),
seps = double(1),
sxi =double(1),
Saa =array(0,c(ka,ka,N)),
Sab =array(0,c(ka,kb,N)),
Sbb =array(0,c(kb,kb,N)),
nobs=as.integer(nobs),
M=as.integer(M),
N=N,
ka=as.integer(ka),
kb=as.integer(kb),
p=as.integer(p),
penalties = penalties,
K = Kt,
minV=control$minimum.variance,
Iterations=control$max.iterations,
tol=control$convergence.tolerance,
dl=control$C.debug,
dp=double(dpl),
ip=integer(ipl)
),class=c('pfda.dual.cc.rawC','pfda.dual.cc','list'))}
} else {
structure(dual.cc.core(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol)
,class=c('pfda.dual.cc.R','pfda.dual.cc','list'))
},
list(tbase=tbase,subject=subject,y=y,z=z))}
}
}
print.pfda.dual.cc<-function(x,...){
cat('Paired Functional Principal Component Model (Continuous/Continuous)\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(attr(x,'name.y'),' has ', NCOL(x$tf),' principal components\n')
cat(attr(x,'name.z'),' has ', NCOL(x$tg),' principal components\n')
cat('penalties are \n');print( penalty.pfda.dual.cc(x))
}
penalty.pfda.dual.cc<-function(object,..)with(object,structure(matrix(penalties,2,2),dimnames=list(c(attr(object,'name.y'),attr(object,'name.z')),c('mean','pc'))))
}
{ # Dual(Binary/Continuous) case
.dual.bc.i<-function(y,z,B,subject,ka,kb,min.v){
a<-.single.b.i(y,B,subject,ka,min.v)
b<-.single.c.i(z,B,subject,kb,min.v)
lambda=crossprod(b$alpha,a$alpha)%*%solve(crossprod(a$alpha))
list(tm=a$tm,tn=b$tm,tf=a$tf,tg=b$tf,alpha=a$alpha,beta=b$alpha,aa=a$aa,bb=b$aa,Da=a$Da,Db=b$Da,lambda=lambda,s.xi=b$sigma)
}
.dual.bc.1a<-function(zi,Bi,wi,tm,tn,tf,tg,s.xi,Saa,Sab){
crossprod(Bi%*%tf%*%Saa,wi-Bi%*%tm)+Sab%*%crossprod(Bi%*%tg,zi-Bi%*%tn)/s.xi
}
.dual.bc.1b<-function(zi,Bi,wi,tm,tn,tf,tg,s.xi,Sab,Sbb){
crossprod(Bi%*%tf%*%Sab,wi-Bi%*%tm)+Sbb%*%crossprod(Bi%*%tg,zi-Bi%*%tn)/s.xi
}
.dual.bc.1cde<-function(zi,Bi,wi,wwi,tm,tn,tf,tg,s.xi,Saa,Sab,Sbb){
S1<-.dual.bc.1a(zi,Bi,0,tm,tn,tf,tg,s.xi,Saa,Sab)
S2<-.dual.bc.1b(zi,Bi,0,tm,tn,tf,tg,s.xi,Sab,Sbb)
S3<-tcrossprod(Saa,Bi%*%tf)
S4<-S3%*%wi
S5<-t(Bi%*%tf%*%Sab)
S6<-S5%*%wi
list(
aa= Saa+tcrossprod(S1)+tcrossprod(S1,S4)+tcrossprod(S4,S1) + S3%*%wwi%*%t(S3)
,
ab= Sab+tcrossprod(S1,S2)+tcrossprod(S4,S2)+tcrossprod(S1,S6)+S3%*%tcrossprod(wwi,S5)
,
bb= Sbb+tcrossprod(S2)+tcrossprod(S6,S2)+tcrossprod(S2,S6)+S5%*%tcrossprod(wwi,S5)
)
}
.dual.bc.1<-function(z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,s.xi){
m<-NROW(B)
n<-nlevels(subject)
ka<-NCOL(tf)
kb<-NCOL(tg)
alpha<-matrix(nrow=n,ncol=ka)
beta <-matrix(nrow=n,ncol=kb)
aa <- array(dim=c(ka,ka,n))
ab <- array(dim=c(ka,kb,n))
bb <- array(dim=c(kb,kb,n))
Saa<-vector('list',n)
Sab<-vector('list',n)
Sbb<-vector('list',n)
for(sn in seq_len(n)){
ix<-sn==as.integer(subject)
Bi<-B[ix,]
s<-.gen.dual.sigmas(Bi,Bi,tf,tg,lambda,Da,Db,1,s.xi)
alpha[sn,]<-.dual.bc.1a(z[ix],Bi,w[ix],tm,tn,tf,tg,s.xi,s$Saa,s$Sab)
beta [sn,]<-.dual.bc.1b(z[ix],Bi,w[ix],tm,tn,tf,tg,s.xi,s$Sab,s$Sbb)
cde<-.dual.bc.1cde(z[ix],Bi,w[ix],ww[[sn]],tm,tn,tf,tg,s.xi,s$Saa,s$Sab,s$Sbb)
aa[,,sn]<-cde$aa
ab[,,sn]<-cde$ab
bb[,,sn]<-cde$bb
Saa[[sn]]<-s$Saa
Sab[[sn]]<-s$Sab
Sbb[[sn]]<-s$Sbb
}
list(alpha=alpha,beta=beta,aa=aa,ab=ab,bb=bb,Saa=Saa,Sab=Sab,Sbb=Sbb)
}
.dual.bc.2<-function(z,B,subject,tn,tg,beta,Sbb){
m<-NROW(B)
Rz<-z-B%*%tn
sum.xi<-0
for(i in seq_len(nlevels(subject))){
ix<- i == as.integer(subject)
sum.xi<-sum.xi+sum(diag(B[ix,]%*%tg%*%tcrossprod(Sbb[[i]],B[ix,]%*%tg)))
Rz[ix]<-Rz[ix]-B[ix,]%*%tg%*%beta[i,]
}
list(s.xi=(sum.xi+crossprod(Rz))/m)
}
.dual.bc.3<-function(z,B,subject,w,tf,tg,alpha,beta,s.xi,lm,ln,K){
list(
tm=.gen.tm(w,B,subject,tf,alpha,1 ,lm,K),
tn=.gen.tm(z,B,subject,tg,beta ,s.xi,ln,K)
)
}
.dual.bc.4<-function(z,B,subject,w,tm,tn,tf,tg,alpha,beta,s.xi,aa,bb,lf,lg,K){
list(
tf= .gen.tf(w,B,subject,tm,tf,alpha,1 ,aa,lf,K),
tg= .gen.tf(z,B,subject,tn,tg,beta ,s.xi,bb,lg,K))
}
.dual.bc.5<-function(aa,ab){
.dual.cc.5(aa,ab)
}
.dual.bc.6<-function(tf,tg,alpha,beta,lambda,aa,bb){
.dual.cc.6(tf,tg,alpha,beta,lambda,aa,bb)
}
.dual.bc.genw<-function(j,kr,yi,wi,zi,Bi,tm,tn,tf,tg,Da,Db,lambda){
#
Bij <- Bi[-j,,drop=F]
wij <- wi[-j]
#
r<-.gen.dual.sigmas(Bij,Bi,tf,tg,lambda,Da,Db,1,1)
mu <- r$Saa %*%crossprod(Bij%*%tf,wij-Bij%*%tm) + r$Sab %*% crossprod(Bi%*%tg, zi-Bi%*%tn)
# mu = .dual.bc.1a(zi,Bij,wij,tm,tn,tf,tg,1.0,r$Saa,r$Sab)
# Ss <- solve(crossprod(Bij%*%tf)+solve(diag(Da,length(Da))))
# pfda_bin_single_generate_w_parms2
a <- crossprod(Bi[j,],tm) + t(Bi[j,])%*%tf%*%mu
s <- 1+ t(Bi[j,])%*%tf%*%r$Saa%*%t(tf)%*%Bi[j,]
a+if(yi[j]){
s*.rtruncnormlower(kr,0,1,-a/s)
} else {
-s*.rtruncnormlower(kr,0,1,a/s)
}
}
.dual.bc.w.1<-function(kr,yi,wi,zi,Bi,tm,tn,tf,tg,Da,Db,lambda){
# R version of
n<-NROW(Bi)
w_sim<-sapply(seq_len(n),.dual.bc.genw,kr,yi,wi,zi,Bi,tm,tn,tf,tg,Da,Db,lambda)
wi<-apply(w_sim,2,mean)
wwi<-crossprod(w_sim)/kr
list(wi=wi,wwi=wwi)
}
.dual.bc.w<-function(y,z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,weight,kr){
for(i in seq_len(nlevels(subject))){
ix <- i==as.integer(subject)
rtn <- .dual.bc.w.1(kr,y[ix],w[ix],z[ix],B[ix,,drop=F],tm,tn,tf,tg,Da,Db,lambda)
w[ix] <- (1-weight)*w[ix]+weight*rtn[[1]]
ww[[i]] <- (1-weight)*ww[[i]]+weight*rtn[[2]]
}
return(list(w=w,ww=ww))
}
dual.bc.core<-function(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol){
{ #initial values
init<-.dual.bc.i(y,z,B,subject,ka,kb,min.v)
tm<-init$tm
tn<-init$tn
tf<-init$tf
tg<-init$tg
alpha<-init$alpha
beta<-init$beta
aa<-init$aa
ab<-init$ab
bb<-init$bb
Da<-init$Da
Db<-init$Db
lambda<-init$lambda
n<-nlevels(subject)
Saa<-vector('list',n); for(i in seq_len(n))Saa[[i]]<-matrix(0,ka,ka)
Sab<-vector('list',n); for(i in seq_len(n))Sab[[i]]<-matrix(0,ka,kb)
Sbb<-vector('list',n); for(i in seq_len(n))Sbb[[i]]<-matrix(0,kb,kb)
s.xi<-as.vector(init$s.xi)
I<-0
w<-as.double(y)
nobs<-table(subject)
ww<-vector('list',n); for(i in seq_len(n))ww[[i]]<-matrix(0,nobs[i],nobs[i])
r0=k0=100 ##TODO fix this to read from control values.
kr=10
.cc<-numeric(0)
}
while(I<-I+1){
{ # step w
if(I<r0){
weight=1
rtn<-.dual.bc.w(y,z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,weight,k0)
} else {
weight=10/(10+I)
rtn<-.dual.bc.w(y,z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,weight,kr)
}
w<-rtn$w
ww<-rtn$ww
}
s1<-.dual.bc.1(z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,s.xi)
s2<-.dual.bc.2(z,B,subject,tn,tg,s1$Saa)
s3<-.dual.bc.3(z,B,subject,w,tf,tg,s1$alpha,s2$s.xi,s1$beta,lm,ln,K)
s4<-.dual.bc.4(z,B,w,s3$tm,s3$tn,tf,tg,s1$alpha,s1$beta,s2$s.xi,s1$aa,s1$bb,lf,lg,K)
s5<-.dual.bc.5(s1$aa,s1$ab)
s6<-.dual.bc.6(s4$tf,s4$tg,s1$alpha,s1$beta,s5$lambda,s1$aa,s1$bb)
{ # convergence
ccl<-numeric(0)
ccl['tm']<-sum(abs(tm-s3$tm))
ccl['tn']<-sum(abs(tn-s3$tn))
ccl['tf']<-sum(abs(tf-s6$tf))
ccl['tg']<-sum(abs(tf-s6$tg))
ccl['Da']<-sum(abs(Da-s6$Da))
ccl['Db']<-sum(abs(Da-s6$Db))
ccl['lambda']<-sum(abs(lambda-s6$lambda))
(cc<-sum(ccl))
{ # reassign values
aa<-s1$aa
ab<-s1$ab
bb<-s1$bb
Saa<-s1$Saa
Sab<-s1$Sab
Sbb<-s1$Sab
s.xi<-s2$s.xi
tm<-s3$tm
tn<-s3$tn
tf<-s6$tf
tg<-s6$tg
Da<-s6$Da
Db<-s6$Db
alpha<-s6$alpha
beta<-s6$beta
lambda<-s6$lambda
}
if(cc<tol)break
if(I>max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
.cc<-c(.cc,cc)
}
}
list(tm=tm,tn=tn,tf=tf,tg=tg,alpha=alpha,beta=beta,lambda=lambda,Da=Da,Db=Db,s.xi=s.xi)
}
dual.bc<-function(y,z,t,subject, knots=NULL, penalties=NULL,df=NULL, k=NULL, control=pfdaControl(),tbase=NULL,subset=NULL){
{ #Startup
eval(.X.subset)
eval(.X.dual.k)
eval(.X.single.knots)
eval(.X.dual.penalties)
}
if(any(is.na(k))){
# stop("not finished with number of principal component optimization.")
if(is.na(k[1])) model.y.single <- single.b(y, t,subject, knots=knots, penalties=penalties[,1],k=NULL, control=control,subset=subset)
if(is.na(k[2])) model.z.single <- single.c(z,NULL,t,subject, knots=knots, penalties=penalties[,2],k=NULL, control=control,subset=subset)
Recall(y,z,t,subject, knots=knots,penalties=penalties,k=c(model.y.single$k,model.z.single$k), control=control)
}
else if(any(is.na(penalties))){
funcall <- match.call()
eval(.X.optimize.penalties) }
else { # Pass to core function
eval(.X.binary.y)
rtn<-modifyList(if(control$useC){
eval(.X.read.parameters)
{ #compute memory requirements
k=max(ka,kb)
ni=max(nobs)
kr=max(control$binary.k0, control$binary.kr)
dpl = M + sum(nobs^2) + N*k^2 + N*p^2 + p + p*k + k + max(
M*(3+ka+kb) + (ni*kr+kr+2+5*ka+kb+10*k^2) ,
#ni*kr + kr + (ni * k + ni + k^2 + k + p + p^2 + k*p + 3*k)+1000,# step W
M + M*k + N*k + 2* k^2 + 2 * k + k *ni ,# step 1/E
p^2+ M+ M*k ,# step 2
M + p^2 + p + k^2 ,# step 3
k + 2*k^2 + 2*p^2 + p + p*max(k,8) ,# step 4
2*p^2 + M + p*N + 8*p) # inits
ipl = 8*p
}
{ structure(.C("dual_bc_core",
y =as.integer(y),
z =as.double(z),
B =Bt,
tm =double(p),
tn =double(p),
tf =matrix(0.0,p,ka),
tg =matrix(0.0,p,kb),
alpha =matrix(0.0,N,ka),
beta =matrix(0.0,N,kb),
lambda =matrix(0.0,kb,ka),
Da =double(ka),
Db =double(kb),
sxi =double(1),
aa =array(0.0,dim=c(ka,ka,N)),
ab =array(0.0,dim=c(ka,kb,N)),
bb =array(0.0,dim=c(kb,kb,N)),
Saa =array(0.0,dim=c(ka,ka,N)),
Sab =array(0.0,dim=c(ka,kb,N)),
Sbb =array(0.0,dim=c(kb,kb,N)),
nobs =as.integer(nobs),
N =N,
M =as.integer(M),
ka =as.integer(ka),
kb =as.integer(kb),
p =as.integer(p),
lm =penalties[1],
ln =penalties[2],
lf =penalties[3],
lg =penalties[4],
K = Kt,
minV =control$minimum.variance,
k0 =control$binary.k0,
kr =control$binary.kr,
Iterations=control$max.iterations,
tol =control$convergence.tolerance,
dl =control$C.debug,
dp =double(dpl),
ip =integer(ipl)
),class=c('pfda.dual.bc.rawC','pfda.dual.bc','list'))}
} else {
structure(dual.cc.core(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol)
,class=c('pfda.dual.bc.R','pfda.dual.bc','list'))
},list(subject=subject))
}
}
.dual.bc.updatePCS<-function(Yi, Zi, Da, Db, lambda, phii, psii, Y.rhoi, Z.rhoi, alpha, beta){
ni = NROW(phii)
# mean functions
muy = Y.rhoi+phii%*%alpha
mua = pnorm(muy)
Y.tilda = muy + (Yi-mua)/dnorm(mua)
Z.tilda = Zi
# blocks for means
A = tcrossprod(Da,phii) #' $ A = D_\alpha \phi_i^T $
B = Da %*% t(lambda) %*% t(psii) #' $ B = D_\alpha \Lambda^T \psi_i^T $
C = lambda %*% A #' $ C = \Lambda D_\alpha \psi_i^T = \Lambda A$
D = Db %*% t(psii) #' $ D = D_\beta \psi_i^T $
#blocks associated with variances
E = diag(1, ni) + (phii %*% A) #' $ E = I + \phi_i A = I + \phi_i D_\alpha \phi_i^T $
F = phii %*% B #' $ F = \phi_i D_\alpha \Lambda^T \psi_i $
G = psii %*% C #' $ G = \psi_i \Lambda D_\alpha \phi_i $
H = diag(1, ni) + psii %*% D #' $ H = I + \psi_i D_\beta \psi_i $
# inversion by block
E1 = solve(E)
T = solve(H-G%*%E1%*%F) #' $ T = (H - G E^{-1} F)^{-1} $
Q = E1+E1%*%F%*%T%*%G%*%E1 #' $ Q = E^{-1} + E^{-1} F T G E^{-1} $
R = -E1%*%F%*%T #' $ R = -E^{-1} F T $
S = -T%*%G%*%E1 #' $ S = -T G E^{-1} $
# results
list(
alpha= (A %*% Q + B %*% S) %*% Y.tilda + (A %*% R + B %*% T) %*% Z.tilda,
beta = (C %*% Q + D %*% S) %*% Y.tilda + (C %*% R + D %*% T) %*% Z.tilda)
}
.dual.bc.estimatePCS<-function(Y, Z, B, subject, tm, tn, tf, tg, Da, Db, lambda, sxi){
Y.rho = B%*%tm
Z.rho = B%*%tn
phi = B%*%tf
psi = B%*%tg
ofun<-function(i){
ix = as.integer(subject)==i
alpha = rep(0,NCOL(tf))
beta = rep(0,NCOL(tg))
dif = 1
while(dif>1e-3){
ab<-.dual.bc.updatePCS(Y[ix], Z[ix], Da, Db, lambda, phi[ix,,drop=F], psi[ix,,drop=F], Y.rho[ix,], Z.rho[ix,], alpha, beta)
dif<-mean(alpha-ab$alpha)+mean(beta - ab$beta)
alpha = ab$alpha
beta = ab$beta
}
list(alpha=alpha,beta=beta)
}
ab<-lapply(seq_len(nlevels(subject)), ofun)
alpha<-matrix(0, nlevels(subject), NCOL(tf))
beta <-matrix(0, nlevels(subject), NCOL(tg))
for(i in nlevels(subject)){
alpha[i,]<-ab[[i]]$alpha
beta[i,] <-ab[[i]]$beta
}
list(alpha=alpha,beta=beta)
}
.dual.bc.n2L<-function(y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, sxi){
ka = length(Da)
kb = length(Db)
N = nlevels(subject)
M = length(y)
phi = B %*% tf
psi = B %*% tg
pi = B %*% tm
rho = B %*% tn
ab = .dual.bc.estimatePCS(y, z, B, subject, tm, tn, tf, tg, Da, Db, lambda)
C = rbind(cbind(diag(Da,ka), diag(Da,ka)%*%t(lambda)),cbind(t(lambda%*%diag(Da,ka)), diag(Db,kb)))
N*(sum(log(Da))+determinant(diag(Db,kb)-lambda %*% tcrossprod(diag(Da,ka), lambda))$modulus) +
M*log(sxi) +
with(ab,sum(sapply(seq_len(N),function(i){
ix = subject==i
mu_alpha=pnorm( pi[ix]+phi[ix,,drop=F]%*%alpha[i,])
mu_beta =pnorm(rho[ix]+psi[ix,,drop=F]%*% beta[i,])
W_alpha.i = (dnorm(qnorm(mu_alpha))^2)/mu_alpha*(1-mu_alpha)
gamma<-c(alpha[i,],beta[i,])
# $ \log|I+\Omega W_i \Omega C| $
A<-diag(1,ka)+crossprod(phi[ix,,drop=F], (W_alpha.i * phi[ix,,drop=F]) %*% diag(Da,ka) )
determinant(A)$modulus +
determinant(diag(1,kb) + crossprod(psi[ix,,drop=F]) - lambda %*% (Da * solve(A)) %*% tcrossprod(diag(Da,ka),lambda) )$modulus +
# determinant(diag(1,ka)+ crossprod(phi[ix,],W_alpha.i*phi[ix,,drop=F])%*%C,log=T)$modulus+
# $\sumj^{n_i} Y_i $
sum(-2*(y[ix]*log(mu_alpha) + (1-y[ix])*log(1-mu_alpha)) + (z[ix]-mu_beta)^2) +
{ # $ \gamma_i^T C^{-1} \gamma_i $
Sib <- solve(diag(Db,kb) - lambda %*% Da %*% lambda)
Siab <- -crossprod(lambda ,Sib)
Sia <- diag(1/Da,ka) - Siab %*% lambda %*% diag(1/Da,ka)
crossprod(gamma,solve(C,gamma))
crossprod(alpha[i,],Sia %*% alpha[i,]) +
crossprod(alpha[i,],Siab %*% beta[i,])*2 +
crossprod(beta[i,] ,Sib %*% beta[i,])
}
})))
}
AIC.pfda.dual.bc<-function(object,...){
with(object,.dual.ge.AIC(
.dual.bc.n2L(y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, sxi),0,
B , B , B , B ,
1 , ka , 1 , kb,
lm , ln , lf , lg,
K , K , K , K ,
1, 1, sxi, sxi)
)
}
print.pfda.dual.bc<-function(x,...){
cat('Paired Functional Principal Component Model (Binary/Continuous)\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(attr(x,'name.y'),' has ', NCOL(x$tf),' principal components\n')
cat(attr(x,'name.z'),' has ', NCOL(x$tg),' principal components\n')
cat('penalties are \n');print( penalty.pfda.dual.bc(x))
}
penalty.pfda.dual.bc<-function(object,..)with(object,structure(matrix(c(lm,ln,lf,lg),2,2),dimnames=list(c(attr(object,'name.y'),attr(object,'name.z')),c('mean','pc'))))
}
{ # Additive aka Calcium Model
.dual.ca.i<-function(y,Z,Bt,Bx,subject,kg,kd,min.v){
tz<- if(!is.null(Z) && NCOL(Z)>=1) solve(crossprod(Z),crossprod(Z,y)) else numeric(0)
a<-.single.c.i(y,Bt,subject,kg,min.v)
b<-.single.c.i(y,Bx,subject,kd,min.v)
lambda <- crossprod(b$alpha,a$alpha)%*%solve(crossprod(a$alpha))
gd<-array(0,dim=c(kg,kd,nlevels(subject)))
R=y - Z%*%tz - Bt%*%a$tm - Bx%*%b$tm
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%a$tf%*%a$alpha[i,] - Bx[ix,]%*%b$tf%*%b$alpha[i,]
}
sigma = as.vector(crossprod(R)/length(R))
for(i in seq_len(nlevels(subject)))gd[,,i]<-tcrossprod(a$alpha[i,],b$alpha[i,])
{ list(
tz=tz,
tt=a$tm,tx=b$tm,
tf=a$tf,tg=b$tf,
gamma=a$alpha,delta=b$alpha,
Dg=a$Da,Dd=b$Da,
lambda=lambda,
sigma=sigma,
gg=array(unlist(a$aa),c(kg,kg,nlevels(subject))),
gd=gd,
dd=array(unlist(b$aa),c(kd,kd,nlevels(subject)))
)}
}
.dual.ca.E.1<-function(Ri,tt,tx,sigma, phi, psi, O){
# phi = Bti %*% tf
# psi = Bxi %*% tg
# O is list of block inverted Dg,C,Dd from .gen.symblock.solve
S = .gen.symblock.solve(
O[[1]]+crossprod(phi)/sigma,
O[[2]]+crossprod(phi,psi)/sigma,
O[[3]]+crossprod(psi)/sigma)
a=crossprod(phi,Ri/sigma)
b=crossprod(psi,Ri/sigma)
g=S[[1]]%*%a+S[[2]]%*%b
d=crossprod(S[[2]],a)+S[[3]]%*%b
list(
mu.gamma = as.vector(g),
mu.delta = as.vector(d),
Sigma=S,
gg = tcrossprod(g)+S[[1]],
gd = tcrossprod(g,d)+S[[2]],
dd = tcrossprod(d)+S[[3]]
)
}
.dual.ca.E<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma){
O<-.gen.symblock.solve(
diag(Dg,length(Dg)),
tcrossprod(diag(Dg,length(Dg)),lambda),
diag(Dd,length(Dd)))
{
kg<-NCOL(tf)
kd<-NCOL(tg)
n<-nlevels(subject)
Sgg<-array(0,dim=c(kg,kg,n))
Sgd<-array(0,dim=c(kg,kd,n))
Sdd<-array(0,dim=c(kd,kd,n))
gg <-array(0,dim=c(kg,kg,n))
gd <-array(0,dim=c(kg,kd,n))
dd <-array(0,dim=c(kd,kd,n))
gamma<-matrix(nrow=n,ncol=kg)
delta<-matrix(nrow=n,ncol=kd)
R<-y-Z%*%tz-Bt%*%tt-Bx%*%tx
phi = Bt%*%tf
psi = Bx%*%tg
}
for(i in seq_len(n)){
ix = i==as.integer(subject)
r<-.dual.ca.E.1(R[ix],tt,tx,sigma,phi[ix,,drop=F],psi[ix,,drop=F],O)
{
Sgg[,,i]<-r$Sigma[[1]]
Sgd[,,i]<-r$Sigma[[2]]
Sdd[,,i]<-r$Sigma[[3]]
gg[,,i]<-r$gg
gd[,,i]<-r$gd
dd[,,i]<-r$dd
gamma[i,]<-r$mu.gamma
delta[i,]<-r$mu.delta
}
}
list(Sgg=Sgg,Sgd=Sgd,Sdd=Sdd,gamma=gamma,delta=delta,gg=gg,gd=gd,dd=dd)
}
.dual.ca.unpenalized<-function(y,Z,Bt,Bx,subject,tt,tx,tf,tg,gamma,delta,sigma){
if(is.null(Z)||NCOL(Z)==0)return(list(tz=numeric(0)))
R=y-Bt%*%tt-Bx%*%tx
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%tf%*%gamma[i,] - Bx[ix,]%*%tg%*%delta[i,]
}
list(tz=solve(crossprod(Z),crossprod(Z,R)))
}
.dual.ca.penalized.1<-function(R, B, sigma, l, K){
solve(crossprod(B)+sigma*l*K,crossprod(B,R))
}
.dual.ca.penalized<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,gamma,delta,sigma,lt,lx,Kt,Kx){
R=y-Z%*%tz
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%tf%*%gamma[i,] - Bx[ix,]%*%tg%*%delta[i,]
}
list( tt= .dual.ca.penalized.1(R-Bx%*%tx,Bt,sigma,lt,Kt),
tx= .dual.ca.penalized.1(R-Bt%*%tt,Bx,sigma,lx,Kx))
}
.dual.ca.princcomp<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,gamma,delta,sigma,gg,gd,dd, lf,lg,Kt,Kx){
R = y - Z%*%tz -Bt%*%tt - Bx%*%tx
# tf
for(j in seq_len(NCOL(tf))){
left=sigma*lf*Kt
right=numeric(NROW(tf))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
left = left + gg[j,j,i]*crossprod(Bt[ix,])
right= right+ crossprod(
Bt[ix,,drop=F],
R[ix]*gamma[i,j]-
Bx[ix,,drop=F]%*%tg%*%gd[j,,i]-
(if(NCOL(tf)>1)Bt[ix,,drop=F]%*%tf[,-j]%*%gg[j,-j,i] else 0))
}
tf[,j]= .u.orthogonalize(tf[,seq_len(j-1),drop=FALSE],solve(left,right))
}
# tg
for(j in seq_len(NCOL(tg))){
left=sigma*lg*Kx
right=numeric(NROW(tg))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
left = left + dd[j,j,i]*crossprod(Bx[ix,])
right= right+ crossprod(
Bx[ix,,drop=F],
R[ix]*delta[i,j]-
Bt[ix,,drop=F]%*%tf%*%gd[,j,i] -
(if(NCOL(tg)>1)Bx[ix,,drop=F]%*%tg[,-j]%*%dd[j,-j,i] else 0))
}
tg[,j]= .u.orthogonalize(tg[,seq_len(j-1),drop=FALSE],solve(left,right))
}
# return value
list(tf = tf, tg = tg)
}
.dual.ca.variances<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,gamma,delta,sigma,gg,gd,dd){
R=y - Z%*%tz - Bt%*%tt - Bx%*%tx
sum.gg<-matrix(0,NCOL(tf),NCOL(tf))
sum.gd<-matrix(0,NCOL(tf),NCOL(tg))
sum.dd<-matrix(0,NCOL(tg),NCOL(tg))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%tf%*%gamma[i,] - Bx[ix,]%*%tg%*%delta[i,]
sum.gg<-sum.gg+gg[,,i]
sum.gd<-sum.gd+gd[,,i]
sum.dd<-sum.dd+dd[,,i]
}
sigma = as.vector(crossprod(R)/length(R))
eg<-eigen(sum.gg/nlevels(subject))
ed<-eigen(sum.dd/nlevels(subject))
Tg<-eg$vectors*rep(sign(tf[1,]%*%eg$vectors),each=nrow(eg$vectors))
Td<-ed$vectors*rep(sign(tg[1,]%*%ed$vectors),each=nrow(ed$vectors))
lambda = t(solve(sum.gg,sum.gd))
list(sigma=sigma,
Dg = eg$values,
Dd = ed$values,
lambda = t(Td)%*%lambda%*%Tg,
tf = tf%*%Tg,
tg = tg%*%Td,
gamma = gamma%*%Tg,
delta = delta%*%Td
)
}
.dual.ca.core<-function(y,Z,Bt,Bx,subject,kg,kd,lt,lx,lf,lg,Kt,Kx,min.v,max.I,tol){
{ # initial values
if(is.null(Z)) Z<-matrix(nrow=NROW(y),ncol=0)
init<-.dual.ca.i(y,Z,Bt,Bx,subject,kg,kd,min.v)
init$I<-0
init$cctrace<-numeric(0)
}
within(init,while(I<-I+1){
E <- .dual.ca.E(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma)
U <- .dual.ca.unpenalized(y,Z,Bt,Bx,subject, tt, tx,tf,tg ,E$gamma,E$delta ,sigma)
P <- .dual.ca.penalized(y ,Z,Bt,Bx,subject,U$tz, tt, tx,tf,tg ,E$gamma,E$delta ,sigma,lt,lx,Kt,Kx)
C <- .dual.ca.princcomp(y ,Z,Bt,Bx,subject,U$tz,P$tt,P$tx,tf,tg ,E$gamma,E$delta ,sigma,E$gg,E$gd,E$dd, lf,lg,Kt,Kx)
V <- .dual.ca.variances(y ,Z,Bt,Bx,subject,U$tz,P$tt,P$tx,C$tf,C$tg,E$gamma,E$delta,sigma,E$gg,E$gd,E$dd)
{ #Convergence
ccl<-numeric(0)
ccl['tz']<-sum(abs((tz-U$tz)/tz))
ccl['tt']<-sum(abs((tt-P$tt)/tt))
ccl['tx']<-sum(abs((tx-P$tx)/tx))
ccl['tf']<-sum(abs((tf-V$tf)/tf))
ccl['tg']<-sum(abs((tg-V$tg)/tg))
ccl['Dg']<-sum(abs((Dg-V$Dg)/Dg))
ccl['Dd']<-sum(abs((Dd-V$Dd)/Dd))
(cc<-sum(ccl,na.rm=TRUE))
{ # reassign values
gg <-E$gg
gd <-E$gd
dd <-E$dd
tz <-U$tz
tt <-P$tt
tx <-P$tx
tf <-V$tf
tg <-V$tg
gamma <-V$gamma
delta <-V$delta
sigma <-V$sigma
Dg <-V$Dg
Dd <-V$Dd
lambda<-V$lambda
}
if(cc<tol)break
if(I>max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
cctrace<-c(cctrace,cc)
}
})
}
.dual.ca.n2L<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma,...){
ll<-0
phi = Bt%*%tf
psi = Bx%*%tg
Dg<-diag(Dg,length(Dg))
Dd<-diag(Dd,length(Dd))
Ry<-y-Bt%*%tt-Bx%*%tx-(if(is.null(Z))0 else Z%*%tz)
for(i in nlevels(subject)){
ix = i==as.integer(subject)
U = phi[ix,,drop=F]%*%tcrossprod(Dg,psi[ix,]%*%lambda)
S = phi[ix,,drop=F]%*%tcrossprod(Dg,phi[ix,,drop=F])+ psi[ix,,drop=F]%*%tcrossprod(Dd,psi[ix,,drop=F]) + U + t(U) + diag(sigma,sum(ix))
ll = ll+determinant(S)$modulus+crossprod(Ry[ix],solve(S)%*%Ry[ix])
}
attr(ll,"log")<-TRUE
ll
}
AIC.pfda.additive<-function(object,...){
with(object,.dual.ge.AIC(.dual.ca.n2L(y,Z,Bt,Bx,subject, tz, tt, tx, tf, tg, lambda, Dg, Dd, sigma),
kz,
Bt, Bt, Bx, Bx,
1 , kg, 1 , kd,
lt, lf, lx, lg,
Kt, Kt, Kx, Kx,
sigma, sigma, sigma, sigma))
# with(object,.dual.ca.AIC.rawC(object, y, Z, Bt, Bx, subject, Kt, Kx))
}
logLik.pfda.additive<-function(object,...,newdata=NULL,n2L=TRUE){
with(object,with(newdata,.dual.ca.n2L(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma,...)))
}
dual.ca<-function(y,Z,t,x,subject,knots=NULL,penalties=NULL,df=NULL,k=NULL,bases=NULL,control=pfdaControl()){
# name.t = deparse(substitute(t))
# name.x = deparse(substitute(x))
{ #setup material
Z = if(is.null(Z))matrix(nrow=length(y),ncol=0) else as.matrix(Z)
if(is.null(bases))bases=new.env()
stopifnot(is.environment(bases))
if(!(exists("tbase",envir=bases,inherits=FALSE) && exists("xbase",envir=bases, inherits=FALSE))){ # knot and base identification
if(is.null(knots)){
if(length(control$nknots)==2) { nkt <- control$nknots[1]; nkx <- control$nknots[2] } else nkt <- nkx <- control$nknots[1]
kt<-expand.knots(unique(quantile(t,seq(0,1,length.out=nkt))))
kx<-expand.knots(unique(quantile(x,seq(0,1,length.out=nkt))))
} else if(length(knots)==2) {
if(name.t %in% names(knots)) kt <-knots[[name.t]]
else if('t' %in% names(knots)) kt<-knots[['t']]
else kt <- knots[[1]]
if(name.x %in% names(knots)) kx <-knots[[name.x]]
else if('x' %in% names(knots)) kx<-knots[['x']]
else kx <- knots[[2]]
} else stop("knots needs to be a list of length 2")
assign('tbase',OBasis(kt),envir=bases)
assign('xbase',OBasis(kx),envir=bases)
knots<-list(kt,kx)
# names(knots)<-c(name.t,name.x)
}
{ # evaluate Bases
Bt = evaluate(bases$tbase,t)
Bx = evaluate(bases$xbase,x)
Kt = OuterProdSecondDerivative(bases$tbase)
Kx = OuterProdSecondDerivative(bases$xbase)
Bx<-Bx[,-1] #Removing the first column of Bx is to improve stability and remove interdependence between t and x
Kx<-Kx[-1,-1]
}
eval(.X.dual.k)
eval(.X.dual.penalties)
}
if(any(is.na(k))){
# stop("identification of number of principle components is not done yet.")
k<-c(1,1)
model0<-RecallWith(k=k)
AIC0<-AIC(model0)
while(TRUE){
modelA<-try(RecallWith(k=c(model0$kg+1,model0$kd)),silent=TRUE)
modelB<-try(RecallWith(k=c(model0$kg,model0$kd+1)),silent=TRUE)
AICA<- if(is(modelA,'pfda.additive')) AIC(modelA) else Inf
AICB<- if(is(modelB,'pfda.additive')) AIC(modelB) else Inf
if(AIC0<AICA && AIC0<AICB) break
else if(AICA<AICB){
model0<-modelA
AIC0<-AICA
} else {
model0<-modelB
AIC0<-AICB
}
}
return(model0)
} else
if (any(is.na(penalties))) {
funcall <- structure(match.call(),envir=parent.frame(),expand.dots = FALSE)
if(is.null(control$optim.start)){
lt <- l.from.df(2.1,Bt,Kt)
lx <- l.from.df(2.1,Bx,Kx)
control$optim.start<-c(lt,lx,lt,lx)
}
eval(.X.optimize.penalties) }
else {
rtn<-if(control$useC){
{ # compute memory and parameters
kz = if(is.null(Z))0L else NCOL(Z)
kg = as.integer(k[1])
kd = as.integer(k[2])
N = nlevels(subject)
nobs= table(subject)
M = length(y)
pt = ncol(Bt)
px = ncol(Bx)
p = max(pt, px)
k = max(kg, kd)
dpl = kz + kg + kd + pt +px + pt*kg + px*kd + kg*kd + max(
pt*pt*N + px*px*N + 2*p^2 + M + M*k+ p*N + 8*p,
2^kg^2 + 2*kd^2 + kg*kd+ max(kg,kd)*10 + 2*kg*kd+(1+kg+kd)*M,
M + kz^2 + 10*kz + M*max(kg, kd),
10*p + p^2 +M,
M + p^2 + p + max( 10*pt , M + N*(kg+kd) + M*max(kg,kd) ),
2*kg^2 + 2*kd^2 + M + 10*k + 9*p + 3*p^2 + N*k + M*k
)
}
structure(.C('dual_ca_core', residuals=y, Z=Z, Bt=Bt, Bx=Bx,
tz=double(kz), tt=double(pt), tx=double(px),
tf=matrix(0,pt,kg), tg=matrix(0,px,kd),
gamma=matrix(0,N,kg), delta=matrix(0,N,kd),
lambda=matrix(0,kd,kg), Dg=double(kg), Dd=double(kd),
sigma=0.0,
gg=array(0,c(kg,kg,N)), gd=array(0,c(kg,kd,N)), dd=array(0,c(kd,kd,N)),
Sgg=array(0,c(kg,kg,N)), Sgd=array(0,c(kg,kd,N)), Sdd=array(0,c(kd,kd,N)),
nobs=nobs, N=N, M=M, kz=kz, kg=kg, kd=kd, pt=pt, px=px,
lt=penalties[1], lx=penalties[2], lf=penalties[3], lg=penalties[4], Kt=Kt, Kx=Kx,
minV=control$minimum.variance, Iterations=control$max.iterations, tol=control$convergence.tolerance,
dl=control$C.debug, dp=double(dpl), ip=integer(max(6*p,kz)))
,class=c('pfda.additive.rawC','pfda.additive','list'))
} else {
structure(.dual.ca.core(y,Z,Bt,Bx,subject,k[1],k[2],penalties[1],penalties[2],penalties[3],penalties[4],Kt,Kx,control$minimum.variance,control$max.iterations,control$convergence.tolerance)
,class=c('pfda.additive.R','pfda.additive','list'))
}
rtn$tbase<-bases$tbase
rtn$xbase<-bases$xbase
rtn$y<-y
rtn$subject<-subject
rtn
}
}
plot.pfda.additive<-function(x,...){
with(x,{
layout(matrix(1:4,nrow=2,ncol=2,byrow=T))
plot(tbase,tt, main=paste("Plot of mean curve for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab=attr(x,'name.y'))
plot(tbase,tf, main=paste("Principle components for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab='')
plot(xbase,c(0,tx), main=paste("Plot of mean curve for ",attr(x,'name.x')),xlab=attr(x,'name.x'),ylab=attr(x,'name.y'))
plot(xbase,rbind(0,tg), main=paste("Principle components for ",attr(x,'name.x')),xlab=attr(x,'name.x'),ylab='')
})
}
print.pfda.additive<-function(x,...){
cat('Additive Principal Component Model\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(attr(x,'name.t'),' has ', NCOL(x$tf),' principal components\n')
cat(attr(x,'name.x'),' has ', NCOL(x$tg),' principal components\n')
cat('penalties are \n');print( penalty.pfda.additive(x))
}
penalty.pfda.additive<-function(object,..)with(object,structure(matrix(c(lt,lx,lf,lg),2,2),dimnames=list(c(attr(object,'name.t'),attr(object,'name.x')),c('mean','pc'))))
persp.pfda.additive<-function(x,col.fun=NULL,nt=11,nx=11,...){
if(is.null(col.fun))col.fun<-heat.colors
with(x,{
tx.grid<-expand.grid(
t=seq(tbase@knots[4],tbase@knots[length(tbase@knots)-3],length=nt),
x=seq(xbase@knots[4],xbase@knots[length(xbase@knots)-3],length=nx))
Btg<-evaluate(tbase,tx.grid$t)
Bxg<-evaluate(xbase,tx.grid$x)[,-1]
yg<-matrix(Btg%*%(tt+tf%*%apply(gamma,2,mean)) + Bxg%*%(tx+tg%*%apply(delta,2,mean)),nx,nt,byrow=T)
{
# this code was copied and modified from the persp example page
z <- yg
nrz <- nrow(yg)
ncz <- ncol(yg)
# Create a function interpolating colors in the range of specified colors
nbcol <- 200
color <- col.fun(nbcol)
zfacet <- yg[-1, -1] + yg[-1, -ncz] + yg[-nrz, -1] + yg[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)
}
persp(yg,zlab=attr(x,"name.y"),xlab=attr(x,"name.x"),ylab=attr(x,"name.t"),
col=color[facetcol], ...)
})
}
}
{ # general
.dual.ge.AIC<-function(n2L,kz,
B1, B2, B3, B4,
k1, k2, k3, k4,
l1, l2, l3, l4,
K1, K2, K3, K4,
s1, s2, s3, s4){
as.vector(n2L)+2*(kz +
k1*.pfda.df(B1,l1,K1,s1) +
k2*.pfda.df(B2,l2,K2,s2) +
k3*.pfda.df(B3,l3,K3,s3) +
k4*.pfda.df(B4,l4,K4,s4)
)
}
pfda<-function(model, data=environment(model), ..., driver){
mf <- pfdaParseFormula(model,data)
if(missing(driver))driver = infer.driver(mf)
structure(with(mf,switch(driver,
single.continuous = structure(single.c(response,additive,splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response,'name'),name.t=attr(splinegroup[[1]],'name')),
single.binary = structure(single.b(response,splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response,'name'),name.t=attr(splinegroup[[1]],'name')),
dual.continuous = structure(dual.cc(response[[1]],response[[2]],splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response[[1]],'name'),name.z=attr(response[[2]],'name'),name.t=attr(splinegroup[[1]],'name')),
dual.mixed = structure(dual.bc(response[[1]],response[[2]],splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response[[1]],'name'),name.z=attr(response[[2]],'name'),name.t=attr(splinegroup[[1]],'name'),name.x=attr(splinegroup[[1]],'name')),
additive = structure(dual.ca(response,additive,splinegroup[[1]],splinegroup[[2]],splinegroup[[3]],...),name.y=attr(response,'name'),name.t=attr(splinegroup[[1]],'name'),name.x=attr(splinegroup[[2]],'name'))
)), formula = model)
}
}
Try the pfda package in your browser
Any scripts or data that you put into this service are public.
pfda documentation built on May 2, 2019, 5 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.
{ # global scope variables
# these variables are defined for the sole purpose of confusing the checker
# Each of these should be masked or ignored by a different frame path.
# This is necessary due to the eval/structure programming method employed here.
subject<-
y<-
Z<-
B<-
Bt<-
Bx<-
fname<-
s.xi<-
s.eps<-
ka<-
kb<-
kz<-
nobs<-
M<-
N<-
p<-
K<-
Kt<-
Kx<-
ln<-
lf<-
lg<-
min.v<-
max.I<-
tol<-
name.t<-
name.x<-
tbase<-
xbase<-
NULL
}
{ # general steps
.gen.tm<-function(y,B,subject,tf,alpha,sigma,l,K){
left<-sigma*l*K
right<-matrix(0,NCOL(B),1)
for(i in seq(nlevels(subject))){
ix <- i == as.integer(subject)
Bi<-B[ix,]
left = left + crossprod(Bi)
right = right + crossprod(Bi,y[ix]-Bi%*%tf%*%alpha[i,])
}
solve(left,right)
}
.gen.tf<-function(y,B,subject,tm,tf,alpha,sigma,aa,lf,K){
n<-nlevels(subject)
k<-NCOL(tf)
for(est in seq(k)){
left <- sigma*lf*K
right <- matrix(0,NCOL(B),1)
for(i in seq(n)){
ix <- i == as.integer(subject)
Bi <- B[ix,]
left = left + crossprod(Bi)*(aa[est,est,i])
sne <- matrix(0,NROW(Bi),1)
for(l in setdiff(seq(k),est))sne = sne + Bi%*%tf[,l]*aa[est,l,i]
right = right + crossprod(Bi, (y[ix]-Bi%*%tm)*alpha[i,est] - sne )
}
tf[,est]<-solve(left,right)
}
tf*rep(sign(tf[1,]),each=NROW(tf))
}
.gen.dual.lambda<-function(sum.aa,sum.ab){
t(sum.ab)%*%solve(sum.aa)
}
.gen.orthog<-function(tf,alpha,sum.aa){
n<-NROW(alpha)
k<-NCOL(tf)
Sa<-sum.aa/n
e<-eigen(tf%*%tcrossprod(Sa,tf))
qf = e$vectors[,seq_len(k),drop=FALSE]
qf = qf*rep(sign(qf[1,]),each=NROW(qf))
Da = x=e$values[seq_len(k)]
R = crossprod(qf,tf)
list(tf=qf,D=Da,d=tcrossprod(alpha,R),transformation=R)
}
.gen.dual.sigmas<-function(By,Bz,tf,tg,lambda,Da,Db,s.eps,s.xi){
Seta <- diag(Db,length(Db))-lambda%*%tcrossprod(diag(Da,length(Da)),lambda)
Setai<- solve(Seta)
phi <- By %*% tf
psi <- Bz %*% tg
Sa <- diag(1/Da,length(Da)) + crossprod(lambda,Setai)%*%lambda + crossprod(phi)/s.eps
Sab1 <- -crossprod(lambda,Setai)
Sb <- Setai+crossprod(psi)/s.xi
Saa <- solve(Sa - Sab1%*%solve(Sb,t(Sab1)))
Sab <- -Saa%*%Sab1%*%solve(Sb)
Sbb <- solve(Sb - crossprod(Sab1,solve(Sa,Sab1)))
list(Saa=Saa,Sab=Sab,Sbb=Sbb)
}
.gen.symblock.solve<-function(A,B,C){
#! Solves matrixof form
#! [[ A B]
#! [t(B) C]]
#! Returning the block inverses
#!
stopifnot(NROW(A)==NROW(B))
stopifnot(NCOL(C)==NCOL(B))
svAB = solve(A,B)
Ci = solve(C-crossprod(B,svAB))
Bi = -svAB%*%Ci
Ai = solve(A)+ svAB%*%tcrossprod(Ci,svAB)
list(Ai,Bi,Ci)
}
.X.handle.z<-expression({ # handle Z
stopifnot(exists('Z',inherits=FALSE))
if(missing(Z)||is.null(Z))
Z = matrix(nrow=length(y),ncol=0)
else {
if(is(Z,"formula"))
Z = model.matrix(Z)
else Z = as.matrix(Z)
}
})
.X.subset<-expression({
stopifnot(exists('subset',inherits=FALSE))
if(!missing(subset) && !is.null(subset)){
if(exists('Z',inherits=FALSE))Z<-subset(Z,subset)
y<-subset(y,subset)
if(exists('z',inherits=FALSE))z<-subset(z,subset)
t<-subset(t,subset)
if(exists('x',inherits=FALSE))x<-subset(x,subset)
subject<-subject[subject,drop=T]
}
})
.X.single.knots<-expression({ # knots identification
if(is.list(knots)){
Bt = knots[[1]]
Kt = knots[[2]]
if(!all(dim(Kt)==dim(Bt)[2]))stop("If knots are provided as a list the first component must be the evaluated coefficient matrix for the basis functions, and the second must be the penalty matrix to be used.")
}else {
if(!exists("tbase",inherits=FALSE) || is.null(tbase)){
if(is.null(knots)){
kt<-expand.knots(unique(quantile(t,seq(0,1,length.out=control$nknots))))
tbase = OBasis(kt)
} else if(is(knots,'numeric')) {
kt<-knots
tbase = OBasis(kt)
} else {
tbase = knots
}
Bt = evaluate(tbase,t)
Kt = penaltyMatrix(tbase)
}
}
})
.X.single.penalties<-expression({ # penalties
if(is.null(penalties)){
penalties<- if(is.null(df)) rep(NA,2) else c(l.from.df(df[1],Bt,Kt),l.from.df(df[2],Bt,Kt))
}
})
.X.single.optimize.npc<-expression({
# funcall$k=1
# if(exists('y'))funcall$y=y
# if(exists('Z'))funcall$y=Z
# if(exists('t'))funcall$y=t
# if(exists('subject'))funcall$y=subject
# model.k0<-eval(funcall,env=attr(funcall,'envir'))
message("optimizing number of principal components using AIC")
k=1
model.k0 <- RecallWith(k=k)
while(TRUE){
# funcall$k<-funcall$k+1
# model.k1<-try(eval(funcall,env=attr(funcall,'envir')),silent=TRUE)
k=k+1
model.k1 <- RecallWith(k=k)
if(class(model.k1)[1]=="try-error")return(model.k0)
else if(AIC(model.k0)<AIC(model.k1))return(model.k0)
else model.k0<-model.k1
}
})
.F.single.optimize.npc<-function(){
message("optimizing number of principal components using AIC")
# localfuncs("RecallWith")
k=1
model.k0 <- tryCatch(RecallWith(k=k,fname=fname),warning=function(w)stop(w$message))
tryCatch(while(TRUE){
k=k+1
model.k1 <- RecallWith(k=k,fname=fname)
if(AIC(model.k0)<AIC(model.k1))return(model.k0)
else model.k0<-model.k1
}, warning=function(w){
if(w$message=="EM-algorithm did not converge")model.k0 else stop(w$message)
}, error=function(e)model.k0)
}
.X.dual.k<-expression({ # resloves k input for dual pc (cc/bc/add)
if(is.null(k)) k<-rep(NA,2)
if(length(k)!=2) k<-rep(as.vector(k),length.out=2)
ka = as.integer(k[1])
kb = as.integer(k[2])
# if(names(k)==NULL) names(k)<-c(name.t,name.x)
})
.X.dual.penalties <-expression({ # resolve penalties input
if(is.null(penalties)){
penalties<- if(exists('Bx',inherits=FALSE))
if(is.null(df)) matrix(NA,2,2) else c(l.from.df(df[1],Bt,Kt),l.from.df(df[2],Bx,Kx),l.from.df(df[3],Bt,Kt),l.from.df(df[4],Bx,Kx))
else
if(is.null(df)) matrix(NA,2,2) else c(l.from.df(df[1],Bt,Kt),l.from.df(df[2],Bt,Kt),l.from.df(df[3],Bt,Kt),l.from.df(df[4],Bt,Kt))
}
if(!class(penalties)=='matrix') penalties<-matrix(penalties,2,2)
# if(is.null(colnames(penalties))) colnames(penalties)<-c(name.t,name.x)
})
.X.read.parameters <- expression({
if(exists('Z',inherits=FALSE)){ kz = if(is.null(Z))0L else NCOL(Z)}
k = as.integer(k)
N = nlevels(subject)
nobs= table(subject)
M = length(y)
p = ncol(Bt)
})
.X.binary.y<-expression({ # enforces the conditions on a binary y
y<-if(is(y,'factor')){
stopifnot(nlevels(y)==2)
as.integer(y)-1
} else if(is(y,'logical')){
as.integer(y)
} else {
stopifnot(identical(as.integer(sort(unique(y))), as.integer(0:1)))
as.integer(y)
}
})
.X.funcall.replacevars<-expression({
funcall$y=y
if(exists("z",inherits = FALSE))funcall$z=z
if(exists("Z",inherits = FALSE))funcall$Z=Z
funcall$t = t
if(exists("x",inherits = FALSE))funcall$x = x
funcall$subject = subject
funcall$k=k
funcall$penalties = penalties
funcall$df = df
if(exists("bases",inherits=FALSE))funcall$bases=bases
funcall$knots = knots
funcall$control = control
})
.X.optimize.penalties<-expression({
eval(.X.funcall.replacevars)
pix<-which(is.na(penalties))
if(control$penalty.method=='CV'){
message("optimizing penalties using cross-validation")
if(is.null(control$folds))control$folds<-cv.folds(nlevels(subject),10)
ews<-function(s,p){
ix = !(subject %in% s)
fc<-funcall
fc$penalties=p
fc$subset=ix
model<-eval(fc,env=attr(fc,'envir'))
logLik(model,newdata=list(y=y[!ix],Z=Z[!ix,,drop=FALSE],Bt=Bt[!ix,,drop=FALSE],Bx=Bx[!ix,,drop=FALSE],subject=subject[!ix,drop=TRUE]),n2L=TRUE)
}
cvf<-function(pen,...){
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
cvl<-try(lapply(control$folds,ews,p=p),TRUE)
if(class(cvl)=="try-error")return(NA)
(sum(unlist(cvl)))
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(1,length(pix))
optimpar<-optim(log(control$optim.start),cvf,method=control$optim.method)
penalties[pix]<-exp(optimpar$par)
funcall$penalties=penalties
eval(funcall,env=attr(funcall,'envir'))
} else
if(control$penalty.method=='AIC') {
message("optimizing penalties using AIC")
aicf<-function(pen){
fc<-funcall
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
fc$penalties<-p
m<-try(eval(fc,env=attr(fc,'envir')),silent=TRUE)
if(class(m)[1]=="try-error") NA else AIC(m)
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(l.from.df(2.1,Bt,Kt),length(pix))
optimpar<-optim(log(control$optim.start),aicf,method=control$optim.method)
penalties[pix] <- exp(optimpar$par)
funcall$penalties=penalties
eval(funcall,env=attr(funcall,'envir'))
}
})
.F.optimize.penalties<-function(){
# eval(.X.funcall.replacevars)
pix<-which(is.na(penalties))
if(control$penalty.method=='CV'){
message("optimizing penalties using cross-validation")
if(is.null(control$folds))control$folds<-cv.folds(nlevels(subject),10)
ews<-function(s,p){
ix = !(subject %in% s)
fc<-funcall
fc$penalties=p
fc$subset=ix
model<-eval(fc,env=attr(fc,'envir'))
logLik(model,newdata=list(y=y[!ix],Z=Z[!ix,,drop=FALSE],Bt=Bt[!ix,,drop=FALSE],Bx=Bx[!ix,,drop=FALSE],subject=subject[!ix,drop=TRUE]),n2L=TRUE)
}
cvf<-function(pen,...){
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
cvl<-try(lapply(control$folds,ews,p=p),TRUE)
if(class(cvl)=="try-error")return(NA)
(sum(unlist(cvl)))
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(1,length(pix))
optimpar<-optim(log(control$optim.start),cvf,method=control$optim.method)
penalties[pix]<-exp(optimpar$par)
funcall$penalties=penalties
eval(funcall,env=attr(funcall,'envir'))
} else
if(control$penalty.method=='AIC') {
message("optimizing penalties using AIC")
aicf<-function(pen){
# fc<-funcall
p<-penalties
p[pix]<-exp(pen)
if(any(is.infinite(p)))return(Inf)
# fc$penalties<-p
tryCatch({
m<-RecallWith(penalties=p,fname=fname)
if(class(m)[1]=="try-error") NA else AIC(m)
}
,error=function(e)NA
,warning=function(w){
if(w$message=="EM-algorithm did not converge") return(NA)
else stop(w$message)
}
)
}
if(is.null(control$optim.method))control$optim.method<-"Nelder-Mead"
if(is.null(control$optim.start))control$optim.start<-rep(l.from.df(2.1,Bt,Kt),length(pix))
optimpar<-optim(log(control$optim.start),aicf,method=control$optim.method)
penalties[pix] <- exp(optimpar$par)
# funcall$penalties=penalties
# eval(funcall,env=attr(funcall,'envir'))
RecallWith(penalties=penalties,fname=fname)
}
}
}
{ # utilities and convenience functions
.u.single.resid<-function(y,B,subject,tm,tf,alpha){
for(i in seq_len(nlevels(subject))){
ix = i == as.integer(subject)
y[ix]<-y[ix]-B[ix,,drop=FALSE]%*%(tm+tf%*%alpha[i,])
}
y
}
.u.orthogonalize<-function(U,v){
# U is assumed to be orthogonal already
u = as.vector(if(ncol(U)>0) v-U%*%as.vector(crossprod(v,U)) else v)
u/sqrt(crossprod(u))
}
cv.folds<-function(n, folds = 10){
split(sample(1:n), rep(1:folds, length = n))
}
positive.first.row<-function(X){
stopifnot(is.matrix(X))
X*rep(sign(X[1,]),each=nrow(X))
}
.pfda.df<-Vectorize(function(B,l,K,sigma=1){
tr<-function(x)sum(diag(x))
while(TRUE){
L <- crossprod(B)+sigma*l*K
if(kappa(L)<1e14)break
else if(l>1) l <- l*.90 else NA #l <- l*1.1
}
tr(solve(L,crossprod(B)))
},'l')
l.from.df<-function(df,B,K,sigma=1)if(is.na(df)) NA else uniroot(function(l).pfda.df(B,l,K,sigma)-df,c(1e-4,1e10))$root
RecallWith<-function(...,fname){
calls <- unlist(lapply(sys.calls(),function(x)deparse(x[[1]])[1]))
if(missing(fname)){
iRW<-max(which(calls=="RecallWith"))
if(iRW==1) stop("RecallWith must be called within another function")
unacceptable.calls<-c(which(calls=="{"),unlist(sapply(c("^eval","^function","^try","TryCatch","do.in.envir","^debug:::","^get","assign","mlocal"),grep,calls)))
acceptable.calls<- setdiff(seq_len(iRW-1),unacceptable.calls)
w = max(acceptable.calls)
fname<-calls[w]# deparse(match.call(call = sys.call(w))[[1]])
} else {
w = max(which(calls==fname))
}
oldArgs<-formals(sys.function(w))
if(is.null(oldArgs))oldArgs<-list()
for(a in names(oldArgs)){
try(oldArgs[[a]]<-get(a,sys.frame(w)),silent=TRUE)
}
do.call(fname,as.list(modifyList(oldArgs,list(...))))
}
}
{ # single version
.single.c.i<-function(y,B,subject,k,min.v){
tm = solve(crossprod(B)+diag(min.v,NCOL(B)),crossprod(B,y))
Ry = y-B%*%tm
N = nlevels(subject)
tfa<-matrix(0,NCOL(B),N)
T=matrix(0,NCOL(B),NCOL(B))
for(i in seq_len(nlevels(subject))){
ix <- i==as.integer(subject)
tfa[,i] = solve(crossprod(B[ix,])+min.v*diag(NCOL(B)),crossprod(B[ix,],Ry[ix]))
Ry[ix]<-Ry[ix]-B[ix,]%*%tfa[,i]
T<-T+tcrossprod(tfa[,i])
}
e<-eigen(T)
tf<-e$vectors[,seq_len(k),drop=FALSE]
tf<-tf*rep(sign(tf[1,]),each=NROW(tf))
alpha<-matrix(0,N,k)
aa<-list()
for(i in seq_len(nlevels(subject))){
alpha[i,]<-solve(crossprod(tf),crossprod(tf,tfa[,i]))
aa[[i]] <- tcrossprod(alpha[i,])
}
sigma=crossprod(.u.single.resid(y,B,subject,tm,tf,alpha))/length(y)
return(list(tm=tm,tf=tf,alpha=alpha,aa=aa,Da=e$values[seq_len(k)],sigma=sigma))
}
.single.c.1.1<-function(yi,Bi,tm,tf,Da,sigma){
Sa = solve(crossprod(Bi%*%tf)+diag(1/Da,length(Da)))
mu = Sa%*%crossprod(Bi%*%tf,yi-Bi%*%tm)/sigma
aa = tcrossprod(mu)+Sa
list(alpha=mu,Sa=Sa,aa=aa)
}
.single.c.1<-function(y,B,subject,tm,tf,Da,sigma){
alpha<-matrix(nrow=nlevels(subject),ncol=NCOL(tf))
aa<-array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
Sa<-array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
r<-.single.c.1.1(y[ix],B[ix,],tm,tf,Da,sigma)
alpha[i,]<-r$alpha
Sa[,,i]<-r$Sa
aa[,,i]<-r$aa
}
list(alpha=alpha,Saa=Sa,aa=aa)
}
.single.c.2<-function(y,B,subject,tm,tf,alpha,Sa){
Ry <- y-B%*%tm
s<-0
for(i in seq_len(nlevels(subject))){
ix<- i==as.integer(subject)
phi=B[ix,,drop=FALSE]%*%tf
s<-s+crossprod(Ry[ix]-phi%*%alpha[i,])+sum(diag(phi%*%Sa[,,i]%*%t(phi)))
}
list(sigma=as.vector(s/length(subject)))
}
.single.c.3<-function(y,B,subject,tf,alpha,sigma,lm,K){
list(tm=.gen.tm(y,B,subject,tf,alpha,sigma,lm,K))
}
.single.c.4<-function(y,B,subject,tm,tf,alpha,sigma,aa,lf,K){
list(tf=.gen.tf(y,B,subject,tf,tm,alpha,sigma,aa,lf,K))
}
.single.c.5<-function(tf,alpha,aa){
saa<-matrix(0,NCOL(tf),NCOL(tf))
for(sub in seq(NROW(alpha))) saa = saa + aa[[sub]]
a<-.gen.orthog(tf,alpha,saa)
list(tf=a$tf,Da=a$D, alpha = a$d)
}
single.c.core<-function(y,B,subject,k,lm,lf,K,min.v,max.I,tol){
{ #Initial values
r<-.single.c.i(y,B,subject,k,min.v)
tm=r$tm
tf=r$tf
alpha=r$alpha
aa=r$aa
Da=r$Da
sigma=r$sigma
I=0
.cc<-numeric(0)
}
while(I<-I+1){
s1<-.single.c.1(y,B,subject,tm,tf,Da,sigma)
s2<-.single.c.2(y,B,subject,tm,tf,s1$alpha,s1$Saa)
s3<-.single.c.3(y,B,subject,tf,s1$alpha,s2$sigma,lm,K)
s4<-.single.c.4(y,B,subject,s3$tm,tf,s1$alpha,s2$sigma,s1$aa,lf,K)
s5<-.single.c.5(s4$tf,s1$alpha,s1$aa)
{ # convergence checks
ccl<-numeric(0)
ccl['tm'] <- sum(abs((tm-s3$tm)/s3$tm))
ccl['tf'] <- sum(abs((tf-s5$tf)/s5$tf))
ccl['Da'] <- sum(abs((Da-s5$Da)/s5$Da))
cc <- sum(ccl)
.cc<-c(.cc,cc)
if(cc<tol)break
if(I>=max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
{ # reassign values to
tm<-s3$tm
tf<-s5$tf
alpha<-s5$alpha
Da<-s5$Da
sigma<-s2$sigma
aa<-s5$aa
Saa<-s1$Saa
}
}
}
list(y=y,Bt=B,subject=subject,tm=tm,tf=tf,alpha=alpha, Da=Da,sigma=sigma,aa=aa,Saa=Saa,lm=lm,lf=lf,K=K,k=k,I=I,cc=cc)
}
.single.c.n2L<-function(y, subject, B, tm, tf, Da, sigma){
phi<-B%*%tf
L <- 0.0
for(i in seq_len(nlevels(subject))){
ix<-i==as(subject,"integer")
Sy<-phi[ix,,drop=FALSE]%*%tcrossprod(diag(x=Da,ncol=length(Da)),phi[ix,,drop=FALSE])+diag(sigma,nrow=sum(ix))
Ry<-y[ix]-B[ix,]%*%tm
D<-determinant(Sy,logarithm=TRUE)
if(D$sign<0)stop("Negative determinant encountered for variance matrix.")
L<-L+as.numeric(D$modulus+crossprod(Ry,solve(Sy,Ry)))
}
L
}
AIC.pfda.single.c<-function(object,...){
with(object,{
as.vector(.single.c.n2L(y, subject, Bt, tm, tf, Da, sigma))+
2*(.pfda.df(Bt,lm,K,sigma)+k*.pfda.df(Bt,lf,K,sigma))
})
}
# AIC.pfda.single.c.R<-function(object,...)
logLik.pfda.single.c.R<-logLik.pfda.single.c.rawC<-function(object,...,newdata=NULL, n2L=TRUE){
r<-with(object,with(newdata,.single.c.n2L(y,subject,B,tm,tf,Da,sigma)))
if(n2L) r else exp(-r/2)
}
single.c<-function(y,Z,t,subject,knots=NULL,penalties=NULL,df=NULL,k=NULL,control=pfdaControl(),subset=NULL,...){
{ # setup
fname = deparse(match.call()[[1L]])
localfuncs(c('.F.single.optimize.npc','.F.optimize.penalties'))
name.t = deparse(substitute(t))
eval(.X.handle.z)
eval(.X.subset)
eval(.X.single.knots)
eval(.X.single.penalties)
}
if(is.null(k)||any(is.na(k))){
# funcall <- match.call()
#stop("identification of number of principle components is not done yet.")
# eval(.X.single.optimize.npc)
.F.single.optimize.npc()
} else
if (any(is.na(penalties))) {
funcall <- match.call()
.F.optimize.penalties()
}
else {
rtn<-if(control$useC){
{ # setup for passing to Compiled code
eval(.X.read.parameters)
{ #Compute Memory Requirements
pfda_computeResid = M*k
pfda_s_i = p*p + M + N*p + p*p + max(
pfda_eigens =8*p,
pfda_computeResid )
single_c_resid = pfda_computeResid
pfda_m1 = M + M*k + 2*k^2
pfda_m2 = M + M*k + p^2
pfda_m3_for_subjectnum = p
pfda_m3_for_estcol = p^2 + pfda_m3_for_subjectnum
pfda_m3_core = M + pfda_m3_for_estcol
pfda_m3 = k^2 + pfda_m3_core
pfda_m5_0 = p^2 + p^2 + p + k + max(
pfda_eigens = 8*p,
pfda_matrix_outer_quadratic_form =k*p)
pfda_m5_1 = k^2 + pfda_m5_0
pfdaSingle_m5 = k^2 + pfda_m5_1
pfdaSingle_e_1 = k + 2*k^2
pfdaSingle_e = M + M*k + k*N + max(
pfda_computeResid ,
pfdaSingle_e_1 )
single_c_E = M + pfdaSingle_e
single_c_unpenalized = M + kz*kz + max(
single_c_resid ,
dsysv_=10*kz)
single_c_penalized = M + pfda_m2
single_c_princcomp = M + p*p*N + pfda_m3_core
single_c_variances = M + max(
pfda_m1 ,
pfdaSingle_m5 )
single_c_core = p*p*N + pfda_s_i + p + p*k + k + max(
single_c_E ,
single_c_unpenalized ,
single_c_penalized ,
single_c_princcomp ,
single_c_variances ,
single_c_resid )
ipl<-8*p
dpl<- single_c_core
}
}
structure(.C('single_c_core', residuals=y, Z=Z, Bt=Bt, tz=double(kz), tm=double(p), tf=matrix(0,p,k), alpha=matrix(0,N,k), Da=double(k), sigma=0, aa=array(0,dim=c(k,k,N)), Saa=array(0,dim=c(k,k,N)), nobs=nobs, N=N, M=M, kz=kz, k=k, p=p, lm=penalties[1], lf=penalties[2], K=Kt, minV=control$minimum.variance, Iterations=control$max.iterations, tol=control$convergence.tolerance, dl=control$C.debug, dp=double(dpl), ip=integer(max(6*p,kz)))
,class=c('pfda.single.c.rawC','pfda.single.c','list'))
} else {
structure(single.c.core(y,Bt,subject,k,lm=penalties[1],lf=penalties[2],Kt,control$minimum.variance,control$max.iterations,control$convergence.tolerance)
,class=c('pfda.single.c.R','pfda.single.c','list'))
}
rtn$tbase<-tbase
rtn$y<-y
rtn$subject<-subject
return(rtn)
}
}
print.pfda.single.c<-function(x,...){
cat('Univariate Functional Principal Component Model\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(NCOL(x$tf),' principal components\n')
cat('penalties are \n');print( penalty.pfda.single.c(x))
}
penalty.pfda.single.c<-function(x){
with(x,c(mean=lm,pc=lf))
}
plot.pfda.single.c<-function(x,...){
with(x,{
layout(matrix(1:2,nrow=2,ncol=1,byrow=T))
plot(tbase,tm, main=paste("Plot of mean curve for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab=attr(x,'name.y'))
plot(tbase,tf, main=paste("Principle components for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab='')
})
}
}
{ # single binaries
.roberts1<-function(c){
a<-(c+sqrt(c^2+4))/2
while(TRUE){
u=runif(1)
z=c-log(u)/a
p=exp(-(z-a)^2/2)
v=runif(1)
if(v<=p)return(z)
}
}
.reject1<-function(c){
while(TRUE){
z=rnorm(1)
if(z>c)return(z)
}
}
.rtruncnormlower<-function(m,mean,sd,lower){
c = (lower-mean)/sd
if(c>0){
replicate(m,.roberts1(c))
} else {
replicate(m,.reject1(c))
}
}
.single.b.i<-function(y,B,subject,k,min.v){
.single.c.i(y,B,subject,k,min.v)[1:5]
}
.single.b.1a.Sa.i<-function(B,tf,Da)solve(crossprod(B%*%tf)+solve(Da))
.single.b.1a.mu.i<-function(B,w, tm, tf, Sa)Sa%*%crossprod(B%*%tf, w-B%*%tm)
.single.b.1b.aa.i<-function(B,w,ww,tm,tf,Sa){
phi=B%*%tf
s1=-Sa%*%crossprod(phi,B%*%tm)
s2=Sa%*%crossprod(phi,tcrossprod(as.matrix(w),s1))
Sa%*%t(phi)%*%ww%*%phi%*%t(Sa) +
s2 + t(s2) +
tcrossprod(s1) +
Sa
}
.single.b.1.i<-function(B,w,ww,tm,tf,Da){
Sa = .single.b.1a.Sa.i(B,tf,Da)
alpha = .single.b.1a.mu.i(B, w, tm, tf, Sa)
aa = .single.b.1b.aa.i(B,w,ww,tm,tf,Sa)
list(alpha=alpha,aa=aa,Sa=Sa)
}
.single.b.1<-function(B,subject,w,ww,tm,tf,Da){
subject<-as.integer(as.factor(subject))
n<-max(subject)
k<-NCOL(tf)
aa<-vector('list',n)
alpha<-matrix(nrow=n,ncol=ncol(tf))
Sa<-array(0,dim=c(k,k,n))
for(i in seq(n)){
ix<-i == as.integer(subject)
r<-.single.b.1.i(B[ix,],w[ix],ww[[i]],tm,tf,Da)
Sa[,,i]<-r$Sa
alpha[i,]<-r$alpha
aa[[i]]<-r$aa
}
list(alpha=alpha,aa=structure(unlist(aa),dim=c(1,1,10)),Sa=Sa)
}
.single.b.2<-function(B,subject,w,tf,alpha,lm,K){
list(tm=.gen.tm(w,B,subject,tf,alpha,1,lm,K))
}
.single.b.3<-function(B,subject,w,tm,tf,alpha,aa,lf,K){
list(tf=.gen.tf(w,B,subject,tm,tf,alpha,1,aa,lf,K))
}
.single.b.4<-function(tf,alpha,aa){
saa<-matrix(0,NCOL(tf),NCOL(tf))
for(sub in seq(NROW(alpha))) saa = saa + aa[[sub]]
a<-.gen.orthog(tf,alpha,saa)
list(tf=a$tf,Da=a$D, alpha = a$d)
}
.single.b.w.genw<-function(yi,wi,Bi,tm,tf,Da,kr,j){
# R version of pfda_bin_s_gen_w in C code
Bij <- Bi[-j,,drop=F]
wij <- wi[-j]
# pfda_bin_single_generate_w_parms1
Ss <- solve(crossprod(Bij%*%tf)+solve(diag(Da,length(Da))))
mu <- Ss%*%t(Bij%*%tf)%*%(wij-Bij%*%tm)
# pfda_bin_single_generate_w_parms2
a <- t(Bi[j,])%*%tm + t(Bi[j,])%*%tf%*%mu
s <- 1+t(Bi[j,])%*%tf%*%Ss%*%t(tf)%*%Bi[j,]
a+if(yi[j]) s*.rtruncnormlower(kr,0,1,-a/s) else -s*.rtruncnormlower(kr,0,1,a/s)
}
.single.b.w.1<-function(yi,Bi,wi,tm,tf,Da,kr){
# R version of pfda_bin_single_approximate_moments_forobs
n<-NROW(Bi)
w_sim<-matrix(sapply(seq_len(n),.single.b.w.genw,kr=kr,yi=yi,Bi=Bi,wi=wi,tm=tm,tf=tf,Da=Da),nrow=kr)
wi<-apply(w_sim,2,mean)
# wwa<-array(apply(w_sim,2,tcrossprod),dim=c(n,n,kr))
# wws<-matrix(0,n,n)
# for(l in seq_len(kr))wws<-wws+wwa[,,l]
wwi<-crossprod(w_sim)/kr
return(list(wi=wi,wwi=wwi))
}
.single.b.w<-function(y,B,subject,w,ww,tm,tf,Da,weight,kr){
for(i in seq_len(nlevels(subject))){
ix <- i==as.integer(subject)
rtn <- .single.b.w.1(y[ix],B[ix,,drop=F],w[ix],tm,tf,Da,kr)
w[ix] <- (1-weight)*w[ix]+weight*rtn[[1]]
ww[[i]] <- (1-weight)*ww[[i]]+weight*rtn[[2]]
}
return(list(w=w,ww=ww))
}
.single.b.core<-function(y,B,subject,k,lm,lf,K,minimum.variance, max.iterations,convergence.tolerance){
{ # setup variables & initial values
M<-NROW(y)
N<-nlevels(subject)
nobs<-table(subject)
w <- y
ww<-vector('list',N)
for(i in seq_len(N)) ww[[i]]<-matrix(0,nobs[i],nobs[i])
Saa<-array(0,dim=c(k,k,N))
rtn<-.single.b.i(y,B,subject,k,minimum.variance)
tm<-rtn$tm
tf<-rtn$tf
Da<-rtn$Da
alpha<-rtn$alpha
old<-vector('list',0)
aa <- rtn$aa #vector('list',N); for(i in seq_len(N)) aa<-matrix(0,k,k)
I=0
.cc<-numeric(0)
}
while(I<-I+1){
{ #step w
r0=k0=100
kr=10
if(I<r0){
weight=1
rtn<-.single.b.w(y,B,subject,w,ww,tm,tf,Da,k0,weight)
} else {
weight=10/(10+I)
rtn<-.single.b.w(y,B,subject,w,ww,tm,tf,Da,kr,weight)
}
w<-rtn$w
ww<-rtn$ww
}
s1 <- .single.b.1(B,subject,w,ww,tm,tf,Da)
s2 <- .single.b.2(B,subject,w,tf,s1$alpha,lm,K)
s3 <- .single.b.3(B,subject,w,s2$tm,tf,s1$alpha,s1$aa,lf,K)
s4 <- .single.b.4(tf,alpha,aa)
{ #check convergence
ccl<-numeric(0)
ccl['tm'] <- sum(abs((tm-s2$tm)/s2$tm))
ccl['tf'] <- sum(abs((tf-s4$tf)/s4$tf))
ccl['Da'] <- sum(abs((Da-s4$Da)/s4$Da))
cc <- sum(ccl)
if(cc<convergence.tolerance)break
if(I>=max.iterations){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
.cc<-c(.cc,cc)
{ #reassign variable
tm <- s2$tm
tf <- s4$tf
Da <- s4$Da
alpha <- s4$alpha
aa <- s1$aa
Saa <- s1$Sa
}
}
}
list(tm=tm,tf=tf,Da=Da,aa=aa,Saa=Saa,I=I,cc=cc)
}
.single.b.updatePCS<-function(Yi, Da, phii, Y.rhoi, alpha){
ni = length(Yi)
muy = Y.rhoi+phii%*%alpha
mua = pnorm(muy)
Y.tilda = muy + (Yi-mua)/dnorm(mua)
A = tcrossprod(diag(Da,length(Da)), phii)
Q = solve(diag(1, ni) + phii%*%A)
list(alpha= as.vector((A %*% Q) %*% Y.tilda))
}
.single.b.estimatePCS<-function(Y, B, subject, tm, tf, Da){
Y.rho = B%*%tm
phi = B%*%tf
ofun<-function(i){
ix = as.integer(subject)==i
alpha = rep(0,NCOL(tf))
dif = 1
while(dif>1e-3){
ab<-.single.b.updatePCS(Y[ix], Da, phi[ix,,drop=F], Y.rho[ix,,drop=F], alpha)
dif<-mean(abs(alpha-ab$alpha))
alpha = ab$alpha
}
alpha
}
alpha<-matrix(unlist(lapply(seq_len(nlevels(subject)), ofun)),nlevels(subject),NCOL(tf),byrow=T)
alpha
}
.single.b.n2L.1<-function(yi,Bi,tm,tf,alpha,Da){
sum(.single.b.n2L.parts(yi,Bi,tm,tf,alpha,Da))
}
.single.b.n2L.parts<-function(yi,Bi,tm,tf,alpha,Da){
phi = Bi%*%tf
mui = pnorm(Bi%*%tm+phi%*%alpha)
nui = mui*(1-mui)
Wi = as.vector(dnorm(mui)^2/nui)
c(sum(log(Da)),
determinant(diag(NROW(Da))+crossprod(phi,diag(Wi)%*%phi%*%diag(Da,length(Da))))$modulus,
crossprod(alpha,solve(diag(Da,length(Da)))%*%alpha),
-2*sum(yi*log(mui)+ (1-yi)*log(1-mui)-Wi-mui)
)
}
.single.b.n2L.bySubject<-function(y, subject, B, tm, tf, Da, alpha){
sapply(seq_len(nlevels(subject)),function(si){
ix=as.integer(subject)==si
.single.b.n2L.1(y[ix], B[ix,], tm, tf, alpha[si,], Da)
})
}
.single.b.n2L<-function(y, subject, B, tm, tf, Da){
alpha<-.single.b.estimatePCS(y, B, subject, tm, tf, Da)
ll<-sum(sapply(seq_len(nlevels(subject)),function(si){
ix=as.integer(subject)==si
.single.b.n2L.1(y[ix], B[ix,], tm, tf, alpha[si,], Da)
}))
}
loglik.pfda.single.b<-function(object,...,newdata=NULL,n2L=TRUE){
y<-object@FittingData[[1]][[1]]
t<-object@FittingData[[2]]
subject<-object@FittingData[[3]]
B<-evaluate(object@Basis,t)
}
AIC.pfda.single.b<-function(object,...){
n2L<-with(object,.single.b.n2L(y, subject, Bt, tm, tf, Da))
n2L + with(object,2*(.pfda.df(Bt,lm,K)+k*.pfda.df(Bt,lf,K)))
}
single.b<-function(y,t,subject, knots=NULL, penalties=NULL, df=NULL, k=NULL, control=pfdaControl(),subset=NULL){
{ # setup
fname = deparse(match.call()[[1L]])
localfuncs(c('.F.single.optimize.npc','.F.optimize.penalties'))
eval(.X.subset)
eval(.X.single.knots)
eval(.X.single.penalties)
}
if(is.null(k)||any(is.na(k))){
# stop('number of principal components optimization not finished yet')
.F.single.optimize.npc()
} else
if(any(is.na(penalties))) {
funcall <- match.call()
.F.optimize.penalties()
}
else {
eval(.X.binary.y)
rtn<-if(control$useC){
{ # setup for passing to Compiled code
k = as.integer(k[1])
N = nlevels(subject)
nobs= table(subject)
M = length(y)
kr = with(control,max(binary.k0,binary.kr))
p = ncol(Bt)
{ #Compute Memory Requirements
ni = max(nobs)
.inits = 2*p^2 + M + p*N + 8*p
.step.W = ni*kr + kr + (ni * k + ni + k^2 + k + p + p^2 + k*p + 3*k)
.step.1.E = M + M*k + N*k + 2* k^2 + 2 * k + k *ni
.step.2 = p^2+ M+ M*k
.step.3 = M + p^2 + p + k^2
.step.4 = k + 2*k^2 + 2*p^2 + p + p*max(k,8)
dpl <- M + sum(nobs^2) + N*k^2 + N*p^2 + p + p*k + k + N + max(.inits, .step.W, .step.1.E, .step.2, .step.3, .step.4)
ipl <- 8*p
}
}
structure(.C('pfda_bin_single', y=y, nobs=nobs, M=M, N=N, k=k, Bt=Bt, p=p, lm=penalties[1], lf=penalties[2], K=Kt, tm=double(p), tf=matrix(0,p,k), Da=double(k), alpha=matrix(0,N,k), Saa=array(0,dim=c(k,k,N)), minV=control$minimum.variance, tol=control$convergence.tolerance, Iterations=control$max.iterations, burninlength=as.integer(control$binary.burnin), burningenerate=as.integer(control$binary.k0), weightedgenerate=as.integer(control$binary.kr), dl=control$C.debug, dp=double(dpl) , p=integer(ipl))
,class=c('pfda.single.b.rawC','pfda.single.b','list'))
} else {
structure(.single.b.core(y,Bt,subject,k,penalties[1],penalties[2],K=Kt,control$minimum.variance, control$max.iterations,control$convergence.tolerance)
,class=c('pfda.single.b.R','pfda.single.b','list'))
}
rtn$tbase<-tbase
rtn$y<-y
rtn$subject<-subject
return(rtn)
}
}
print.pfda.single.b<-function(x,...){
cat('Univariate Binary Functional Principal Component Model\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(NCOL(x$tf),' principal components\n')
cat('penalties are \n');print( penalty.pfda.single.b(x))
}
penalty.pfda.single.b<-function(x){
with(x,c(mean=lm,pc=lf))
}
}
{ # Dual (Continuous/Continuous) case
.dual.cc.i<-function(y,z,B,subject,ka,kb,min.v){
a<-.single.c.i(y,B,subject,ka,min.v)
b<-.single.c.i(z,B,subject,kb,min.v)
lambda=crossprod(b$alpha,a$alpha)%*%solve(crossprod(a$alpha))
ab<-array(0,dim=c(ka,kb,nlevels(subject)))
for(i in seq_len(nlevels(subject)))ab[,,i]<-tcrossprod(a$alpha[i,],b$alpha[i,])
{ list(tm=a$tm,tn=b$tm,
tf=a$tf,tg=b$tf,
alpha=a$alpha,beta=b$alpha,
Da=a$Da,Db=b$Da,
lambda=lambda,
s.eps=a$sigma,s.xi=b$sigma,
aa=a$aa,ab=ab,bb=b$aa
)}
}
.dual.cc.1.1<-function(yi,zi,Bi,tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi,Setai){
phi = Bi%*%tf
psi = Bi%*%tg
Saa = diag(1/Da,length(Da)) + t(lambda)%*%Setai%*%lambda+crossprod(phi)/s.eps
Sab = crossprod(-lambda,Setai)
Sbb = Setai+crossprod(psi)/s.xi
Sbbi = solve(Sbb)
Zaa = solve(Saa-Sab%*%Sbbi%*%t(Sab))
Zab = -Zaa%*%Sab%*%Sbbi
Zbb = solve(Sbb-t(Sab)%*%solve(Saa)%*%Sab)
mu.a = Zaa%*%crossprod(phi,yi-Bi%*%tm)/s.eps + Zab%*%crossprod(psi,zi-Bi%*%tn)
mu.b = t(Zab)%*%crossprod(phi,yi-Bi%*%tm)/s.eps + Zbb%*%crossprod(psi,zi-Bi%*%tn)
aa = tcrossprod(mu.a)+Zaa
ab = tcrossprod(mu.a,mu.b)+Zab
bb = tcrossprod(mu.b)+Zbb
list(alpha=mu.a, beta=mu.b, Saa=Zaa, Sab=Zab, Sbb = Zbb,aa=aa,ab=ab,bb=bb)
}
.dual.cc.1<-function(y,z,B,subject,tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi){
Setai = solve(Db-lambda%*%Da%*%t(lambda))
#TODO
{
Saa = array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
Sab = array(0,dim=c(NCOL(tf),NCOL(tg),nlevels(subject)))
Sbb = array(0,dim=c(NCOL(tg),NCOL(tg),nlevels(subject)))
alpha = matrix(nrow=nlevels(subject),ncol=NCOL(tf))
beta = matrix(nrow=nlevels(subject),ncol=NCOL(tg))
aa = array(0,dim=c(NCOL(tf),NCOL(tf),nlevels(subject)))
ab = array(0,dim=c(NCOL(tf),NCOL(tg),nlevels(subject)))
bb = array(0,dim=c(NCOL(tg),NCOL(tg),nlevels(subject)))
}
for(i in seq_len(nlevels(subject))){
ix = i == as.integer(subject)
r<-.dual.cc.1.1(y[ix],z[ix],B[ix,],tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi,Setai)
{ # assign values
Saa[,,i] <-r$Saa
Sab[,,i] <-r$Sab
Sbb[,,i] <-r$Sbb
alpha[i,]<-r$alpha
beta[i,] <-r$beta
aa[,,i] <-r$aa
ab[,,i] <-r$ab
bb[,,i] <-r$bb
}
}
list(alpha=alpha,beta=beta,Saa=Saa,Sab=Sab,Sbb=Sbb,aa=aa,ab=ab,bb=bb)
}
.dual.cc.2<-function(y,z,B,subject,tm,tn,tf,tg,alpha,beta,Saa,Sbb){
list(
s.eps = .single.c.2(y,B,subject,tm,tf,alpha,Saa)$sigma,
s.xi = .single.c.2(z,B,subject,tn,tg,beta ,Sbb)$sigma)
}
.dual.cc.3<-function(y,z,B,subject,tf,tg,alpha,beta,s.eps,s.xi,lm,ln,K){
list(
tm = .gen.tm(y,B,subject,tf,alpha,s.eps,lm,K),
tn = .gen.tm(y,B,subject,tg,beta ,s.xi ,ln,K))
}
.dual.cc.4<-function(y,z,B,subject,tm,tn,tf,tg,alpha,beta,s.eps,s.xi,aa,bb,lf,lg,K){
list(
tf = .gen.tf(y,B,subject,tm,tf,alpha,s.eps,aa,lf,K),
tg = .gen.tf(z,B,subject,tn,tg,beta ,s.xi ,bb,lg,K))
}
.dual.cc.5<-function(aa,ab){
sum.aa<-matrix(0,dim(aa)[1],dim(aa)[2])
sum.ab<-matrix(0,dim(ab)[1],dim(ab)[2])
for(i in seq_len(dim(aa)[3])){
sum.aa = sum.aa + aa[,,i]
sum.ab = sum.ab + ab[,,i]
}
list(lambda = .gen.dual.lambda(sum.aa,sum.ab))
}
.dual.cc.6<-function(tf,tg,alpha,beta,lambda,aa,bb){
sum.aa<-matrix(0,dim(aa)[1],dim(aa)[2])
sum.bb<-matrix(0,dim(bb)[1],dim(bb)[2])
for(sub in seq(NROW(alpha))){
sum.aa = sum.aa + aa[,,sub]
sum.bb = sum.bb + bb[,,sub]
}
a<-.gen.orthog(tf,alpha,sum.aa)
b<-.gen.orthog(tg,beta ,sum.bb)
list(tf = a$tf, tg = b$tf, Da=a$D, Db=b$D, alpha = a$d, beta = b$d,
lambda = b$transformation%*%lambda%*%solve(a$transformation))
}
.dual.cc.n2L<-function(t, y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, seps, sxi){
ka<-length(Da)
kb<-length(Db)
phi<-B%*%tf
psi<-B%*%tg
ll<-0
for(i in seq_len(nlevels(subject))){
ind<-i==as(subject,"integer")
ni<-sum(ind)
Sy<-phi[ind,,drop=FALSE]%*%diag(x=Da,ka,ka)%*%t(phi[ind,,drop=FALSE])+diag(seps,sum(ind))
Sz<-psi[ind,,drop=FALSE]%*%diag(x=Db,kb,kb)%*%t(psi[ind,,drop=FALSE])+diag(sxi,sum(ind))
Syz<-phi[ind,,drop=FALSE]%*%diag(Da,ka,ka)%*%t(lambda)%*%t(psi[ind,,drop=FALSE])
Ry<-y[ind]-B[ind,]%*%tm
Rz<-z[ind]-B[ind,]%*%tn
Cmat<-rbind(cbind(Sy,Syz),cbind(t(Syz),Sz))
Cimat<-solve(Cmat)
Sy.inv <-Cimat[seq_len(ni),seq_len(ni)] #solve(Sy-crossprod(Syz, solve(Sz,t(Syz))))
Sz.inv <- Cimat[seq_len(ni)+ni,seq_len(ni)+ni] #Szi+Szi%*%t(Syz)%*%Sy.inv%*%Syx%*%Szi #solve(Sz-tcrossprod(Syz,solve(Sy,t(Syz))))
Syz.inv<- Cimat[seq_len(ni),seq_len(ni)+ni] #-solve(Sy, Syz)%*%Sz.inv
Dety<-determinant(Sy,logarithm=TRUE)
Detzinv<-determinant(Sz.inv,logarithm=TRUE)
if((Dety$sign<0)|(Detzinv$sign<0))stop("Negative determinant encountered for variance matrix.")
ll<- ll + Dety$modulus-Detzinv$modulus
ll<- ll + crossprod(Ry,Sy.inv%*%Ry)
ll<- ll + crossprod(Rz,Sz.inv%*%Rz)
ll<- ll + 2*crossprod(Ry,Syz.inv%*%Rz)
}
ll
}
AIC.pfda.dual.cc<-function(object,...){
with(object,.dual.ge.AIC(
.dual.cc.n2L(t, y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, seps, sxi),0,
B , B , B , B ,
1 , ka , 1 , kb,
penalties[1] , penalties[2] , penalties[3] , penalties[4],
K , K , K , K ,
seps, seps, sxi, sxi)
)
# .dual.cc.AIC.rawC(object, t, y, z, B, subject, lm, ln, lf, lg, K))
}
logLik.pfda.dual.cc<-function(object,...,newdata=NULL,n2L=TRUE){
with(object,with(newdata,.dual.cc.n2L(t,y,z,B,subject,tm,tn,tf,tg,lambda,Da,Db,seps,sxi)))
}
dual.cc.core<-function(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol){
{ #initial values
init<-.dual.cc.i(y,z,B,subject,ka,kb,min.v)
tm<-init$tm
tn<-init$tn
tf<-init$tf
tg<-init$tg
alpha<-init$alpha
beta<-init$beta
aa<-init$aa
ab<-init$ab
bb<-init$bb
Da<-init$Da
Db<-init$Db
lambda<-init$lambda
n<-nlevels(subject)
Saa<-vector('list',n); for(i in seq_len(n))Saa[[i]]<-matrix(0,ka,ka)
Sab<-vector('list',n); for(i in seq_len(n))Sab[[i]]<-matrix(0,ka,kb)
Sbb<-vector('list',n); for(i in seq_len(n))Sbb[[i]]<-matrix(0,kb,kb)
s.eps<-as.vector(init$s.eps)
s.xi<-as.vector(init$s.xi)
I<-0
nobs<-table(subject)
.cc<-numeric(0)
}
while(I<-I+1){
s1<-.dual.cc.1(y,z,B,subject,tm,tn,tf,tg,lambda,Da,Db,s.eps,s.xi)
s2<-.dual.cc.2(y,z,B,subject,tm,tn,tf,tg,s1$alpha,s1$beta,s1$Saa,s1$Sbb)
s3<-.dual.cc.3(y,z,B,subject,tf,tg,alpha,beta,s2$s.eps,s2$s.xi,lm,ln,K)
s4<-.dual.cc.4(y,z,B,subject,s3$tm,s3$tn,tf,tg,s1$alpha,s1$beta,s2$s.eps,s2$s.xi,s1$aa,s1$bb,lf,lg,K)
s5<-.dual.cc.5(s1$aa,s1$ab)
s6<-.dual.cc.6(s4$tf,s4$tg,s1$alpha,s1$beta,s5$lambda,s1$aa,s1$bb)
{ # convergence
ccl<-numeric(0)
ccl['tm']<-sum(abs(tm-s3$tm))
ccl['tn']<-sum(abs(tn-s3$tn))
ccl['tf']<-sum(abs(tf-s6$tf))
ccl['tg']<-sum(abs(tf-s6$tg))
ccl['Da']<-sum(abs(Da-s6$Da))
ccl['Db']<-sum(abs(Da-s6$Db))
ccl['lambda']<-sum(abs(lambda-s6$lambda))
(cc<-sum(ccl))
{ # reassign values
aa<-s1$aa
ab<-s1$ab
bb<-s1$bb
Saa<-s1$Saa
Sab<-s1$Sab
Sbb<-s1$Sab
s.eps<-s2$s.eps
s.xi<-s2$s.xi
tm<-s3$tm
tn<-s3$tn
tf<-s6$tf
tg<-s6$tg
Da<-s6$Da
Db<-s6$Db
alpha<-s6$alpha
beta<-s6$beta
lambda<-s6$lambda
}
if(cc<tol)break
if(I>max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
.cc<-c(.cc,cc)
}
}
list(tm=tm,tn=tn,tf=tf,tg=tg,alpha=alpha,beta=beta,lambda=lambda,Da=Da,Db=Db,s.eps=s.eps,s.xi=s.xi)
}
dual.cc<-function(y,z,t,subject, knots=NULL, penalties=NULL,df=NULL, k=NULL, control=pfdaControl(),tbase=NULL, subset=NULL){
fname = deparse(match.call()[[1L]])
eval(.X.subset)
eval(.X.dual.k)
eval(.X.single.knots)
eval(.X.dual.penalties)
if(any(is.na(k))){
# stop("not finished with number of principal component optimization.")
if(is.na(ka)){
model.y.single <- single.c(y,NULL,t,subject, knots=knots, penalties=penalties[,1],k=NULL, control=control,subset=subset)
ka <- model.y.single$k
}
if(is.na(kb)){
model.z.single <- single.c(z,NULL,t,subject, knots=knots, penalties=penalties[,2],k=NULL, control=control,subset=subset)
kb<-model.z.single$k
}
Recall(y,z,t,subject, knots=knots,penalties=penalties,k=c(ka,kb), control=control)
}
else if(any(is.na(penalties))){
# funcall <- match.call()
# eval(.X.optimize.penalties)}
localfuncs('.F.optimize.penalties')
.F.optimize.penalties()}
else { # Pass to core function
rtn<-{ modifyList(
if(control$useC){
eval(.X.read.parameters)
{ #Compute Memory Requirements
k<-max(ka,kb)
dpl_1 <- M + M*ka + 2*k^2
dpl_2 <- p^2 + M + M*k
dpl_3 <- M + p^2 + p + k^2
dpl_4 <- k^2 + ka*kb
dpl_5_1 <- k + 2*p^2 + p +
max( outer_qf = k*p , eigens = 8*p)
dpl_5_2 <- ka^2 + kb*ka
dpl_5 <- ka^2 + kb^2 +
max(dpl_5_1 ,dpl_5_2 )
dpl_E_1 <- kb^2 + ka^2 +
max( outer_qf = ka*kb,
sym_inv = 2* kb^2,
inner_qf = ka*kb )
dpl_E_2_1 <- 2*k^2
dpl_E_2_2 <- 3*k^2
dpl_E_2 <- kb^2 + max(dpl_E_2_1,dpl_E_2_2,ka*kb)
dpl_E_3_1 <- 2*k*max(nobs)
dpl_E_3 <- dpl_E_3_1
dpl_E <- 2*M + M*ka + M*kb + ka^2 + kb^2 + ka*kb +
max(dpl_E_1, dpl_E_2, dpl_E_3)
dpl<- N*(p**2) + p*2 + p*ka+p*kb + ka*kb + ka + kb +
max(dpl_1, dpl_2, dpl_3, dpl_4, dpl_5, dpl_E )
ipl<-8*p
}
{ structure(.C("pfdaDual",
residuals.y=as.double(y), residuals.z=as.double(z),
B=(Bt),
tm=double(p), tn=double(p),
tf=matrix(0,p,ka), tg=matrix(0,p,kb),
alpha =matrix(0,N,ka),
beta =matrix(0,N,kb),
lambda =matrix(0,kb,ka),
Da =double(ka),
Db =double(kb),
seps = double(1),
sxi =double(1),
Saa =array(0,c(ka,ka,N)),
Sab =array(0,c(ka,kb,N)),
Sbb =array(0,c(kb,kb,N)),
nobs=as.integer(nobs),
M=as.integer(M),
N=N,
ka=as.integer(ka),
kb=as.integer(kb),
p=as.integer(p),
penalties = penalties,
K = Kt,
minV=control$minimum.variance,
Iterations=control$max.iterations,
tol=control$convergence.tolerance,
dl=control$C.debug,
dp=double(dpl),
ip=integer(ipl)
),class=c('pfda.dual.cc.rawC','pfda.dual.cc','list'))}
} else {
structure(dual.cc.core(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol)
,class=c('pfda.dual.cc.R','pfda.dual.cc','list'))
},
list(tbase=tbase,subject=subject,y=y,z=z))}
}
}
print.pfda.dual.cc<-function(x,...){
cat('Paired Functional Principal Component Model (Continuous/Continuous)\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(attr(x,'name.y'),' has ', NCOL(x$tf),' principal components\n')
cat(attr(x,'name.z'),' has ', NCOL(x$tg),' principal components\n')
cat('penalties are \n');print( penalty.pfda.dual.cc(x))
}
penalty.pfda.dual.cc<-function(object,..)with(object,structure(matrix(penalties,2,2),dimnames=list(c(attr(object,'name.y'),attr(object,'name.z')),c('mean','pc'))))
}
{ # Dual(Binary/Continuous) case
.dual.bc.i<-function(y,z,B,subject,ka,kb,min.v){
a<-.single.b.i(y,B,subject,ka,min.v)
b<-.single.c.i(z,B,subject,kb,min.v)
lambda=crossprod(b$alpha,a$alpha)%*%solve(crossprod(a$alpha))
list(tm=a$tm,tn=b$tm,tf=a$tf,tg=b$tf,alpha=a$alpha,beta=b$alpha,aa=a$aa,bb=b$aa,Da=a$Da,Db=b$Da,lambda=lambda,s.xi=b$sigma)
}
.dual.bc.1a<-function(zi,Bi,wi,tm,tn,tf,tg,s.xi,Saa,Sab){
crossprod(Bi%*%tf%*%Saa,wi-Bi%*%tm)+Sab%*%crossprod(Bi%*%tg,zi-Bi%*%tn)/s.xi
}
.dual.bc.1b<-function(zi,Bi,wi,tm,tn,tf,tg,s.xi,Sab,Sbb){
crossprod(Bi%*%tf%*%Sab,wi-Bi%*%tm)+Sbb%*%crossprod(Bi%*%tg,zi-Bi%*%tn)/s.xi
}
.dual.bc.1cde<-function(zi,Bi,wi,wwi,tm,tn,tf,tg,s.xi,Saa,Sab,Sbb){
S1<-.dual.bc.1a(zi,Bi,0,tm,tn,tf,tg,s.xi,Saa,Sab)
S2<-.dual.bc.1b(zi,Bi,0,tm,tn,tf,tg,s.xi,Sab,Sbb)
S3<-tcrossprod(Saa,Bi%*%tf)
S4<-S3%*%wi
S5<-t(Bi%*%tf%*%Sab)
S6<-S5%*%wi
list(
aa= Saa+tcrossprod(S1)+tcrossprod(S1,S4)+tcrossprod(S4,S1) + S3%*%wwi%*%t(S3)
,
ab= Sab+tcrossprod(S1,S2)+tcrossprod(S4,S2)+tcrossprod(S1,S6)+S3%*%tcrossprod(wwi,S5)
,
bb= Sbb+tcrossprod(S2)+tcrossprod(S6,S2)+tcrossprod(S2,S6)+S5%*%tcrossprod(wwi,S5)
)
}
.dual.bc.1<-function(z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,s.xi){
m<-NROW(B)
n<-nlevels(subject)
ka<-NCOL(tf)
kb<-NCOL(tg)
alpha<-matrix(nrow=n,ncol=ka)
beta <-matrix(nrow=n,ncol=kb)
aa <- array(dim=c(ka,ka,n))
ab <- array(dim=c(ka,kb,n))
bb <- array(dim=c(kb,kb,n))
Saa<-vector('list',n)
Sab<-vector('list',n)
Sbb<-vector('list',n)
for(sn in seq_len(n)){
ix<-sn==as.integer(subject)
Bi<-B[ix,]
s<-.gen.dual.sigmas(Bi,Bi,tf,tg,lambda,Da,Db,1,s.xi)
alpha[sn,]<-.dual.bc.1a(z[ix],Bi,w[ix],tm,tn,tf,tg,s.xi,s$Saa,s$Sab)
beta [sn,]<-.dual.bc.1b(z[ix],Bi,w[ix],tm,tn,tf,tg,s.xi,s$Sab,s$Sbb)
cde<-.dual.bc.1cde(z[ix],Bi,w[ix],ww[[sn]],tm,tn,tf,tg,s.xi,s$Saa,s$Sab,s$Sbb)
aa[,,sn]<-cde$aa
ab[,,sn]<-cde$ab
bb[,,sn]<-cde$bb
Saa[[sn]]<-s$Saa
Sab[[sn]]<-s$Sab
Sbb[[sn]]<-s$Sbb
}
list(alpha=alpha,beta=beta,aa=aa,ab=ab,bb=bb,Saa=Saa,Sab=Sab,Sbb=Sbb)
}
.dual.bc.2<-function(z,B,subject,tn,tg,beta,Sbb){
m<-NROW(B)
Rz<-z-B%*%tn
sum.xi<-0
for(i in seq_len(nlevels(subject))){
ix<- i == as.integer(subject)
sum.xi<-sum.xi+sum(diag(B[ix,]%*%tg%*%tcrossprod(Sbb[[i]],B[ix,]%*%tg)))
Rz[ix]<-Rz[ix]-B[ix,]%*%tg%*%beta[i,]
}
list(s.xi=(sum.xi+crossprod(Rz))/m)
}
.dual.bc.3<-function(z,B,subject,w,tf,tg,alpha,beta,s.xi,lm,ln,K){
list(
tm=.gen.tm(w,B,subject,tf,alpha,1 ,lm,K),
tn=.gen.tm(z,B,subject,tg,beta ,s.xi,ln,K)
)
}
.dual.bc.4<-function(z,B,subject,w,tm,tn,tf,tg,alpha,beta,s.xi,aa,bb,lf,lg,K){
list(
tf= .gen.tf(w,B,subject,tm,tf,alpha,1 ,aa,lf,K),
tg= .gen.tf(z,B,subject,tn,tg,beta ,s.xi,bb,lg,K))
}
.dual.bc.5<-function(aa,ab){
.dual.cc.5(aa,ab)
}
.dual.bc.6<-function(tf,tg,alpha,beta,lambda,aa,bb){
.dual.cc.6(tf,tg,alpha,beta,lambda,aa,bb)
}
.dual.bc.genw<-function(j,kr,yi,wi,zi,Bi,tm,tn,tf,tg,Da,Db,lambda){
#
Bij <- Bi[-j,,drop=F]
wij <- wi[-j]
#
r<-.gen.dual.sigmas(Bij,Bi,tf,tg,lambda,Da,Db,1,1)
mu <- r$Saa %*%crossprod(Bij%*%tf,wij-Bij%*%tm) + r$Sab %*% crossprod(Bi%*%tg, zi-Bi%*%tn)
# mu = .dual.bc.1a(zi,Bij,wij,tm,tn,tf,tg,1.0,r$Saa,r$Sab)
# Ss <- solve(crossprod(Bij%*%tf)+solve(diag(Da,length(Da))))
# pfda_bin_single_generate_w_parms2
a <- crossprod(Bi[j,],tm) + t(Bi[j,])%*%tf%*%mu
s <- 1+ t(Bi[j,])%*%tf%*%r$Saa%*%t(tf)%*%Bi[j,]
a+if(yi[j]){
s*.rtruncnormlower(kr,0,1,-a/s)
} else {
-s*.rtruncnormlower(kr,0,1,a/s)
}
}
.dual.bc.w.1<-function(kr,yi,wi,zi,Bi,tm,tn,tf,tg,Da,Db,lambda){
# R version of
n<-NROW(Bi)
w_sim<-sapply(seq_len(n),.dual.bc.genw,kr,yi,wi,zi,Bi,tm,tn,tf,tg,Da,Db,lambda)
wi<-apply(w_sim,2,mean)
wwi<-crossprod(w_sim)/kr
list(wi=wi,wwi=wwi)
}
.dual.bc.w<-function(y,z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,weight,kr){
for(i in seq_len(nlevels(subject))){
ix <- i==as.integer(subject)
rtn <- .dual.bc.w.1(kr,y[ix],w[ix],z[ix],B[ix,,drop=F],tm,tn,tf,tg,Da,Db,lambda)
w[ix] <- (1-weight)*w[ix]+weight*rtn[[1]]
ww[[i]] <- (1-weight)*ww[[i]]+weight*rtn[[2]]
}
return(list(w=w,ww=ww))
}
dual.bc.core<-function(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol){
{ #initial values
init<-.dual.bc.i(y,z,B,subject,ka,kb,min.v)
tm<-init$tm
tn<-init$tn
tf<-init$tf
tg<-init$tg
alpha<-init$alpha
beta<-init$beta
aa<-init$aa
ab<-init$ab
bb<-init$bb
Da<-init$Da
Db<-init$Db
lambda<-init$lambda
n<-nlevels(subject)
Saa<-vector('list',n); for(i in seq_len(n))Saa[[i]]<-matrix(0,ka,ka)
Sab<-vector('list',n); for(i in seq_len(n))Sab[[i]]<-matrix(0,ka,kb)
Sbb<-vector('list',n); for(i in seq_len(n))Sbb[[i]]<-matrix(0,kb,kb)
s.xi<-as.vector(init$s.xi)
I<-0
w<-as.double(y)
nobs<-table(subject)
ww<-vector('list',n); for(i in seq_len(n))ww[[i]]<-matrix(0,nobs[i],nobs[i])
r0=k0=100 ##TODO fix this to read from control values.
kr=10
.cc<-numeric(0)
}
while(I<-I+1){
{ # step w
if(I<r0){
weight=1
rtn<-.dual.bc.w(y,z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,weight,k0)
} else {
weight=10/(10+I)
rtn<-.dual.bc.w(y,z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,weight,kr)
}
w<-rtn$w
ww<-rtn$ww
}
s1<-.dual.bc.1(z,B,subject,w,ww,tm,tn,tf,tg,lambda,Da,Db,s.xi)
s2<-.dual.bc.2(z,B,subject,tn,tg,s1$Saa)
s3<-.dual.bc.3(z,B,subject,w,tf,tg,s1$alpha,s2$s.xi,s1$beta,lm,ln,K)
s4<-.dual.bc.4(z,B,w,s3$tm,s3$tn,tf,tg,s1$alpha,s1$beta,s2$s.xi,s1$aa,s1$bb,lf,lg,K)
s5<-.dual.bc.5(s1$aa,s1$ab)
s6<-.dual.bc.6(s4$tf,s4$tg,s1$alpha,s1$beta,s5$lambda,s1$aa,s1$bb)
{ # convergence
ccl<-numeric(0)
ccl['tm']<-sum(abs(tm-s3$tm))
ccl['tn']<-sum(abs(tn-s3$tn))
ccl['tf']<-sum(abs(tf-s6$tf))
ccl['tg']<-sum(abs(tf-s6$tg))
ccl['Da']<-sum(abs(Da-s6$Da))
ccl['Db']<-sum(abs(Da-s6$Db))
ccl['lambda']<-sum(abs(lambda-s6$lambda))
(cc<-sum(ccl))
{ # reassign values
aa<-s1$aa
ab<-s1$ab
bb<-s1$bb
Saa<-s1$Saa
Sab<-s1$Sab
Sbb<-s1$Sab
s.xi<-s2$s.xi
tm<-s3$tm
tn<-s3$tn
tf<-s6$tf
tg<-s6$tg
Da<-s6$Da
Db<-s6$Db
alpha<-s6$alpha
beta<-s6$beta
lambda<-s6$lambda
}
if(cc<tol)break
if(I>max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
.cc<-c(.cc,cc)
}
}
list(tm=tm,tn=tn,tf=tf,tg=tg,alpha=alpha,beta=beta,lambda=lambda,Da=Da,Db=Db,s.xi=s.xi)
}
dual.bc<-function(y,z,t,subject, knots=NULL, penalties=NULL,df=NULL, k=NULL, control=pfdaControl(),tbase=NULL,subset=NULL){
{ #Startup
eval(.X.subset)
eval(.X.dual.k)
eval(.X.single.knots)
eval(.X.dual.penalties)
}
if(any(is.na(k))){
# stop("not finished with number of principal component optimization.")
if(is.na(k[1])) model.y.single <- single.b(y, t,subject, knots=knots, penalties=penalties[,1],k=NULL, control=control,subset=subset)
if(is.na(k[2])) model.z.single <- single.c(z,NULL,t,subject, knots=knots, penalties=penalties[,2],k=NULL, control=control,subset=subset)
Recall(y,z,t,subject, knots=knots,penalties=penalties,k=c(model.y.single$k,model.z.single$k), control=control)
}
else if(any(is.na(penalties))){
funcall <- match.call()
eval(.X.optimize.penalties) }
else { # Pass to core function
eval(.X.binary.y)
rtn<-modifyList(if(control$useC){
eval(.X.read.parameters)
{ #compute memory requirements
k=max(ka,kb)
ni=max(nobs)
kr=max(control$binary.k0, control$binary.kr)
dpl = M + sum(nobs^2) + N*k^2 + N*p^2 + p + p*k + k + max(
M*(3+ka+kb) + (ni*kr+kr+2+5*ka+kb+10*k^2) ,
#ni*kr + kr + (ni * k + ni + k^2 + k + p + p^2 + k*p + 3*k)+1000,# step W
M + M*k + N*k + 2* k^2 + 2 * k + k *ni ,# step 1/E
p^2+ M+ M*k ,# step 2
M + p^2 + p + k^2 ,# step 3
k + 2*k^2 + 2*p^2 + p + p*max(k,8) ,# step 4
2*p^2 + M + p*N + 8*p) # inits
ipl = 8*p
}
{ structure(.C("dual_bc_core",
y =as.integer(y),
z =as.double(z),
B =Bt,
tm =double(p),
tn =double(p),
tf =matrix(0.0,p,ka),
tg =matrix(0.0,p,kb),
alpha =matrix(0.0,N,ka),
beta =matrix(0.0,N,kb),
lambda =matrix(0.0,kb,ka),
Da =double(ka),
Db =double(kb),
sxi =double(1),
aa =array(0.0,dim=c(ka,ka,N)),
ab =array(0.0,dim=c(ka,kb,N)),
bb =array(0.0,dim=c(kb,kb,N)),
Saa =array(0.0,dim=c(ka,ka,N)),
Sab =array(0.0,dim=c(ka,kb,N)),
Sbb =array(0.0,dim=c(kb,kb,N)),
nobs =as.integer(nobs),
N =N,
M =as.integer(M),
ka =as.integer(ka),
kb =as.integer(kb),
p =as.integer(p),
lm =penalties[1],
ln =penalties[2],
lf =penalties[3],
lg =penalties[4],
K = Kt,
minV =control$minimum.variance,
k0 =control$binary.k0,
kr =control$binary.kr,
Iterations=control$max.iterations,
tol =control$convergence.tolerance,
dl =control$C.debug,
dp =double(dpl),
ip =integer(ipl)
),class=c('pfda.dual.bc.rawC','pfda.dual.bc','list'))}
} else {
structure(dual.cc.core(y,z,B,subject,ka,kb,lm,ln,lf,lg,K,min.v,max.I,tol)
,class=c('pfda.dual.bc.R','pfda.dual.bc','list'))
},list(subject=subject))
}
}
.dual.bc.updatePCS<-function(Yi, Zi, Da, Db, lambda, phii, psii, Y.rhoi, Z.rhoi, alpha, beta){
ni = NROW(phii)
# mean functions
muy = Y.rhoi+phii%*%alpha
mua = pnorm(muy)
Y.tilda = muy + (Yi-mua)/dnorm(mua)
Z.tilda = Zi
# blocks for means
A = tcrossprod(Da,phii) #' $ A = D_\alpha \phi_i^T $
B = Da %*% t(lambda) %*% t(psii) #' $ B = D_\alpha \Lambda^T \psi_i^T $
C = lambda %*% A #' $ C = \Lambda D_\alpha \psi_i^T = \Lambda A$
D = Db %*% t(psii) #' $ D = D_\beta \psi_i^T $
#blocks associated with variances
E = diag(1, ni) + (phii %*% A) #' $ E = I + \phi_i A = I + \phi_i D_\alpha \phi_i^T $
F = phii %*% B #' $ F = \phi_i D_\alpha \Lambda^T \psi_i $
G = psii %*% C #' $ G = \psi_i \Lambda D_\alpha \phi_i $
H = diag(1, ni) + psii %*% D #' $ H = I + \psi_i D_\beta \psi_i $
# inversion by block
E1 = solve(E)
T = solve(H-G%*%E1%*%F) #' $ T = (H - G E^{-1} F)^{-1} $
Q = E1+E1%*%F%*%T%*%G%*%E1 #' $ Q = E^{-1} + E^{-1} F T G E^{-1} $
R = -E1%*%F%*%T #' $ R = -E^{-1} F T $
S = -T%*%G%*%E1 #' $ S = -T G E^{-1} $
# results
list(
alpha= (A %*% Q + B %*% S) %*% Y.tilda + (A %*% R + B %*% T) %*% Z.tilda,
beta = (C %*% Q + D %*% S) %*% Y.tilda + (C %*% R + D %*% T) %*% Z.tilda)
}
.dual.bc.estimatePCS<-function(Y, Z, B, subject, tm, tn, tf, tg, Da, Db, lambda, sxi){
Y.rho = B%*%tm
Z.rho = B%*%tn
phi = B%*%tf
psi = B%*%tg
ofun<-function(i){
ix = as.integer(subject)==i
alpha = rep(0,NCOL(tf))
beta = rep(0,NCOL(tg))
dif = 1
while(dif>1e-3){
ab<-.dual.bc.updatePCS(Y[ix], Z[ix], Da, Db, lambda, phi[ix,,drop=F], psi[ix,,drop=F], Y.rho[ix,], Z.rho[ix,], alpha, beta)
dif<-mean(alpha-ab$alpha)+mean(beta - ab$beta)
alpha = ab$alpha
beta = ab$beta
}
list(alpha=alpha,beta=beta)
}
ab<-lapply(seq_len(nlevels(subject)), ofun)
alpha<-matrix(0, nlevels(subject), NCOL(tf))
beta <-matrix(0, nlevels(subject), NCOL(tg))
for(i in nlevels(subject)){
alpha[i,]<-ab[[i]]$alpha
beta[i,] <-ab[[i]]$beta
}
list(alpha=alpha,beta=beta)
}
.dual.bc.n2L<-function(y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, sxi){
ka = length(Da)
kb = length(Db)
N = nlevels(subject)
M = length(y)
phi = B %*% tf
psi = B %*% tg
pi = B %*% tm
rho = B %*% tn
ab = .dual.bc.estimatePCS(y, z, B, subject, tm, tn, tf, tg, Da, Db, lambda)
C = rbind(cbind(diag(Da,ka), diag(Da,ka)%*%t(lambda)),cbind(t(lambda%*%diag(Da,ka)), diag(Db,kb)))
N*(sum(log(Da))+determinant(diag(Db,kb)-lambda %*% tcrossprod(diag(Da,ka), lambda))$modulus) +
M*log(sxi) +
with(ab,sum(sapply(seq_len(N),function(i){
ix = subject==i
mu_alpha=pnorm( pi[ix]+phi[ix,,drop=F]%*%alpha[i,])
mu_beta =pnorm(rho[ix]+psi[ix,,drop=F]%*% beta[i,])
W_alpha.i = (dnorm(qnorm(mu_alpha))^2)/mu_alpha*(1-mu_alpha)
gamma<-c(alpha[i,],beta[i,])
# $ \log|I+\Omega W_i \Omega C| $
A<-diag(1,ka)+crossprod(phi[ix,,drop=F], (W_alpha.i * phi[ix,,drop=F]) %*% diag(Da,ka) )
determinant(A)$modulus +
determinant(diag(1,kb) + crossprod(psi[ix,,drop=F]) - lambda %*% (Da * solve(A)) %*% tcrossprod(diag(Da,ka),lambda) )$modulus +
# determinant(diag(1,ka)+ crossprod(phi[ix,],W_alpha.i*phi[ix,,drop=F])%*%C,log=T)$modulus+
# $\sumj^{n_i} Y_i $
sum(-2*(y[ix]*log(mu_alpha) + (1-y[ix])*log(1-mu_alpha)) + (z[ix]-mu_beta)^2) +
{ # $ \gamma_i^T C^{-1} \gamma_i $
Sib <- solve(diag(Db,kb) - lambda %*% Da %*% lambda)
Siab <- -crossprod(lambda ,Sib)
Sia <- diag(1/Da,ka) - Siab %*% lambda %*% diag(1/Da,ka)
crossprod(gamma,solve(C,gamma))
crossprod(alpha[i,],Sia %*% alpha[i,]) +
crossprod(alpha[i,],Siab %*% beta[i,])*2 +
crossprod(beta[i,] ,Sib %*% beta[i,])
}
})))
}
AIC.pfda.dual.bc<-function(object,...){
with(object,.dual.ge.AIC(
.dual.bc.n2L(y, z, B, subject, tm, tn, tf, tg, lambda, Da, Db, sxi),0,
B , B , B , B ,
1 , ka , 1 , kb,
lm , ln , lf , lg,
K , K , K , K ,
1, 1, sxi, sxi)
)
}
print.pfda.dual.bc<-function(x,...){
cat('Paired Functional Principal Component Model (Binary/Continuous)\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(attr(x,'name.y'),' has ', NCOL(x$tf),' principal components\n')
cat(attr(x,'name.z'),' has ', NCOL(x$tg),' principal components\n')
cat('penalties are \n');print( penalty.pfda.dual.bc(x))
}
penalty.pfda.dual.bc<-function(object,..)with(object,structure(matrix(c(lm,ln,lf,lg),2,2),dimnames=list(c(attr(object,'name.y'),attr(object,'name.z')),c('mean','pc'))))
}
{ # Additive aka Calcium Model
.dual.ca.i<-function(y,Z,Bt,Bx,subject,kg,kd,min.v){
tz<- if(!is.null(Z) && NCOL(Z)>=1) solve(crossprod(Z),crossprod(Z,y)) else numeric(0)
a<-.single.c.i(y,Bt,subject,kg,min.v)
b<-.single.c.i(y,Bx,subject,kd,min.v)
lambda <- crossprod(b$alpha,a$alpha)%*%solve(crossprod(a$alpha))
gd<-array(0,dim=c(kg,kd,nlevels(subject)))
R=y - Z%*%tz - Bt%*%a$tm - Bx%*%b$tm
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%a$tf%*%a$alpha[i,] - Bx[ix,]%*%b$tf%*%b$alpha[i,]
}
sigma = as.vector(crossprod(R)/length(R))
for(i in seq_len(nlevels(subject)))gd[,,i]<-tcrossprod(a$alpha[i,],b$alpha[i,])
{ list(
tz=tz,
tt=a$tm,tx=b$tm,
tf=a$tf,tg=b$tf,
gamma=a$alpha,delta=b$alpha,
Dg=a$Da,Dd=b$Da,
lambda=lambda,
sigma=sigma,
gg=array(unlist(a$aa),c(kg,kg,nlevels(subject))),
gd=gd,
dd=array(unlist(b$aa),c(kd,kd,nlevels(subject)))
)}
}
.dual.ca.E.1<-function(Ri,tt,tx,sigma, phi, psi, O){
# phi = Bti %*% tf
# psi = Bxi %*% tg
# O is list of block inverted Dg,C,Dd from .gen.symblock.solve
S = .gen.symblock.solve(
O[[1]]+crossprod(phi)/sigma,
O[[2]]+crossprod(phi,psi)/sigma,
O[[3]]+crossprod(psi)/sigma)
a=crossprod(phi,Ri/sigma)
b=crossprod(psi,Ri/sigma)
g=S[[1]]%*%a+S[[2]]%*%b
d=crossprod(S[[2]],a)+S[[3]]%*%b
list(
mu.gamma = as.vector(g),
mu.delta = as.vector(d),
Sigma=S,
gg = tcrossprod(g)+S[[1]],
gd = tcrossprod(g,d)+S[[2]],
dd = tcrossprod(d)+S[[3]]
)
}
.dual.ca.E<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma){
O<-.gen.symblock.solve(
diag(Dg,length(Dg)),
tcrossprod(diag(Dg,length(Dg)),lambda),
diag(Dd,length(Dd)))
{
kg<-NCOL(tf)
kd<-NCOL(tg)
n<-nlevels(subject)
Sgg<-array(0,dim=c(kg,kg,n))
Sgd<-array(0,dim=c(kg,kd,n))
Sdd<-array(0,dim=c(kd,kd,n))
gg <-array(0,dim=c(kg,kg,n))
gd <-array(0,dim=c(kg,kd,n))
dd <-array(0,dim=c(kd,kd,n))
gamma<-matrix(nrow=n,ncol=kg)
delta<-matrix(nrow=n,ncol=kd)
R<-y-Z%*%tz-Bt%*%tt-Bx%*%tx
phi = Bt%*%tf
psi = Bx%*%tg
}
for(i in seq_len(n)){
ix = i==as.integer(subject)
r<-.dual.ca.E.1(R[ix],tt,tx,sigma,phi[ix,,drop=F],psi[ix,,drop=F],O)
{
Sgg[,,i]<-r$Sigma[[1]]
Sgd[,,i]<-r$Sigma[[2]]
Sdd[,,i]<-r$Sigma[[3]]
gg[,,i]<-r$gg
gd[,,i]<-r$gd
dd[,,i]<-r$dd
gamma[i,]<-r$mu.gamma
delta[i,]<-r$mu.delta
}
}
list(Sgg=Sgg,Sgd=Sgd,Sdd=Sdd,gamma=gamma,delta=delta,gg=gg,gd=gd,dd=dd)
}
.dual.ca.unpenalized<-function(y,Z,Bt,Bx,subject,tt,tx,tf,tg,gamma,delta,sigma){
if(is.null(Z)||NCOL(Z)==0)return(list(tz=numeric(0)))
R=y-Bt%*%tt-Bx%*%tx
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%tf%*%gamma[i,] - Bx[ix,]%*%tg%*%delta[i,]
}
list(tz=solve(crossprod(Z),crossprod(Z,R)))
}
.dual.ca.penalized.1<-function(R, B, sigma, l, K){
solve(crossprod(B)+sigma*l*K,crossprod(B,R))
}
.dual.ca.penalized<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,gamma,delta,sigma,lt,lx,Kt,Kx){
R=y-Z%*%tz
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%tf%*%gamma[i,] - Bx[ix,]%*%tg%*%delta[i,]
}
list( tt= .dual.ca.penalized.1(R-Bx%*%tx,Bt,sigma,lt,Kt),
tx= .dual.ca.penalized.1(R-Bt%*%tt,Bx,sigma,lx,Kx))
}
.dual.ca.princcomp<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,gamma,delta,sigma,gg,gd,dd, lf,lg,Kt,Kx){
R = y - Z%*%tz -Bt%*%tt - Bx%*%tx
# tf
for(j in seq_len(NCOL(tf))){
left=sigma*lf*Kt
right=numeric(NROW(tf))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
left = left + gg[j,j,i]*crossprod(Bt[ix,])
right= right+ crossprod(
Bt[ix,,drop=F],
R[ix]*gamma[i,j]-
Bx[ix,,drop=F]%*%tg%*%gd[j,,i]-
(if(NCOL(tf)>1)Bt[ix,,drop=F]%*%tf[,-j]%*%gg[j,-j,i] else 0))
}
tf[,j]= .u.orthogonalize(tf[,seq_len(j-1),drop=FALSE],solve(left,right))
}
# tg
for(j in seq_len(NCOL(tg))){
left=sigma*lg*Kx
right=numeric(NROW(tg))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
left = left + dd[j,j,i]*crossprod(Bx[ix,])
right= right+ crossprod(
Bx[ix,,drop=F],
R[ix]*delta[i,j]-
Bt[ix,,drop=F]%*%tf%*%gd[,j,i] -
(if(NCOL(tg)>1)Bx[ix,,drop=F]%*%tg[,-j]%*%dd[j,-j,i] else 0))
}
tg[,j]= .u.orthogonalize(tg[,seq_len(j-1),drop=FALSE],solve(left,right))
}
# return value
list(tf = tf, tg = tg)
}
.dual.ca.variances<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,gamma,delta,sigma,gg,gd,dd){
R=y - Z%*%tz - Bt%*%tt - Bx%*%tx
sum.gg<-matrix(0,NCOL(tf),NCOL(tf))
sum.gd<-matrix(0,NCOL(tf),NCOL(tg))
sum.dd<-matrix(0,NCOL(tg),NCOL(tg))
for(i in seq_len(nlevels(subject))){
ix = i==as.integer(subject)
R[ix]=R[ix] - Bt[ix,]%*%tf%*%gamma[i,] - Bx[ix,]%*%tg%*%delta[i,]
sum.gg<-sum.gg+gg[,,i]
sum.gd<-sum.gd+gd[,,i]
sum.dd<-sum.dd+dd[,,i]
}
sigma = as.vector(crossprod(R)/length(R))
eg<-eigen(sum.gg/nlevels(subject))
ed<-eigen(sum.dd/nlevels(subject))
Tg<-eg$vectors*rep(sign(tf[1,]%*%eg$vectors),each=nrow(eg$vectors))
Td<-ed$vectors*rep(sign(tg[1,]%*%ed$vectors),each=nrow(ed$vectors))
lambda = t(solve(sum.gg,sum.gd))
list(sigma=sigma,
Dg = eg$values,
Dd = ed$values,
lambda = t(Td)%*%lambda%*%Tg,
tf = tf%*%Tg,
tg = tg%*%Td,
gamma = gamma%*%Tg,
delta = delta%*%Td
)
}
.dual.ca.core<-function(y,Z,Bt,Bx,subject,kg,kd,lt,lx,lf,lg,Kt,Kx,min.v,max.I,tol){
{ # initial values
if(is.null(Z)) Z<-matrix(nrow=NROW(y),ncol=0)
init<-.dual.ca.i(y,Z,Bt,Bx,subject,kg,kd,min.v)
init$I<-0
init$cctrace<-numeric(0)
}
within(init,while(I<-I+1){
E <- .dual.ca.E(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma)
U <- .dual.ca.unpenalized(y,Z,Bt,Bx,subject, tt, tx,tf,tg ,E$gamma,E$delta ,sigma)
P <- .dual.ca.penalized(y ,Z,Bt,Bx,subject,U$tz, tt, tx,tf,tg ,E$gamma,E$delta ,sigma,lt,lx,Kt,Kx)
C <- .dual.ca.princcomp(y ,Z,Bt,Bx,subject,U$tz,P$tt,P$tx,tf,tg ,E$gamma,E$delta ,sigma,E$gg,E$gd,E$dd, lf,lg,Kt,Kx)
V <- .dual.ca.variances(y ,Z,Bt,Bx,subject,U$tz,P$tt,P$tx,C$tf,C$tg,E$gamma,E$delta,sigma,E$gg,E$gd,E$dd)
{ #Convergence
ccl<-numeric(0)
ccl['tz']<-sum(abs((tz-U$tz)/tz))
ccl['tt']<-sum(abs((tt-P$tt)/tt))
ccl['tx']<-sum(abs((tx-P$tx)/tx))
ccl['tf']<-sum(abs((tf-V$tf)/tf))
ccl['tg']<-sum(abs((tg-V$tg)/tg))
ccl['Dg']<-sum(abs((Dg-V$Dg)/Dg))
ccl['Dd']<-sum(abs((Dd-V$Dd)/Dd))
(cc<-sum(ccl,na.rm=TRUE))
{ # reassign values
gg <-E$gg
gd <-E$gd
dd <-E$dd
tz <-U$tz
tt <-P$tt
tx <-P$tx
tf <-V$tf
tg <-V$tg
gamma <-V$gamma
delta <-V$delta
sigma <-V$sigma
Dg <-V$Dg
Dd <-V$Dd
lambda<-V$lambda
}
if(cc<tol)break
if(I>max.I){
warning('Maximum number of iterations exceeded, convergence not obtained.')
break
}
cctrace<-c(cctrace,cc)
}
})
}
.dual.ca.n2L<-function(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma,...){
ll<-0
phi = Bt%*%tf
psi = Bx%*%tg
Dg<-diag(Dg,length(Dg))
Dd<-diag(Dd,length(Dd))
Ry<-y-Bt%*%tt-Bx%*%tx-(if(is.null(Z))0 else Z%*%tz)
for(i in nlevels(subject)){
ix = i==as.integer(subject)
U = phi[ix,,drop=F]%*%tcrossprod(Dg,psi[ix,]%*%lambda)
S = phi[ix,,drop=F]%*%tcrossprod(Dg,phi[ix,,drop=F])+ psi[ix,,drop=F]%*%tcrossprod(Dd,psi[ix,,drop=F]) + U + t(U) + diag(sigma,sum(ix))
ll = ll+determinant(S)$modulus+crossprod(Ry[ix],solve(S)%*%Ry[ix])
}
attr(ll,"log")<-TRUE
ll
}
AIC.pfda.additive<-function(object,...){
with(object,.dual.ge.AIC(.dual.ca.n2L(y,Z,Bt,Bx,subject, tz, tt, tx, tf, tg, lambda, Dg, Dd, sigma),
kz,
Bt, Bt, Bx, Bx,
1 , kg, 1 , kd,
lt, lf, lx, lg,
Kt, Kt, Kx, Kx,
sigma, sigma, sigma, sigma))
# with(object,.dual.ca.AIC.rawC(object, y, Z, Bt, Bx, subject, Kt, Kx))
}
logLik.pfda.additive<-function(object,...,newdata=NULL,n2L=TRUE){
with(object,with(newdata,.dual.ca.n2L(y,Z,Bt,Bx,subject,tz,tt,tx,tf,tg,lambda,Dg,Dd,sigma,...)))
}
dual.ca<-function(y,Z,t,x,subject,knots=NULL,penalties=NULL,df=NULL,k=NULL,bases=NULL,control=pfdaControl()){
# name.t = deparse(substitute(t))
# name.x = deparse(substitute(x))
{ #setup material
Z = if(is.null(Z))matrix(nrow=length(y),ncol=0) else as.matrix(Z)
if(is.null(bases))bases=new.env()
stopifnot(is.environment(bases))
if(!(exists("tbase",envir=bases,inherits=FALSE) && exists("xbase",envir=bases, inherits=FALSE))){ # knot and base identification
if(is.null(knots)){
if(length(control$nknots)==2) { nkt <- control$nknots[1]; nkx <- control$nknots[2] } else nkt <- nkx <- control$nknots[1]
kt<-expand.knots(unique(quantile(t,seq(0,1,length.out=nkt))))
kx<-expand.knots(unique(quantile(x,seq(0,1,length.out=nkt))))
} else if(length(knots)==2) {
if(name.t %in% names(knots)) kt <-knots[[name.t]]
else if('t' %in% names(knots)) kt<-knots[['t']]
else kt <- knots[[1]]
if(name.x %in% names(knots)) kx <-knots[[name.x]]
else if('x' %in% names(knots)) kx<-knots[['x']]
else kx <- knots[[2]]
} else stop("knots needs to be a list of length 2")
assign('tbase',OBasis(kt),envir=bases)
assign('xbase',OBasis(kx),envir=bases)
knots<-list(kt,kx)
# names(knots)<-c(name.t,name.x)
}
{ # evaluate Bases
Bt = evaluate(bases$tbase,t)
Bx = evaluate(bases$xbase,x)
Kt = OuterProdSecondDerivative(bases$tbase)
Kx = OuterProdSecondDerivative(bases$xbase)
Bx<-Bx[,-1] #Removing the first column of Bx is to improve stability and remove interdependence between t and x
Kx<-Kx[-1,-1]
}
eval(.X.dual.k)
eval(.X.dual.penalties)
}
if(any(is.na(k))){
# stop("identification of number of principle components is not done yet.")
k<-c(1,1)
model0<-RecallWith(k=k)
AIC0<-AIC(model0)
while(TRUE){
modelA<-try(RecallWith(k=c(model0$kg+1,model0$kd)),silent=TRUE)
modelB<-try(RecallWith(k=c(model0$kg,model0$kd+1)),silent=TRUE)
AICA<- if(is(modelA,'pfda.additive')) AIC(modelA) else Inf
AICB<- if(is(modelB,'pfda.additive')) AIC(modelB) else Inf
if(AIC0<AICA && AIC0<AICB) break
else if(AICA<AICB){
model0<-modelA
AIC0<-AICA
} else {
model0<-modelB
AIC0<-AICB
}
}
return(model0)
} else
if (any(is.na(penalties))) {
funcall <- structure(match.call(),envir=parent.frame(),expand.dots = FALSE)
if(is.null(control$optim.start)){
lt <- l.from.df(2.1,Bt,Kt)
lx <- l.from.df(2.1,Bx,Kx)
control$optim.start<-c(lt,lx,lt,lx)
}
eval(.X.optimize.penalties) }
else {
rtn<-if(control$useC){
{ # compute memory and parameters
kz = if(is.null(Z))0L else NCOL(Z)
kg = as.integer(k[1])
kd = as.integer(k[2])
N = nlevels(subject)
nobs= table(subject)
M = length(y)
pt = ncol(Bt)
px = ncol(Bx)
p = max(pt, px)
k = max(kg, kd)
dpl = kz + kg + kd + pt +px + pt*kg + px*kd + kg*kd + max(
pt*pt*N + px*px*N + 2*p^2 + M + M*k+ p*N + 8*p,
2^kg^2 + 2*kd^2 + kg*kd+ max(kg,kd)*10 + 2*kg*kd+(1+kg+kd)*M,
M + kz^2 + 10*kz + M*max(kg, kd),
10*p + p^2 +M,
M + p^2 + p + max( 10*pt , M + N*(kg+kd) + M*max(kg,kd) ),
2*kg^2 + 2*kd^2 + M + 10*k + 9*p + 3*p^2 + N*k + M*k
)
}
structure(.C('dual_ca_core', residuals=y, Z=Z, Bt=Bt, Bx=Bx,
tz=double(kz), tt=double(pt), tx=double(px),
tf=matrix(0,pt,kg), tg=matrix(0,px,kd),
gamma=matrix(0,N,kg), delta=matrix(0,N,kd),
lambda=matrix(0,kd,kg), Dg=double(kg), Dd=double(kd),
sigma=0.0,
gg=array(0,c(kg,kg,N)), gd=array(0,c(kg,kd,N)), dd=array(0,c(kd,kd,N)),
Sgg=array(0,c(kg,kg,N)), Sgd=array(0,c(kg,kd,N)), Sdd=array(0,c(kd,kd,N)),
nobs=nobs, N=N, M=M, kz=kz, kg=kg, kd=kd, pt=pt, px=px,
lt=penalties[1], lx=penalties[2], lf=penalties[3], lg=penalties[4], Kt=Kt, Kx=Kx,
minV=control$minimum.variance, Iterations=control$max.iterations, tol=control$convergence.tolerance,
dl=control$C.debug, dp=double(dpl), ip=integer(max(6*p,kz)))
,class=c('pfda.additive.rawC','pfda.additive','list'))
} else {
structure(.dual.ca.core(y,Z,Bt,Bx,subject,k[1],k[2],penalties[1],penalties[2],penalties[3],penalties[4],Kt,Kx,control$minimum.variance,control$max.iterations,control$convergence.tolerance)
,class=c('pfda.additive.R','pfda.additive','list'))
}
rtn$tbase<-bases$tbase
rtn$xbase<-bases$xbase
rtn$y<-y
rtn$subject<-subject
rtn
}
}
plot.pfda.additive<-function(x,...){
with(x,{
layout(matrix(1:4,nrow=2,ncol=2,byrow=T))
plot(tbase,tt, main=paste("Plot of mean curve for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab=attr(x,'name.y'))
plot(tbase,tf, main=paste("Principle components for",attr(x,'name.t')),xlab=attr(x,'name.t'),ylab='')
plot(xbase,c(0,tx), main=paste("Plot of mean curve for ",attr(x,'name.x')),xlab=attr(x,'name.x'),ylab=attr(x,'name.y'))
plot(xbase,rbind(0,tg), main=paste("Principle components for ",attr(x,'name.x')),xlab=attr(x,'name.x'),ylab='')
})
}
print.pfda.additive<-function(x,...){
cat('Additive Principal Component Model\n')
cat('Formula: ', deparse(attr(x,'formula')),'\n')
cat(attr(x,'name.t'),' has ', NCOL(x$tf),' principal components\n')
cat(attr(x,'name.x'),' has ', NCOL(x$tg),' principal components\n')
cat('penalties are \n');print( penalty.pfda.additive(x))
}
penalty.pfda.additive<-function(object,..)with(object,structure(matrix(c(lt,lx,lf,lg),2,2),dimnames=list(c(attr(object,'name.t'),attr(object,'name.x')),c('mean','pc'))))
persp.pfda.additive<-function(x,col.fun=NULL,nt=11,nx=11,...){
if(is.null(col.fun))col.fun<-heat.colors
with(x,{
tx.grid<-expand.grid(
t=seq(tbase@knots[4],tbase@knots[length(tbase@knots)-3],length=nt),
x=seq(xbase@knots[4],xbase@knots[length(xbase@knots)-3],length=nx))
Btg<-evaluate(tbase,tx.grid$t)
Bxg<-evaluate(xbase,tx.grid$x)[,-1]
yg<-matrix(Btg%*%(tt+tf%*%apply(gamma,2,mean)) + Bxg%*%(tx+tg%*%apply(delta,2,mean)),nx,nt,byrow=T)
{
# this code was copied and modified from the persp example page
z <- yg
nrz <- nrow(yg)
ncz <- ncol(yg)
# Create a function interpolating colors in the range of specified colors
nbcol <- 200
color <- col.fun(nbcol)
zfacet <- yg[-1, -1] + yg[-1, -ncz] + yg[-nrz, -1] + yg[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)
}
persp(yg,zlab=attr(x,"name.y"),xlab=attr(x,"name.x"),ylab=attr(x,"name.t"),
col=color[facetcol], ...)
})
}
}
{ # general
.dual.ge.AIC<-function(n2L,kz,
B1, B2, B3, B4,
k1, k2, k3, k4,
l1, l2, l3, l4,
K1, K2, K3, K4,
s1, s2, s3, s4){
as.vector(n2L)+2*(kz +
k1*.pfda.df(B1,l1,K1,s1) +
k2*.pfda.df(B2,l2,K2,s2) +
k3*.pfda.df(B3,l3,K3,s3) +
k4*.pfda.df(B4,l4,K4,s4)
)
}
pfda<-function(model, data=environment(model), ..., driver){
mf <- pfdaParseFormula(model,data)
if(missing(driver))driver = infer.driver(mf)
structure(with(mf,switch(driver,
single.continuous = structure(single.c(response,additive,splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response,'name'),name.t=attr(splinegroup[[1]],'name')),
single.binary = structure(single.b(response,splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response,'name'),name.t=attr(splinegroup[[1]],'name')),
dual.continuous = structure(dual.cc(response[[1]],response[[2]],splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response[[1]],'name'),name.z=attr(response[[2]],'name'),name.t=attr(splinegroup[[1]],'name')),
dual.mixed = structure(dual.bc(response[[1]],response[[2]],splinegroup[[1]],splinegroup[[2]],...),name.y=attr(response[[1]],'name'),name.z=attr(response[[2]],'name'),name.t=attr(splinegroup[[1]],'name'),name.x=attr(splinegroup[[1]],'name')),
additive = structure(dual.ca(response,additive,splinegroup[[1]],splinegroup[[2]],splinegroup[[3]],...),name.y=attr(response,'name'),name.t=attr(splinegroup[[1]],'name'),name.x=attr(splinegroup[[2]],'name'))
)), formula = model)
}
}
Try the pfda 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.