| @@ -1,123 +1,191 @@ | ||
| # | ||
| # growth : A Library of Normal Distribution Growth Curve Models | ||
| # Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey | ||
| # | ||
| # This program is free software; you can redistribute it and/or modify | ||
| # it under the terms of the GNU General Public Licence as published by | ||
| # the Free Software Foundation; either version 2 of the Licence, or | ||
| # (at your option) any later version. | ||
| # | ||
| # This program is distributed in the hope that it will be useful, | ||
| # but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| # GNU General Public Licence for more details. | ||
| # | ||
| # You should have received a copy of the GNU General Public Licence | ||
| # along with this program; if not, write to the Free Software | ||
| # Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. | ||
| # | ||
| # SYNOPSIS | ||
| # | ||
| # pergram(y) | ||
| # corgram(y, wt=1, add=FALSE, lty=1, xlim=NULL, ylim=NULL, xlab=NULL, | ||
| # ylab=NULL, main=NULL) | ||
| # | ||
| # DESCRIPTION | ||
| # | ||
| # Functions to compute and plot correlograms and periodograms | ||
|
|
||
| ### periodogram functions | ||
| ### | ||
| pergram <- function(y){ | ||
| ll <- length(y) | ||
| len <- trunc(ll/2) | ||
| len1 <- if(ll==2*len)len+1 else len+2 | ||
| fc <- fft(y) | ||
| z <- cbind(2*pi*(1:len)/ll, | ||
| (Re(fc[ll:len1])^2+Im(fc[ll:len1])^2)*ll/4/pi) | ||
| class(z) <- "pergram" | ||
| z} | ||
|
|
||
| #pergram <- function(y){ | ||
| # ll <- length(y) | ||
| # len <- trunc(ll/2) | ||
| # fc <- matrix(0,ncol=2,nrow=ll) | ||
| # for(i in 1:ll){ | ||
| # a <- 2*pi*i*(0:(ll-1))/ll | ||
| # fc[i,1] <- sum(y*sin(a)) | ||
| # fc[i,2] <- sum(y*cos(a))} | ||
| # invisible(cbind(2*pi*(1:len)/ll,(fc[,1]^2+fc[,2]^2)*ll/4/pi))} | ||
|
|
||
| plot.pergram <- function(y, add=FALSE, lty=1, xlab="Frequency", | ||
| ylab="Periodogram", main="Periodogram", ylim=c(0,max(po[,2]))){ | ||
| if(inherits(y,"pergram"))po <- y | ||
| else po <- pergram(y) | ||
| if(add)lines(po[,1],po[,2],lty=lty) | ||
| else plot(po[,1],po[,2],xlim=c(0,3.6),xlab=xlab,ylab=ylab,type="l", | ||
| lty=lty,ylim=ylim) | ||
| invisible(po)} | ||
|
|
||
| plot.cum <- function(z, ...) | ||
| UseMethod("plot.cum") | ||
|
|
||
| plot.cum.pergram <- function(y, xlab="Frequency", ylab="Periodogram", | ||
| main="Cumulative periodogram", | ||
| ylim=c(0,max(cpo+1.358/(a+0.12+0.11/a)))){ | ||
| if(inherits(y,"pergram")){ | ||
| len <- 2*dim(y)[1] | ||
| po <- y} | ||
| else { | ||
| len <- length(y) | ||
| po <- pergram(y)} | ||
| cpo <- cumsum(po[,2]) | ||
| cpo <- cpo/max(cpo) | ||
| a <- sqrt(len) | ||
| pa <- 2*pi*(1:length(cpo))/len | ||
| x <- (1:dim(po)[1])/dim(po)[1] | ||
| plot(po[,1],cpo,xlim=c(0,3.6),xlab=xlab,ylab=ylab,type="l", | ||
| ylim=ylim,main=main) | ||
| lines(po[,1],x+1.358/(a+0.12+0.11/a),lty=3) | ||
| lines(po[,1],x-1.358/(a+0.12+0.11/a),lty=3)} | ||
| ### | ||
| ### correlogram function | ||
| ### | ||
| corgram <- function(y, wt=1, maxlag=NULL, partial=FALSE, add=FALSE, lty=1, | ||
| xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL, main=NULL, ...){ | ||
| if(any(wt!=0&&wt!=1))stop("weigts must be 0 or 1") | ||
| len <- length(y) | ||
| if(length(wt)==1)wt <- rep(1,len) | ||
| if(any(is.na(y)))wt[is.na(y)] <- 0 | ||
| wt <- wt>0 | ||
| num <- sum(wt) | ||
| limit <- 2/sqrt(num) | ||
| if(is.null(maxlag))maxlag <- trunc(num/4) | ||
| ave <- sum(y*wt,na.rm=TRUE)/num | ||
| var <- sum((y-ave)^2*wt,na.rm=TRUE)/(num-1) | ||
| co <- rep(1,maxlag+1) | ||
| for(lag in 1:maxlag) | ||
| co[lag+1] <- sum((y[1:(len-lag)]-ave)*(y[(lag+1):len]-ave)* | ||
| wt[1:(len-lag)]*wt[(lag+1):len],na.rm=TRUE)/ | ||
| var/(sum(wt[1:(len-lag)]*wt[(lag+1):len],na.rm=TRUE)-1) | ||
| if(partial){ | ||
| a <- tmp <- rep(0,maxlag) | ||
| a[1] <- tmp[1] <- co[2] | ||
| for(lag in 2:maxlag){ | ||
| a[lag] <- tmp[lag] <- (co[lag+1]-sum(a[1:(lag-1)]* | ||
| co[lag:2]))/(1-sum(a[1:(lag-1)]*co[2:lag])) | ||
| a[1:(lag-1)] <- a[1:(lag-1)]-a[lag]*a[(lag-1):1]} | ||
| co <- tmp} | ||
| start <- if(partial)1 else 0 | ||
| if(is.null(xlim))xlim <- c(0,maxlag) | ||
| if(is.null(ylim))ylim <- c(if(min(co)<0)min(c(min(co),-limit))else 0,1) | ||
| if(is.null(xlab))xlab <- "Lag" | ||
| if(is.null(ylab))ylab <- expression(rho) | ||
| if(is.null(main))main <- if(partial)"Partial Correlogram" | ||
| else"Correlogram" | ||
| if(add)lines(start:maxlag,co,lty=lty) | ||
| else { | ||
| plot(start:maxlag,co,type="l",lty=lty,xlim=xlim,ylim=ylim, | ||
| xlab=xlab,ylab=ylab,main=main,...) | ||
| abline(h=limit,lty=3) | ||
| if(min(co)<0)abline(h=-limit,lty=3) | ||
| abline(h=0)} | ||
| invisible(cbind(start:maxlag,co))} | ||
| # | ||
| # growth : A Library of Normal Distribution Growth Curve Models | ||
| # Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey | ||
| # | ||
| # This program is free software; you can redistribute it and/or modify | ||
| # it under the terms of the GNU General Public Licence as published by | ||
| # the Free Software Foundation; either version 2 of the Licence, or | ||
| # (at your option) any later version. | ||
| # | ||
| # This program is distributed in the hope that it will be useful, | ||
| # but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| # GNU General Public Licence for more details. | ||
| # | ||
| # You should have received a copy of the GNU General Public Licence | ||
| # along with this program; if not, write to the Free Software | ||
| # Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. | ||
| # | ||
| # SYNOPSIS | ||
| # | ||
| # pergram(y) | ||
| # corgram(y, wt=1, add=FALSE, lty=1, xlim=NULL, ylim=NULL, xlab=NULL, | ||
| # ylab=NULL, main=NULL) | ||
| # | ||
| # DESCRIPTION | ||
| # | ||
| # Functions to compute and plot correlograms and periodograms | ||
|
|
||
| ### periodogram functions | ||
| ### | ||
|
|
||
|
|
||
| #' Calculate and Plot a Periodogram | ||
| #' | ||
| #' \code{pergram} calculates the values of a periodogram, \code{plot.pergram} | ||
| #' plots it, and \code{plot.cum.pergram} plots the corresponding cumulative | ||
| #' periodogram. | ||
| #' | ||
| #' | ||
| #' @aliases pergram plot.pergram plot_cum plot_cum.pergram | ||
| #' @param y A time series vector. | ||
| #' @param add If TRUE, adds a new periodogram to an existing plot. | ||
| #' @param x Plotting parameters | ||
| #' @param lty Plotting parameters | ||
| #' @param xlab Plotting parameters | ||
| #' @param ylab Plotting parameters | ||
| #' @param main Plotting parameters | ||
| #' @param ylim Plotting parameters | ||
| #' @param ... Plotting parameters | ||
| #' @return \code{pergram} prints and returns a two-column matrix of class, | ||
| #' \code{pergram}, containing the periodogram. | ||
| #' @author J.K. Lindsey | ||
| #' @keywords hplot | ||
| #' @examples | ||
| #' | ||
| #' y <- rnorm(100) | ||
| #' print(z <- pergram(y)) | ||
| #' plot(z) | ||
| #' plot_cum(z) | ||
| #' | ||
| #' @export | ||
| pergram <- function(y){ | ||
| ll <- length(y) | ||
| len <- trunc(ll/2) | ||
| len1 <- if(ll==2*len)len+1 else len+2 | ||
| fc <- fft(y) | ||
| z <- cbind(2*pi*(1:len)/ll, | ||
| (Re(fc[ll:len1])^2+Im(fc[ll:len1])^2)*ll/4/pi) | ||
| class(z) <- "pergram" | ||
| z} | ||
|
|
||
| #pergram <- function(y){ | ||
| # ll <- length(y) | ||
| # len <- trunc(ll/2) | ||
| # fc <- matrix(0,ncol=2,nrow=ll) | ||
| # for(i in 1:ll){ | ||
| # a <- 2*pi*i*(0:(ll-1))/ll | ||
| # fc[i,1] <- sum(y*sin(a)) | ||
| # fc[i,2] <- sum(y*cos(a))} | ||
| # invisible(cbind(2*pi*(1:len)/ll,(fc[,1]^2+fc[,2]^2)*ll/4/pi))} | ||
| #' @describeIn pergram Plot method | ||
| #' @export | ||
| plot.pergram <- function(x, add=FALSE, lty=1, xlab="Frequency", | ||
| ylab="Periodogram", main="Periodogram", ylim=c(0,max(po[,2])), ...){ | ||
| y <- x | ||
| if(inherits(y,"pergram"))po <- y | ||
| else po <- pergram(y) | ||
| if(add)lines(po[,1],po[,2],lty=lty) | ||
| else plot(po[,1],po[,2],xlim=c(0,3.6),xlab=xlab,ylab=ylab,type="l", | ||
| lty=lty,ylim=ylim) | ||
| invisible(po)} | ||
|
|
||
| #' @export | ||
| plot_cum <- function(x, ...) UseMethod("plot_cum") | ||
|
|
||
| #' @describeIn pergram Plot_cum method | ||
| #' @export | ||
| plot_cum.pergram <- function(x, xlab="Frequency", ylab="Periodogram", | ||
| main="Cumulative periodogram", | ||
| ylim=c(0,max(cpo+1.358/(a+0.12+0.11/a))), ...){ | ||
| y <- x | ||
| if(inherits(y,"pergram")){ | ||
| len <- 2*dim(y)[1] | ||
| po <- y} | ||
| else { | ||
| len <- length(y) | ||
| po <- pergram(y)} | ||
| cpo <- cumsum(po[,2]) | ||
| cpo <- cpo/max(cpo) | ||
| a <- sqrt(len) | ||
| pa <- 2*pi*(1:length(cpo))/len | ||
| x <- (1:dim(po)[1])/dim(po)[1] | ||
| plot(po[,1],cpo,xlim=c(0,3.6),xlab=xlab,ylab=ylab,type="l", | ||
| ylim=ylim,main=main) | ||
| lines(po[,1],x+1.358/(a+0.12+0.11/a),lty=3) | ||
| lines(po[,1],x-1.358/(a+0.12+0.11/a),lty=3)} | ||
| ### | ||
| ### correlogram function | ||
| ### | ||
|
|
||
|
|
||
| #' Calculate and Plot a Correlogram | ||
| #' | ||
| #' \code{corgram} calculates the values of a correlogram (autocorrelation | ||
| #' function or ACF) and plots it. | ||
| #' | ||
| #' | ||
| #' @param y A time series vector. | ||
| #' @param maxlag Maximum number of lags for which the correlation is to be | ||
| #' calculated. | ||
| #' @param partial If TRUE, the partial autocorrelation function (PACF) is | ||
| #' plotted. | ||
| #' @param wt Indicator vector with zeros for values to be ignored. | ||
| #' @param add If TRUE, adds a new correlogram to an existing plot. | ||
| #' @param xlim Plotting parameters | ||
| #' @param lty Plotting parameters | ||
| #' @param xlab Plotting parameters | ||
| #' @param ylab Plotting parameters | ||
| #' @param main Plotting parameters | ||
| #' @param ylim Plotting parameters | ||
| #' @param ... Plotting parameters | ||
| #' @return \code{corgram} returns a two-column matrix containing the (partial) | ||
| #' correlogram coordinates. | ||
| #' @author J.K. Lindsey | ||
| #' @keywords hplot | ||
| #' @examples | ||
| #' | ||
| #' y <- rnorm(100) | ||
| #' corgram(y) | ||
| #' corgram(y, partial=TRUE) | ||
| #' @export | ||
| corgram <- function(y, wt=1, maxlag=NULL, partial=FALSE, add=FALSE, lty=1, | ||
| xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL, main=NULL, ...){ | ||
| if(any(wt!=0&&wt!=1))stop("weights must be 0 or 1") | ||
| len <- length(y) | ||
| if(length(wt)==1)wt <- rep(1,len) | ||
| if(any(is.na(y)))wt[is.na(y)] <- 0 | ||
| wt <- wt>0 | ||
| num <- sum(wt) | ||
| limit <- 2/sqrt(num) | ||
| if(is.null(maxlag))maxlag <- trunc(num/4) | ||
| ave <- sum(y*wt,na.rm=TRUE)/num | ||
| var <- sum((y-ave)^2*wt,na.rm=TRUE)/(num-1) | ||
| co <- rep(1,maxlag+1) | ||
| for(lag in 1:maxlag) | ||
| co[lag+1] <- sum((y[1:(len-lag)]-ave)*(y[(lag+1):len]-ave)* | ||
| wt[1:(len-lag)]*wt[(lag+1):len],na.rm=TRUE)/ | ||
| var/(sum(wt[1:(len-lag)]*wt[(lag+1):len],na.rm=TRUE)-1) | ||
| if(partial){ | ||
| a <- tmp <- rep(0,maxlag) | ||
| a[1] <- tmp[1] <- co[2] | ||
| for(lag in 2:maxlag){ | ||
| a[lag] <- tmp[lag] <- (co[lag+1]-sum(a[1:(lag-1)]* | ||
| co[lag:2]))/(1-sum(a[1:(lag-1)]*co[2:lag])) | ||
| a[1:(lag-1)] <- a[1:(lag-1)]-a[lag]*a[(lag-1):1]} | ||
| co <- tmp} | ||
| start <- if(partial)1 else 0 | ||
| if(is.null(xlim))xlim <- c(0,maxlag) | ||
| if(is.null(ylim))ylim <- c(if(min(co)<0)min(c(min(co),-limit))else 0,1) | ||
| if(is.null(xlab))xlab <- "Lag" | ||
| if(is.null(ylab))ylab <- expression(rho) | ||
| if(is.null(main))main <- if(partial)"Partial Correlogram" | ||
| else"Correlogram" | ||
| if(add)lines(start:maxlag,co,lty=lty) | ||
| else { | ||
| plot(start:maxlag,co,type="l",lty=lty,xlim=xlim,ylim=ylim, | ||
| xlab=xlab,ylab=ylab,main=main,...) | ||
| abline(h=limit,lty=3) | ||
| if(min(co)<0)abline(h=-limit,lty=3) | ||
| abline(h=0)} | ||
| invisible(cbind(start:maxlag,co))} |
| @@ -1,95 +1,135 @@ | ||
| # | ||
| # growth : A Library of Normal Distribution Growth Curve Models | ||
| # Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey | ||
| # | ||
| # This program is free software; you can redistribute it and/or modify | ||
| # it under the terms of the GNU General Public Licence as published by | ||
| # the Free Software Foundation; either version 2 of the Licence, or | ||
| # (at your option) any later version. | ||
| # | ||
| # This program is distributed in the hope that it will be useful, | ||
| # but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| # GNU General Public Licence for more details. | ||
| # | ||
| # You should have received a copy of the GNU General Public Licence | ||
| # along with this program; if not, write to the Free Software | ||
| # Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. | ||
| # | ||
| # SYNOPSIS | ||
| # | ||
| # rmaov(response, tvcov=NULL, ccov=NULL, analysis=TRUE) | ||
| # | ||
| # DESCRIPTION | ||
| # | ||
| # Function to calculate the ANOVA table for a balanced split-plot design | ||
| # Copyright (C) 1999 Ralf Goertz (ralf.goertz@uni-jena.de) | ||
|
|
||
| rmaov <- function(response, tvcov=NULL, ccov=NULL, analysis=TRUE){ | ||
| if(!is.matrix(response)||!is.numeric(response))stop("response must be a matrix") | ||
| if(is.null(tvcov)){ | ||
| nfac <- 1 | ||
| nrep <- nlev <- dim(response)[2] | ||
| rfacname <- "repfac"} | ||
| else if(is.vector(tvcov,mode="numeric")||is.factor(tvcov)){ | ||
| nfac <- 1 | ||
| nrep <- nlev <- length(unique(tvcov)) | ||
| rfacname <- paste(deparse(substitute(tvcov))) | ||
| tvcov <- as.data.frame(tvcov)} | ||
| else if(is.data.frame(tvcov)||is.matrix(tvcov)){ | ||
| nfac <- dim(tvcov)[2] | ||
| if(!is.null(colnames(tvcov)))rfacname <- colnames(tvcov) | ||
| else rfacname <- paste("repfac",1:dim(tvcov)[2],sep="") | ||
| nlev <- NULL | ||
| nrep <- 1 | ||
| for(i in 1:dim(tvcov)[2]){ | ||
| nlev <- c(nlev,length(unique(tvcov[,i]))) | ||
| nrep <- nrep*nlev[i]}} | ||
| else stop("If specified, tvcov must be a vector, matrix, or dataframe containing the repeated factors") | ||
| if(is.vector(ccov,mode="numeric")||is.factor(ccov)){ | ||
| tmp <- paste(deparse(substitute(ccov))) | ||
| ccov <- list(ccov) | ||
| names(ccov) <- tmp} | ||
| else if(is.matrix(ccov))ccov <- data.frame(ccov) | ||
| else if(is.list(ccov)&&!is.data.frame(ccov))for(i in 1:length(ccov)) | ||
| if(is.matrix(ccov[[i]]))ccov[[i]] <- as.vector(t(ccov[[i]])) | ||
| else if(!is.null(ccov))print("ccov must be a vector, matrix, dataframe, or list of time-constant covariates") | ||
| ncases <- dim(response)[2]*dim(response)[1] | ||
| if(nrep!=dim(response)[2]) | ||
| stop("Number of columns in response must match product of levels") | ||
| block <- gl(dim(response)[1],nrep) | ||
| nc <- n <- dim(response)[2] | ||
| response <- as.vector(t(response)) | ||
| for(i in 1:nfac){ | ||
| n <- n/nlev[i] | ||
| com <- paste("fm",as.character(i),"<-gl(nlev[i],n,ncases",sep="") | ||
| if(!is.null(levels(tvcov[,i]))) | ||
| com <- paste(com,",labels=levels(tvcov[,i])",sep="") | ||
| com <- paste(com,")\n",sep="") | ||
| eval(parse(text=com))} | ||
| if(!is.null(ccov)){ | ||
| tccov <- ccov | ||
| if(length(tccov[[1]])!=ncases)for(i in 1:length(tccov)){ | ||
| if(!is.factor(tccov[[i]])) | ||
| tccov[[i]] <- as.factor(tccov[[i]]) | ||
| ccov[[i]] <- gl(nlevels(tccov[[i]]),1,ncases,labels=levels(tccov[[i]])) | ||
| ccov[[i]][(rep(1:length(tccov[[i]]),rep(nc,length(tccov[[i]])))-1)*nc+rep(1:nc,length(tccov[[i]]))] <- tccov[[i]][rep(1:length(tccov[[i]]),rep(nc,length(tccov[[i]])))]}} | ||
| com <- paste("b<-data.frame(response=response,block=block") | ||
| for(i in 1:nfac)com <- paste(com,",",rfacname[i],"=fm",as.character(i),sep="") | ||
| if(!is.null(ccov)){ | ||
| if(is.null(names(ccov))) | ||
| names(ccov) <- paste("ccov",1:length(ccov),sep="") | ||
| for(i in 1:length(ccov)) | ||
| com <- paste(com,",",names(ccov)[i],"=ccov$",names(ccov)[i],sep="")} | ||
| com <- paste(com,")\n",sep="") | ||
| eval(parse(text=com)) | ||
| if(!analysis)return(b) | ||
| efacterm <- facterm <- paste(rfacname,collapse="*") | ||
| if(!is.null(ccov)) | ||
| for(i in 1:length(ccov)) | ||
| facterm <- paste(facterm,names(ccov)[i],sep="*") | ||
| com <- paste("res<-aov(response~",facterm,"+Error(block/(",efacterm,")),data=b)\n",sep="") | ||
| eval(parse(text=com)) | ||
| res} | ||
|
|
||
|
|
||
| # | ||
| # growth : A Library of Normal Distribution Growth Curve Models | ||
| # Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey | ||
| # | ||
| # This program is free software; you can redistribute it and/or modify | ||
| # it under the terms of the GNU General Public Licence as published by | ||
| # the Free Software Foundation; either version 2 of the Licence, or | ||
| # (at your option) any later version. | ||
| # | ||
| # This program is distributed in the hope that it will be useful, | ||
| # but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| # GNU General Public Licence for more details. | ||
| # | ||
| # You should have received a copy of the GNU General Public Licence | ||
| # along with this program; if not, write to the Free Software | ||
| # Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. | ||
| # | ||
| # SYNOPSIS | ||
| # | ||
| # rmaov(response, tvcov=NULL, ccov=NULL, analysis=TRUE) | ||
| # | ||
| # DESCRIPTION | ||
| # | ||
| # Function to calculate the ANOVA table for a balanced split-plot design | ||
| # Copyright (C) 1999 Ralf Goertz (ralf.goertz@uni-jena.de) | ||
|
|
||
|
|
||
|
|
||
| #' Split-plot ANOVA Model | ||
| #' | ||
| #' \code{rmaov} performs the classical balanced split-plot ANOVA, with | ||
| #' \code{summary} providing the table. This is the so-called repeated measures | ||
| #' ANOVA. | ||
| #' | ||
| #' For unbalanced data, \code{\link[growth]{elliptic}} will perform the | ||
| #' analysis for one or two levels of nesting. | ||
| #' | ||
| #' | ||
| #' @param response A matrix or dataframe of response values with units as rows | ||
| #' and repeated measures as columns. | ||
| #' @param tvcov A numeric vector or factor variable defining the clusters. If | ||
| #' there are several levels of nesting, a matrix or dataframe with columns of | ||
| #' such variables defining the nested clusters starting with the highest level | ||
| #' (that is, from slowest to fastest varying). If not provided, each response | ||
| #' value of a unit is assumed to belong to a different cluster (that is, one | ||
| #' factor with \code{ncol(response)} levels is assumed). | ||
| #' @param ccov A vector or factor variable for one inter-subject covariate or a | ||
| #' matrix, dataframe, or list of several such variables. | ||
| #' @param analysis If FALSE, the design matrix is set up, but the analysis is | ||
| #' not performed. | ||
| #' @return The fitted model is returned. | ||
| #' @author Ralf Goertz (ralf.goertz@@uni-jena.de) | ||
| #' @seealso \code{\link[growth]{carma}}, \code{\link[growth]{elliptic}}, | ||
| #' \code{\link{lm}}, \code{\link[growth]{potthoff}}. | ||
| #' @keywords models | ||
| #' @examples | ||
| #' | ||
| #' # vision data for 7 individuals, with response a 7x8 matrix | ||
| #' # two levels of nesting: 4 levels of power for each eye | ||
| #' y <- matrix(rnorm(56),ncol=8) | ||
| #' tvc <- data.frame(eye=c(rep(1,4),rep(2,4)),power=c(1:4,1:4)) | ||
| #' summary(rmaov(y, tvc)) | ||
| #' | ||
| #' @export rmaov | ||
| rmaov <- function(response, tvcov=NULL, ccov=NULL, analysis=TRUE){ | ||
| if(!is.matrix(response)||!is.numeric(response))stop("response must be a matrix") | ||
| if(is.null(tvcov)){ | ||
| nfac <- 1 | ||
| nrep <- nlev <- dim(response)[2] | ||
| rfacname <- "repfac"} | ||
| else if(is.vector(tvcov,mode="numeric")||is.factor(tvcov)){ | ||
| nfac <- 1 | ||
| nrep <- nlev <- length(unique(tvcov)) | ||
| rfacname <- paste(deparse(substitute(tvcov))) | ||
| tvcov <- as.data.frame(tvcov)} | ||
| else if(is.data.frame(tvcov)||is.matrix(tvcov)){ | ||
| nfac <- dim(tvcov)[2] | ||
| if(!is.null(colnames(tvcov)))rfacname <- colnames(tvcov) | ||
| else rfacname <- paste("repfac",1:dim(tvcov)[2],sep="") | ||
| nlev <- NULL | ||
| nrep <- 1 | ||
| for(i in 1:dim(tvcov)[2]){ | ||
| nlev <- c(nlev,length(unique(tvcov[,i]))) | ||
| nrep <- nrep*nlev[i]}} | ||
| else stop("If specified, tvcov must be a vector, matrix, or dataframe containing the repeated factors") | ||
| if(is.vector(ccov,mode="numeric")||is.factor(ccov)){ | ||
| tmp <- paste(deparse(substitute(ccov))) | ||
| ccov <- list(ccov) | ||
| names(ccov) <- tmp} | ||
| else if(is.matrix(ccov))ccov <- data.frame(ccov) | ||
| else if(is.list(ccov)&&!is.data.frame(ccov))for(i in 1:length(ccov)) | ||
| if(is.matrix(ccov[[i]]))ccov[[i]] <- as.vector(t(ccov[[i]])) | ||
| else if(!is.null(ccov))print("ccov must be a vector, matrix, dataframe, or list of time-constant covariates") | ||
| ncases <- dim(response)[2]*dim(response)[1] | ||
| if(nrep!=dim(response)[2]) | ||
| stop("Number of columns in response must match product of levels") | ||
| block <- gl(dim(response)[1],nrep) | ||
| nc <- n <- dim(response)[2] | ||
| response <- as.vector(t(response)) | ||
| for(i in 1:nfac){ | ||
| n <- n/nlev[i] | ||
| com <- paste("fm",as.character(i),"<-gl(nlev[i],n,ncases",sep="") | ||
| if(!is.null(levels(tvcov[,i]))) | ||
| com <- paste(com,",labels=levels(tvcov[,i])",sep="") | ||
| com <- paste(com,")\n",sep="") | ||
| eval(parse(text=com))} | ||
| if(!is.null(ccov)){ | ||
| tccov <- ccov | ||
| if(length(tccov[[1]])!=ncases)for(i in 1:length(tccov)){ | ||
| if(!is.factor(tccov[[i]])) | ||
| tccov[[i]] <- as.factor(tccov[[i]]) | ||
| ccov[[i]] <- gl(nlevels(tccov[[i]]),1,ncases,labels=levels(tccov[[i]])) | ||
| ccov[[i]][(rep(1:length(tccov[[i]]),rep(nc,length(tccov[[i]])))-1)*nc+rep(1:nc,length(tccov[[i]]))] <- tccov[[i]][rep(1:length(tccov[[i]]),rep(nc,length(tccov[[i]])))]}} | ||
| b <- NULL | ||
| com <- paste("b<-data.frame(response=response,block=block") | ||
| for(i in 1:nfac)com <- paste(com,",",rfacname[i],"=fm",as.character(i),sep="") | ||
| if(!is.null(ccov)){ | ||
| if(is.null(names(ccov))) | ||
| names(ccov) <- paste("ccov",1:length(ccov),sep="") | ||
| for(i in 1:length(ccov)) | ||
| com <- paste(com,",",names(ccov)[i],"=ccov$",names(ccov)[i],sep="")} | ||
| com <- paste(com,")\n",sep="") | ||
| eval(parse(text=com)) | ||
| if(!analysis)return(b) | ||
| efacterm <- facterm <- paste(rfacname,collapse="*") | ||
| if(!is.null(ccov)) | ||
| for(i in 1:length(ccov)) | ||
| facterm <- paste(facterm,names(ccov)[i],sep="*") | ||
| res <- NULL | ||
| com <- paste("res<-aov(response~",facterm,"+Error(block/(",efacterm,")),data=b)\n",sep="") | ||
| eval(parse(text=com)) | ||
| res} | ||
|
|
||
|
|
| @@ -0,0 +1,10 @@ | ||
|
|
||
| <!-- README.md is generated from README.Rmd. Please edit README.Rmd --> | ||
| [](https://travis-ci.org/swihart/growth) [](https://cran.r-project.org/package=growth)  | ||
|
|
||
| `growth` R package | ||
| ================== | ||
|
|
||
| This package is intended to be the developmental version to the CRAN version of [Jim Lindsey's growth](http://www.commanster.eu/rcode.html). The .zip files listed on his homepage have been listed as version 1.0 since 2005. For the subsequent maintenance on this github and CRAN, we will start at 1.1. | ||
|
|
||
| To compare this version with the static v1.0 files on [Jim Lindsey's Homepage](http://www.commanster.eu/rcode.html), it may be useful to use [the compare page for this repo's two branches](https://github.com/swihart/growth/compare/jim-lindsey-homepage-version-1.0...master?diff=split&name=master). |
| @@ -0,0 +1,69 @@ | ||
| # growth R package | ||
| Bruce Swihart | ||
| Feb 2019 | ||
|
|
||
| ## Re-Submission 1 | ||
|
|
||
| * fixed an `Found no calls to: ‘R_registerRoutines’, ‘R_useDynamicSymbols’` NOTE. | ||
| * corrected Fatal error: length > 1 in coercion to logical (raised on R-devel check) | ||
|
|
||
| ## Test environments | ||
| * local OS X install: R version 3.5.2 (2018-12-20) -- "Eggshell Igloo" | ||
| * Ubuntu 14.04.5 LTS (on travis-ci): R version 3.5.2 (2017-01-27) | ||
| * win-builder: R Under development (unstable) (2019-01-31 r76038) | ||
|
|
||
| ## R CMD check results | ||
| There were no ERRORs or WARNINGs or NOTEs. | ||
|
|
||
|
|
||
| ## Downstream dependencies | ||
| There are currently no downstream dependencies for this package. | ||
|
|
||
|
|
||
|
|
||
| # growth R package | ||
| Bruce Swihart | ||
| December 2016 | ||
|
|
||
| ## Test environments | ||
| * local OS X install: R version 3.3.2 (2016-10-31) | ||
| * ubuntu 12.04 (on travis-ci): R version 3.3.1 (2016-06-21) | ||
| * win-builder: R Under development (unstable) (2016-12-05 r71733) | ||
|
|
||
| ## R CMD check results | ||
| There were no ERRORs or WARNINGs. | ||
|
|
||
| There was 1 NOTE: | ||
|
|
||
| * 1st submission (Package was archived on CRAN, see Miscellaneous) | ||
|
|
||
| ## Downstream dependencies | ||
| There are currently no downstream dependencies for this package. | ||
|
|
||
| ## Miscellaneous | ||
| As per the CRAN Repository Policy Version Revision 3577, | ||
|
|
||
| * Explain any change in the maintainer’s email address and if possible send confirmation from the previous address (by a separate email to CRAN@R-project.org) or explain why it is not possible. | ||
|
|
||
| This Package was archived on CRAN. I am resurrecting it. | ||
| I have Jim Lindsey's permission to be maintainer of his packages on CRAN. Currently he has them in .zip files on his homepage: http://www.commanster.eu/rcode.html. He does not have access to the original email address he used when this package was on CRAN, and thus cannot send a separate confirmation email. However, you can contact him at James Lindsey <jlindsey@gen.unimaas.nl>. | ||
|
|
||
| He sent this confirmation message on 2016-11-24: | ||
|
|
||
| ---------- Forwarded message ---------- | ||
| From: James Lindsey <jlindsey@gen.unimaas.nl> | ||
| Date: Thu, Nov 24, 2016 at 9:26 AM | ||
| Subject: transferring maintainer role to Bruce Swihart for 'rmutil', 'repeated', 'gnlm', 'growth', 'event', 'stable' | ||
| To: CRAN@r-project.org, ligges@statistik.tu-dortmund.de, Kurt.Hornik@wu.ac.at | ||
| Cc: bruce.swihart@gmail.com | ||
|
|
||
|
|
||
| Dear CRAN, | ||
|
|
||
| Bruce Swihart is taking on the role of maintainer for my R-packages which have been available on my homepage, http://www.commanster.eu/rcode.html. Some of these R-packages were on CRAN but have since been archived. | ||
| I know it is CRAN repository policy to email confirmation from the previous maintainer's email address (jlindsey@luc.ac.be), but alas I can't due to LUC changing its name. | ||
|
|
||
| Please accept this email as confirmation of the maintainer role changing for R packages 'rmutil', 'repeated', 'gnlm', 'growth', 'event', 'stable'. | ||
|
|
||
| Regards, | ||
| Jim Lindsey |
| @@ -0,0 +1,33 @@ | ||
| #include <R_ext/RS.h> | ||
| #include <stdlib.h> // for NULL | ||
| #include <R_ext/Rdynload.h> | ||
|
|
||
| /* FIXME: | ||
| Check these declarations against the C/Fortran source code. | ||
| */ | ||
|
|
||
| /* .Fortran calls */ | ||
| extern void F77_NAME(back)(void *, void *); | ||
| extern void F77_NAME(kalman)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); | ||
| extern void F77_NAME(plra)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); | ||
| extern void F77_NAME(resid)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); | ||
| extern void F77_NAME(resid2)(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); | ||
| extern void F77_NAME(roots)(void *, void *, void *); | ||
| extern void F77_NAME(ttvert)(void *, void *); | ||
|
|
||
| static const R_FortranMethodDef FortranEntries[] = { | ||
| {"back", (DL_FUNC) &F77_NAME(back), 2}, | ||
| {"kalman", (DL_FUNC) &F77_NAME(kalman), 25}, | ||
| {"plra", (DL_FUNC) &F77_NAME(plra), 43}, | ||
| {"resid", (DL_FUNC) &F77_NAME(resid), 24}, | ||
| {"resid2", (DL_FUNC) &F77_NAME(resid2), 16}, | ||
| {"roots", (DL_FUNC) &F77_NAME(roots), 3}, | ||
| {"ttvert", (DL_FUNC) &F77_NAME(ttvert), 2}, | ||
| {NULL, NULL, 0} | ||
| }; | ||
|
|
||
| void R_init_growth(DllInfo *dll) | ||
| { | ||
| R_registerRoutines(dll, NULL, NULL, FortranEntries, NULL); | ||
| R_useDynamicSymbols(dll, FALSE); | ||
| } |