Nothing
fdaLm <- function(formula,data,design,
boxcox=NULL,G=1,lambda=1,nlSearch=TRUE,
K.order=1,D.order=NULL,Fleft="tied",Fright="tied",
left=NULL,right=NULL
# ,lag.max=25
){
# --------------------------------------
# Take call
# --------------------------------------
mf.call <- match.call(expand.dots=FALSE)
# --------------------------------------
# Extract observation and identifier
# --------------------------------------
# Extract observation vector
mf <- mf.call[c(1,match(c("formula","data"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=1,rhs=0)
mf$drop.unused.levels <- TRUE
Ymat <- eval(mf,parent.frame())[,1]
Ntotal <- length(Ymat)
# Extract sample identifier
if (length(Formula::Formula(formula))[1]==1) {
# Only one sample
MM <- 1
} else {
# Look for sample identifier inside 'data'
mf <- mf.call[c(1,match(c("formula","data"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=2,rhs=0)
mf$drop.unused.levels <- TRUE
id.data <- try(as.factor(eval(mf,parent.frame())[,1]),silent=TRUE)
if (!inherits(id.data,"try-error")) {
id.names.data <- levels(id.data)[unique(id.data)]
MM.data <- length(id.names.data)
}
# Look for sample identifier inside 'design'
mf <- mf.call[c(1,match(c("formula","design"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=2,rhs=0)
mf$drop.unused.levels <- TRUE
names(mf)[match("design",names(mf),0)] <- "data"
id.design <- try(as.factor(eval(mf,parent.frame())[,1]),silent=TRUE)
if (!inherits(id.design,"try-error")) {
id.names.design <- levels(id.design)[unique(id.design)]
MM.design <- length(id.names.design)
}
# Validate identifiers
if (inherits(id.data,"try-error")) {
if (inherits(id.design,"try-error")) {
stop("Sample identifier not found")
} else {
MM <- MM.design
id.names <- id.names.design
}
} else {
if (inherits(id.design,"try-error")) {
MM <- MM.data
id.names <- id.names.data
} else {
if (MM.data!=MM.design) stop("Mismatch between length of sample identifiers inside 'data' and 'design'")
if (any(is.na(match(id.names.data,id.names.design)))) stop("Mismatch between names of sample identifiers inside 'data' and 'design'")
MM <- MM.data
id.names <- id.names.data
}
}
}
# ------------------------------------
# Validate observation and identifier
# ------------------------------------
NN <- Ntotal %/% MM
if ((Ntotal %% MM)!=0) stop("Observation length not a multiple of sample numbers")
if (NN < 3) stop("Sample length must be at least 3")
# Copy orignal if needed in Box-Cox analysis
if (is.numeric(boxcox)) {
if (!all(Ymat > 0)) stop("In Box-Cox analysis observations must be positive")
geometricMean <- exp(mean(log(Ymat)))
Yoriginal <- matrix(Ymat,NN,MM) # matrix() to make actual copy
}
# Reshape observation vector as a matrix
attr(Ymat,"dim") <- c(NN,MM)
# ------------------------------------------
# Make design matrices for fixed effect
# ------------------------------------------
# Checks whether all variables are inside 'design' data frame
# Might be invoked later!!!
#if (is.element("design",names(mf.call))) {
# mf <- mf.call[c(1,match("design",names(mf.call)))]
# mf[[1]] <- as.name("colnames")
# names(mf)[2] <- "x"
# if (all(is.element(as.character(attr(terms(formula(Formula(formula),lhs=0,rhs=1)),"variables"))[-1],eval(mf,parent.frame())))) warning("All fixed effect design variables found in 'design' data frame")
#}
# Look for design variables inside 'design'
look.in.data <- FALSE
if (is.element("design",names(mf.call))) {
mf <- mf.call[c(1,match(c("formula","design"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=0,rhs=1)
mf$drop.unused.levels <- TRUE
names(mf)[match("design",names(mf),0)] <- "data"
mf <- try(suppressWarnings(eval(mf,parent.frame())),silent=TRUE)
if (inherits(mf,"try-error")) {
warning("All fixed effect design variables not found in 'design' and/or parent environment. I will look in 'data' frame")
look.in.data <- TRUE
}
} else {
look.in.data <- TRUE
}
if (look.in.data) {
# Not all design variables found inside 'design'
# Look for design variables inside 'data'
mf <- mf.call[c(1,match(c("formula","data"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=0,rhs=1)
mf$drop.unused.levels <- TRUE
mf <- try(suppressWarnings(eval(mf,parent.frame())),silent=TRUE)
if (inherits(mf,"try-error")) stop("All fixed effect design variables not found in 'data' and/or parent environment")
}
# If identifiers in both 'data' and 'design' with same length as data frame,
# then reorder data frame
if ((MM > 1) && (!inherits(id.data,"try-error")) && (!inherits(id.design,"try-error")) && (dim(mf)[1]==MM)) {
mf <- mf[match(id.names.data,id.names.design),,drop=FALSE]
}
# Make design matrix
GammaMat <- model.matrix(Formula::Formula(formula),data=mf,lhs=0,rhs=1)
beta.names <- colnames(GammaMat)
p0 <- length(beta.names)
Nfixed <- dim(GammaMat)[1]
attributes(GammaMat) <- NULL
attr(GammaMat,"dim") <- c(Nfixed,p0)
# ------------------------------------------
# Make design matrices for random effect
# ------------------------------------------
# Design matrix for random effects
if (length(Formula::Formula(formula))[2]==1) {
# No random effects
q0 <- 0
Zmat <- matrix(0,MM,q0)
} else {
# Look for design variables inside 'design'
look.in.data <- FALSE
if (is.element("design",names(mf.call))) {
mf <- mf.call[c(1,match(c("formula","design"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=0,rhs=2)
mf$drop.unused.levels <- TRUE
names(mf)[match("design",names(mf),0)] <- "data"
mf <- try(suppressWarnings(eval(mf,parent.frame())),silent=TRUE)
if (inherits(mf,"try-error")) {
warning("All fixed random design variables not found in 'design' and/or parent environment. I will look in 'data' frame")
look.in.data <- TRUE
}
} else {
look.in.data <- TRUE
}
if (look.in.data) {
# Not all design variables found inside 'design'
# Look for design variables inside 'data'
mf <- mf.call[c(1,match(c("formula","data"),names(mf.call),0))]
mf[[1]] <- as.name("model.frame")
mf$formula <- formula(Formula::Formula(formula),lhs=0,rhs=2)
mf$drop.unused.levels <- TRUE
mf <- try(suppressWarnings(eval(mf,parent.frame())),silent=TRUE)
if (inherits(mf,"try-error")) stop("All random effect design variables not found in 'data' and/or parent frame")
}
# If identifiers in both 'data' and 'design' with same length as data frame,
# then reorder data frame
if ((MM > 1) && (!inherits(id.data,"try-error")) && (!inherits(id.design,"try-error")) && (dim(mf)[1]==MM)) {
mf <- mf[match(id.names.data,id.names.design),,drop=FALSE]
}
# Make design matrix
Zmat <- model.matrix(Formula::Formula(formula),data=mf,lhs=0,rhs=2)
u.names <- colnames(Zmat)
q0 <- length(u.names)
Nrandom <- dim(Zmat)[1]
attributes(Zmat) <- NULL
attr(Zmat,"dim") <- c(Nrandom,q0)
}
# ---------------------------------------------
# Interpretation of intercepts if necessary
# ---------------------------------------------
if ((p0==1) && (beta.names=="(Intercept)")) {
if ((q0==0) || ((q0==1) && (u.names=="(Intercept)"))) {
# Interpret fixed effect as global intercept
Nfixed <- Ntotal
GammaMat <- matrix(1,Nfixed,1)
} else {
# Interpret fixed effect as the random effect
Nfixed <- Nrandom
GammaMat <- matrix(1,Nfixed,1)
}
}
if ((q0==1) && (u.names=="(Intercept)")) {
if ((p0==0) || ((p0==1) && (beta.names=="(Intercept)"))) {
# Interpret random effect as global intercept
Nrandom <- Ntotal
Zmat <- matrix(1,Nrandom,1)
} else {
# Interpret fixed effect as the random effect
Nrandom <- Nfixed
Zmat <- matrix(1,Nrandom,1)
}
}
# ----------------------------------------------
# Validate dimensions and resolve ANOVA's
# ----------------------------------------------
if ((p0>0) && (q0>0) && (Nfixed!=Nrandom)) stop("Mismatch between rows in design matrices for fixed and random effects")
if ((p0>0) && (!is.element(Nfixed,c(0,MM,Ntotal)))) stop("Mismatch between sample identifier and design frame for fixed effect")
if ((q0>0) && (!is.element(Nrandom,c(0,MM,Ntotal)))) stop("Mismatch between sample identifier and design frame for random effect")
anova.type <- "global"
if ((p0>0) && (Nfixed==MM)) anova.type <- "marginal"
if ((q0>0) && (Nrandom==MM)) anova.type <- "marginal"
# ----------------------------
# Interpret parameters
# ----------------------------
# Interpret lambda parameter
lambda[lambda<0] <- 0
coefsMat <- matrix(0,sum(lambda>0),3)
if (lambda[1]==0) stop("Leading lambda coefficient must be strictly positive")
if (K.order %% 2 == 1) {
if (length(lambda)>2) stop("For odd K.order lambda must have length 1 or 2")
lambda.name <- paste("-D^",2*K.order,sep="")
coefsMat[1,3] <- -1
if (length(lambda)>1) {
lambda.name <- c(lambda.name,"I")
if (lambda[2]>0) coefsMat[2,1] <- 1
}
} else {
if (length(lambda) > 3) stop("For even K.order lambda must have length 1, 2 or 3")
lambda.name <- paste("D^",2*K.order,sep="")
coefsMat[1,3] <- 1
if (length(lambda)>1) {
if (K.order %% 4 == 2) {
lambda.name <- c(lambda.name,paste("-D^",K.order,sep=""))
if (lambda[2]>0) coefsMat[2,2] <- -1
} else {
lambda.name <- c(lambda.name,paste("D^",K.order,sep=""))
if (lambda[2]>0) coefsMat[2,2] <- 1
}
}
if (length(lambda)>2) {
lambda.name <- c(lambda.name,"I")
if (lambda[3]>0) coefsMat[2+(lambda[2]>0),1] <- 1
}
}
# Initialize D.order is necessary
if (is.null(D.order)) {
D.order <- K.order
} else {
if ((nlSearch) && (D.order < K.order)) stop("Parameter search requested with D.order < K.order")
}
# Initialize boundary conditions
if (is.character(Fleft) && (Fleft=="tied")) Fleft <- cbind(diag(K.order),matrix(0,K.order,K.order))
if (is.character(Fleft) && (Fleft=="open")) Fleft <- cbind(matrix(0,K.order,1),diag(K.order),matrix(0,K.order,K.order-1))
if (!is.matrix(Fleft)) stop("Fleft misspecified")
if (any(dim(Fleft)!=c(K.order,2*K.order))) stop(paste("Matrix of left boundary conditions does not have dimension=(",K.order,",",2*K.order,")",sep=""))
if (is.character(Fright) && (Fright=="tied")) Fright <- cbind(diag(K.order),matrix(0,K.order,K.order))
if (is.character(Fright) && (Fright=="open")) Fright <- cbind(matrix(0,K.order,1),diag(K.order),matrix(0,K.order,K.order-1))
if (!is.matrix(Fright)) stop("Fright misspecified")
if (any(dim(Fright)!=c(K.order,2*K.order))) stop(paste("Matrix of right boundary conditions does not have dimension=(",K.order,",",2*K.order,")",sep=""))
# Initialize limits if necessary
if (is.null(left)) left <- 0
if (is.null(right)) right <- NN
if (left >= right) stop("Left limit must be strictly less than right limit")
# ----------------------------
# Initialize variables
# ----------------------------
# Initialize variables to contain results
if (anova.type=="global") {
attr(GammaMat,"dim") <- c(NN,MM*p0)
attr(Zmat,"dim") <- c(NN,MM*q0)
Ymat <- cbind(Ymat,GammaMat,Zmat)
projMat <- matrix(0,NN,MM*(1+p0+q0)*(1+D.order))
betaHat <- rep(0,p0)
uBLUP <- rep(0,q0)
}
if (anova.type=="marginal") {
Yres <- matrix(0,NN,MM)
projMat <- matrix(0,NN,MM*(1+D.order))
betaHat <- matrix(0,NN,p0)
uBLUP <- matrix(0,NN,q0)
uBLUPinvG <- matrix(0,NN,q0)
proj.uBLUPinvG <- matrix(0,NN,q0)
condRes.uBLUPinvG <- rep(0,NN*q0)
xBLUP.uBLUPinvG <- matrix(0,NN*q0,1)
logLik.uBLUPinvG <- rep(0,5)
dummyVec <- rep(0,0)
dummyMat <- matrix(0,0,0)
}
condRes <- rep(0,NN*MM)
xBLUP <- matrix(0,NN*MM,1+D.order)
logLikVec <- rep(0,5+D.order)
Cbeta <- matrix(0,p0,p0)
Cu <- matrix(0,q0,q0)
# ---------------------------
# Likelihood function
# ---------------------------
minus2logLik <- function(param,return.sigma=FALSE) {
# -------------------------------------------------------------------------
# Extract parameters:
# Present implementation of variance parameters, might be improved later:
# 1) If q0!=0 then G is assumed to be proportional to identity matrix
# 2) Parameters modeled on log scale
# -------------------------------------------------------------------------
if (is.numeric(boxcox)) {
mu <- param[1]
if (q0 > 0) {
G <- exp(param[2])*diag(q0)
coefs <- as.vector(exp(param[-c(1,2)])%*%coefsMat)
} else {
G <- diag(0)
coefs <- as.vector(exp(param[-1])%*%coefsMat)
}
} else {
if (q0 > 0) {
G <- exp(param[1])*diag(q0)
coefs <- as.vector(exp(param[-1])%*%coefsMat)
} else {
G <- diag(0)
coefs <- as.vector(exp(param)%*%coefsMat)
}
}
# -----------------------------
# Box-Cox transformation
# -----------------------------
if (is.numeric(boxcox)) {
# Do Box-Cox transformation on untransformed observations
.Call("boxcoxTransform",
mu,
geometricMean,
Yoriginal,
Ymat,
PACKAGE="fdaMixed")
}
# --------------------
# Find eigenvalues
# --------------------
roots <- findRoots(c(1,0,0)+(right-left)/NN*coefs,K.order)
# ---------------------
# Trace integral
# ---------------------
innerIntegral <- function(u) {
# Find eigenvalues
roots <- findRoots(cbind(u+(right-left)/NN*coefs[1],
(right-left)/NN*coefs[2],
(right-left)/NN*coefs[3]),K.order)
# Make call
return(as.vector(.Call("fdaTrace",
left,right,coefs[3],
Re(roots$left),Im(roots$left),
Re(roots$right),Im(roots$right),
Fleft,Fright,NN,
PACKAGE="fdaMixed")))
}
# --------------------------------------------
# Marginal ANOVA (multivariate calibration)
# --------------------------------------------
if (anova.type=="marginal") {
# Call marginal ANOVA's
.Call("marginalANOVA",
G,
Ymat,
GammaMat,Zmat,
Yres,
betaHat,uBLUP,uBLUPinvG,
Cbeta,Cu,
PACKAGE="fdaMixed")
# Compute squared length of projected random effects
if (q0==0) {
quad.uBLUP <- 0
} else {
# Call FDA projection
.Call("fdaEngine",
left,right,(right-left)/NN*coefs[3],
Re(roots$left),Im(roots$left),Re(roots$right),Im(roots$right),
Fleft,Fright,
G,
uBLUPinvG,
proj.uBLUPinvG,dummyVec,dummyVec,xBLUP.uBLUPinvG,condRes.uBLUPinvG,
dummyMat,dummyMat,logLik.uBLUPinvG,
PACKAGE="fdaMixed")
# Second quadratic term
quad.uBLUP <- sum(as.vector(uBLUP)*condRes.uBLUPinvG)
}
# ----------------------------------------------
# Fit serial correlated effects
# ----------------------------------------------
# Make call
.Call("fdaEngine",
left,right,(right-left)/NN*coefs[3],
Re(roots$left),Im(roots$left),Re(roots$right),Im(roots$right),
Fleft,Fright,
G,
Yres,
projMat,dummyVec,dummyVec,xBLUP,condRes,dummyMat,dummyMat,logLikVec,
PACKAGE="fdaMixed")
# ---------------------------
# Estimate noise variance
# ---------------------------
if (K.order %% 2 == 1) {
sigma2hat <- sum(c(logLikVec[3],quad.uBLUP,
coefs[1]*logLikVec[5],
-coefs[3]*logLikVec[5+K.order])
)/(NN*(MM-p0))
} else {
sigma2hat <- sum(c(logLikVec[3],quad.uBLUP,
coefs[1]*logLikVec[5],
-coefs[2]*logLikVec[5+K.order/2],
coefs[3]*logLikVec[5+K.order])
)/(NN*(MM-p0))
}
# ---------------------------
# Compute -2*log(likelihood)
# ---------------------------
res <- NN*(MM-p0)*log(sigma2hat)+
(MM-p0)*integrate(innerIntegral,0,1,stop.on.error=FALSE)$value+
NN*(MM-p0)
if (q0!=0) res <- res + NN*determinant(diag(q0)+t(Zmat)%*%Zmat%*%G)$modulus
if (p0!=0) res <- res - NN*determinant(Cbeta)$modulus
attributes(res) <- NULL
}
# ----------------------------------------------
# Global ANOVA: Fit all effects simultaneously
# ----------------------------------------------
if (anova.type=="global") {
# Make call
.Call("fdaEngine",
left,right,(right-left)/NN*coefs[3],
Re(roots$left),Im(roots$left),Re(roots$right),Im(roots$right),
Fleft,Fright,
G,
Ymat,
projMat,betaHat,uBLUP,xBLUP,condRes,Cbeta,Cu,logLikVec,
PACKAGE="fdaMixed")
# ---------------------------
# Estimate noise variance
# ---------------------------
if (K.order %% 2 == 1) {
sigma2hat <- sum(c(logLikVec[3:4],coefs[1]*logLikVec[5],
-coefs[3]*logLikVec[5+K.order])
)/(NN*MM-p0)
} else {
sigma2hat <- sum(c(logLikVec[3:4],coefs[1]*logLikVec[5],
-coefs[2]*logLikVec[5+K.order/2],
coefs[3]*logLikVec[5+K.order])
)/(NN*MM-p0)
}
# ---------------------------
# Compute -2*log(likelihood)
# ---------------------------
res <- (NN*MM-p0)*log(sigma2hat)+
MM*integrate(innerIntegral,0,1,stop.on.error=FALSE)$value+
logLikVec[1]+logLikVec[2]+NN*MM-p0
}
# -------------------
# Return result
# -------------------
if (return.sigma) {
return(list(logLik=res,sigma2hat=sigma2hat))
} else {
return(res)
}
# End objective function: -2*log(likelihood)
}
# ----------------------------------------
# Do non-linear optimization if requested
# ----------------------------------------
if (q0 > 0) {
if (is.numeric(boxcox)) {
param <- c(boxcox[1],log(G),log(lambda[lambda>0]))
if (nlSearch) {
res <- nlminb(param,minus2logLik)
counts <- 1+sum(c(1,length(param))*res$evaluations)
param <- res$par
} else {
counts <- 1
}
boxcox <- param[1]
G <- exp(param[2])*diag(q0)
lambda[lambda>0] <- exp(param[-c(1,2)])
} else {
param <- c(log(G),log(lambda[lambda>0]))
if (nlSearch) {
res <- nlminb(param,minus2logLik)
counts <- 1+sum(c(1,length(param))*res$evaluations)
param <- res$par
} else {
counts <- 1
}
boxcox <- "not done"
G <- exp(param[1])*diag(q0)
lambda[lambda>0] <- exp(param[-1])
}
} else {
if (is.numeric(boxcox)) {
param <- c(boxcox[1],log(lambda[lambda>0]))
if (nlSearch) {
res <- nlminb(param,minus2logLik)
counts <- 1+sum(c(1,length(param))*res$evaluations)
param <- res$par
} else {
counts <- 1
}
boxcox <- param[1]
G <- diag(q0)
lambda[lambda>0] <- exp(param[-1])
} else {
param <- log(lambda[lambda>0])
if (nlSearch) {
res <- nlminb(param,minus2logLik)
counts <- 1+sum(c(1,length(param))*res$evaluations)
param <- res$par
} else {
counts <- 1
}
boxcox <- "not done"
G <- diag(q0)
lambda[lambda>0] <- exp(param)
}
}
# Extract L-coefficients and estimate error variance
# NB: Variable 'res' is used in final return, and hence should not be distorted
res <- minus2logLik(param,return.sigma=TRUE)
# -----------------------------------------
# Compute conditional variance
# (experimental code, might be improved)
# Non-correct formula used below. Hence excluded in Version 0.1. To be updated later.
# -----------------------------------------
# # Make coefficients
# coefs <- lambda[lambda>0]%*%coefsMat
#
# # Find variance of conditional residuals
# if (anova.type=="global") {
# # Diagonal of Green function corresponding to R*solve(I+R)
# invAdiag <- rep(0,NN)
# acfVec <- rep(0,min(1+lag.max,NN-2))
# roots <- findRoots(c(1,0,0)+(right-left)/NN*coefs,K.order)
# .Call("diagGreen",
# left,right,(right-left)/NN*coefs[3],
# Re(roots$left),Im(roots$left),Re(roots$right),Im(roots$right),
# Fleft,Fright,
# invAdiag,acfVec,
# PACKAGE="fdaMixed")
#
# # Variance of conditional residuals
# condResVar <- outer(invAdiag,rep(res$sigma2hat,MM))
# acfVec <- acfVec/acfVec[1]
#
# warning("Variance of conditional residuals not fully implemented for global ANOVA")
# warning("Auto covariance function might be wrong for global ANOVA")
# }
#
# if (anova.type=="marginal") {
# # Diagonal of Green function corresponding to R*solve(I+R)
# invAdiag <- rep(0,NN)
# acfVec <- rep(0,min(1+lag.max,NN-2))
# roots <- findRoots(c(1,0,0)+(right-left)/NN*coefs,K.order)
# .Call("diagGreen",
# left,right,(right-left)/NN*coefs[3],
# Re(roots$left),Im(roots$left),Re(roots$right),Im(roots$right),
# Fleft,Fright,
# invAdiag,acfVec,
# PACKAGE="fdaMixed")
#
# # Make matrix terms
# if (q0==0) {
# mat1 <- diag(MM)
# mat2 <- diag(MM)
# } else {
# mat1 <- diag(MM)-Zmat%*%Cu%*%t(Zmat)
# mat2 <- diag(MM)+Zmat%*%G%*%t(Zmat)
# }
# if (p0==0) {
# mat3 <- diag(MM)
# } else {
# mat3 <- GammaMat%*%Cbeta%*%t(GammaMat)
# }
#
# # Variance of conditional residuals
# condResVar <- res$sigma2hat*acfVec[1]*
# (mat1%*%mat2%*%mat1
# -mat1%*%mat2%*%mat1%*%mat3%*%mat1
# -mat1%*%mat3%*%mat1%*%mat2%*%mat1
# +mat1%*%mat3%*%mat1%*%mat2%*%mat1%*%mat3%*%mat1)
# if (MM > 0) {
# colnames(condResVar) <- id.names
# rownames(condResVar) <- id.names
# }
#
# # Rescale auto covariance as auto correlation
# invAdiag <- invAdiag/acfVec[1]
# acfVec <- acfVec/acfVec[1]
# }
# --------------------
# Reshape results
# --------------------
# Reshape xBLUP and conditional residuals
attr(xBLUP,"dim") <- c(NN,MM,1+D.order)
attr(condRes,"dim") <- c(NN,MM)
# Add appropriate name attributes
names(lambda) <- lambda.name
if (p0 > 0) {
if (anova.type=="global") {
names(betaHat) <- beta.names
} else {
colnames(betaHat) <- beta.names
}
rownames(Cbeta) <- beta.names
colnames(Cbeta) <- beta.names
}
if (q0 > 0) {
if (anova.type=="global") {
names(uBLUP) <- u.names
} else {
colnames(uBLUP) <- u.names
}
rownames(G) <- u.names
colnames(G) <- u.names
rownames(Cu) <- u.names
colnames(Cu) <- u.names
}
if (MM > 1) {
colnames(xBLUP) <- id.names
colnames(condRes) <- id.names
# colnames(condResVar) <- id.names # NB: Non-commented if condResVar is included
}
# --------------------
# Return results
# --------------------
# if ((q0>0) && (p0>0)) warning("uVar is not corrected for estimation of fixed effect.\n To be updated later.")
return(list(logLik=res$logLik,ANOVA=anova.type,nlSearch=nlSearch,counts=counts,
boxcoxHat=boxcox,Ghat=G,lambdaHat=lambda,sigma2hat=res$sigma2hat,
betaHat=betaHat,uBLUP=uBLUP,xBLUP=xBLUP,
condRes=condRes,
# condResVar=condResVar,
# invAdiag=invAdiag,acf=acfVec,
betaVar=Cbeta
# ,uVar=Cu
))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.