## download nlme package source, unpack it; look in R/corStruct.R library(nlme) ?corClasses ## says corMatrix(), coef() are required ## download ape package, unpack it: look in R/PGLS.R, at corBrownian, corMartins etc. ## not quite sure how to manage constrained/unconstrained parameters: where are ## they used vs not used? ## copy safe_phi in corStruct.c ptrans_fromconstr <- function(x) { log((1+x)/(1-x)) } ptrans_toconstr <- function(x) { ex <- exp(-sign(x)*x) sign(x)*(1-ex)/(1+ex) } test <- seq(-0.9,0.9,by=0.1) stopifnot(all.equal(ptrans_toconstr(ptrans_fromconstr(test)),test)) corAnte <- function(value=0, form = ~1, nmax=20, fixed=FALSE) { ## possibly check values: here -1=1)) { stop("Parameters in Ante structure must be between -1 and 1") } ## transform values to unconstrained scale value <- ptrans_fromconstr(value) attr(value, "formula") <- form attr(value, "fixed") <- fixed class(value) <- c("corAnte", "corStruct") value } ## could get away without Initialize for simpler models, e.g. the ## corAR1 constructor starts with default value=0 and Initialize is ## unnecessary, but here we need it to determine the appropriate ## length of the parameter vector Initialize.corAnte <- function(object, data, ...) { form <- formula(object) ## obtaining the groups information, if any if (!is.null(getGroupsFormula(form))) { attr(object, "groups") <- getGroups(object, form, data = data) attr(object, "Dim") <- Dims <- Dim(object, attr(object, "groups")) } else { # no groups attr(object, "Dim") <- Dims <- Dim(object, as.factor(rep(1, nrow(data)))) } ## obtaining the covariate(s) attr(object, "covariate") <- getCovariate(object, data = data) ## copied from Initialize.corSymm; assignment trashes attributes (including class), so save and restore oldAttr <- attributes(object) object <- numeric(Dims$maxLen-1) attributes(object) <- oldAttr object } corMatrix.corAnte <- function(object, covariate = getCovariate(object), corr = TRUE, debug=TRUE, ...) { ## sort out dimensions -- copied from corAR1 if (data.class(covariate) == "list") { if (is.null(names(covariate))) { names(covariate) <- 1:length(covariate) } corD <- Dim(object, rep(names(covariate), unlist(lapply(covariate, length)))) cm <- lapply(corD$len,corm_ante,rho=as.vector(object)) } else { corD <- Dim(object, rep(1, length(covariate))) ## corMatrix.corAR1 / AR1_mat in corStruct.c ## appear to use constrained parameters ## ptrans_toconstr(as.vector(object))) cm <- corm_ante(corD$len,rho=ptrans_toconstr(as.vector(object))) if (!corr) { cm <- t(solve(chol(cm))) ## ugh! } } if (debug) cat("in corMatrix.corAnte: (unconstr)",object[1:3],"(constr)",ptrans_toconstr(object[1:3]),"\n") return(cm) } ## wrapper coef.corAnte <- function(object, unconstrained=TRUE, ...) { if (unconstrained) { if (attr(object, "fixed")) { return(numeric(0)) } else { return(as.vector(object)) } } aux <- ptrans_toconstr(as.vector(object)) names(aux) <- paste("Phi",seq(length(aux)),sep=".") aux } ## inner function (actual computation of correlation matrix) corm_ante <- function(n,rho) { m <- matrix(1,n,n) ## C_{ij} = \prod_{k=i}^{(j-1)} rho_k) ## fill upper triangle ## there is almost certainly a better way to do this via cumprod() for (i in seq(n-1)) for (j in 2:n) { if (j>i) { m[i,j] <- prod(rho[i:(j-1)]) } } m[lower.tri(m)] <- t(m)[lower.tri(m)] ## symmetrize m } ## test corm_ante(3,rho=seq(0.1,0.2,by=0.1)) ## should be: ## C(1,2) = C(2,1) = rho_1 = 0.1 ## C(1,3) = C(3,1) = rho_1*rho_2 = 0.02 ## C(2,3) = C(3,2) = rho_2 = 0.02 rho <- seq(0.9,0.1,by=-0.1) ## rho <- rep(0,9) m <- corm_ante(length(rho)+1,rho) N <- 20 ## generate data library(MASS) library(reshape) set.seed(1001) z0 <- mvrnorm(N,mu=rep(0,nrow(m)),Sigma=m,empirical=TRUE) z <- melt(z0) round(cor(z0),2) round(m,2) ## compute raw correlations ?? dim(cast(z,...~X1+X2)) library(nlme) ## debug(nlme:::corMatrix.corAR1) ## debug(nlme:::corFactor.corAR1) ## debug(nlme:::Initialize.corAR1) ## debug(nlme:::Initialize.corARMA) ## debug(nlme:::Initialize.corStruct) ## debug(corMatrix.corAnte) ## debug(Initialize.corAnte) ## debug(corAnte) ## debug(nlme:::Initialize.glsStruct) ## gls(value~1,data=z,correlation=corAR1(form=~X1|X2)) g1 <- gls(value~1,data=z,correlation=corAnte(form=~X2|X1)) g1 coef(g1$modelStruct$corStruct,unconstrained=TRUE) coef(g1$modelStruct$corStruct,unconstrained=FALSE) ## OK, I still have something backwards here