diff --git a/R/act_plot.R b/R/act_plot.R index df6ab63..3512794 100644 --- a/R/act_plot.R +++ b/R/act_plot.R @@ -1,106 +1,106 @@ -#' trelliscope auxillary function -#' -#' Plot day observations. -#' -#' @param id individual ID. -#' @param id_Nday the ith day observation of the individual. -#' @param act_ori day observation. -#' @param act_ind individual mean observation. -#' @param act_all global mean observation. -#' @param act_max the upper limit of activity (y-axis). -#' @param band smoothing band paramter. -#' -#' @keywords internal -#' -#' @return individual to day observations. - -#### activity plot -act_plot <- function(id,id_Nday,act_ori,act_ind,act_all,act_max,band,ori=FALSE,lw=TRUE) { - timeOfDay <- function(num) { - hour <- floor(num/60) - min <- num - 60*hour - min <- paste(0,min,sep="") - min <- substr(min,nchar(min)-1,nchar(min)) - return(paste(hour,":",min,sep="")) - } - timeOfDay <- Vectorize(timeOfDay,vectorize.args="num") - xtime <- timeOfDay(1:1440) - - len <- length(act_ind) - xseq <- seq(24/len,24,by=24/len) - d <- band*len - lws <- function(x,ff) lowess(x[c((len-d+1):len,1:len,1:d)],f=ff,iter=1)$y[(d+1):(d+len)] - vec <- lws(act_ori,band) - vec_ind <- lws(act_ind,band) - vec_all <- lws(act_all,band) - #vec <- act_ori - title_name <- paste("ID = ",id," #Day = ",id_Nday,sep="") - - if(lw==TRUE & ori==FALSE) { - rbokeh::figure(title=title_name,legend_location = "top_left", - width=600, height=300, - #xlab="Time of the Day", ylab="Activity Count", - xlim=c(0,24), ylim=c(0, max(c(act_max,act_ind,act_all))+50), - tools=c("pan","wheel_zoom","box_zoom","resize", - "box_select","reset","save","help")) %>% - #ly_points(xseq, act_ori, color="black", size=1, hover= "Time @xtime", alpha=0.5) %>% - ly_lines(xseq,act_ori,color="black",alpha=0.2,legend=FALSE) %>% - ly_points(xseq,vec,color="white",size=2,alpha=0, - hover=data.frame(Time=xtime,Activity=floor(vec)),legend=FALSE) %>% - ly_lines(xseq,vec,color="black",width=2,legend="activity lowess") %>% - ly_lines(xseq,vec_ind,width=2,color="blue",legend="individual mean lowess") %>% - ly_lines(xseq,vec_all,color="red",width=2,legend="global mean lowess") %>% - x_axis(label="Time of the Day", - desired_num_ticks = 10, num_minor_ticks= 2) %>% - y_axis(label="Activity Count") %>% - theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", - axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% - theme_title(text_font_size="12pt",text_font="Helvetica") %>% - theme_legend(label_text_font_size = "8pt") - } else if (lw==TRUE & ori==TRUE) { - rbokeh::figure(title=title_name,legend_location = "top_left", - width=600, height=300, - #xlab="Time of the Day", ylab="Activity Count", - xlim=c(0,24), ylim=c(0, max(c(act_max,act_ind,act_all))+50), - tools=c("pan","wheel_zoom","box_zoom","resize", - "box_select","reset","save","help")) %>% - #ly_points(xseq, act_ori, color="black", size=1, hover= "Time @xtime", alpha=0.5) %>% - ly_lines(xseq,act_ori,color="black",alpha=0.2,legend=FALSE) %>% - ly_points(xseq,vec,color="white",size=2,alpha=0, - hover=data.frame(Time=xtime,Activity=floor(vec)),legend=FALSE) %>% - ly_lines(xseq,act_ind,width=2,color="blue",alpha=0.2,legend="individual mean") %>% - ly_lines(xseq,act_all,color="red",alpha=0.2,width=2,legend="global mean") %>% - ly_lines(xseq,vec,color="black",width=2,legend="activity lowess") %>% - ly_lines(xseq,vec_ind,width=2,color="blue",legend="individual mean lowess") %>% - ly_lines(xseq,vec_all,color="yellow",width=2,legend="global mean lowess") %>% - x_axis(label="Time of the Day", - desired_num_ticks = 10, num_minor_ticks= 2) %>% - y_axis(label="Activity Count") %>% - theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", - axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% - theme_title(text_font_size="12pt",text_font="Helvetica") %>% - theme_legend(label_text_font_size = "8pt") - } else if(lw==FALSE & ori==TRUE) { - rbokeh::figure(title=title_name,legend_location = "top_left", - width=600, height=300, - #xlab="Time of the Day", ylab="Activity Count", - xlim=c(0,24), ylim=c(0, max(c(act_max,act_ind,act_all))+50), - tools=c("pan","wheel_zoom","box_zoom","resize", - "box_select","reset","save","help")) %>% - #ly_points(xseq, act_ori,color ="black", size=1, hover= "Time @xtime", alpha=0.5) %>% - ly_lines(xseq,act_ori,color="black",alpha=0.2,legend=FALSE) %>% - ly_points(xseq,vec,color="white",size=2,alpha=0, - hover=data.frame(Time=xtime,Activity=floor(vec)),legend=FALSE) %>% - ly_lines(xseq,act_ind,width=2,color="blue",alpha=0.2,legend="individual mean") %>% - ly_lines(xseq,act_all,color="red",alpha=0.2,width=2,legend="global mean") %>% - x_axis(label="Time of the Day", - desired_num_ticks = 10, num_minor_ticks= 2) %>% - y_axis(label="Activity Count") %>% - theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", - axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% - theme_title(text_font_size="12pt",text_font="Helvetica") %>% - theme_legend(label_text_font_size = "8pt") - } else { - stop("Error: at least plot one of the original and lowess!") - } -} +#' trelliscope auxillary function +#' +#' Plot day observations. +#' +#' @param id individual ID. +#' @param id_Nday the ith day observation of the individual. +#' @param act_ori day observation. +#' @param act_ind individual mean observation. +#' @param act_all global mean observation. +#' @param act_max the upper limit of activity (y-axis). +#' @param band smoothing band paramter. +#' +#' @keywords internal +#' +#' @return individual to day observations. + +#### activity plot +act_plot <- function(id,id_Nday,act_ori,act_ind,act_all,act_max,band,ori=FALSE,lw=TRUE) { + timeOfDay <- function(num) { + hour <- floor(num/60) + min <- num - 60*hour + min <- paste(0,min,sep="") + min <- substr(min,nchar(min)-1,nchar(min)) + return(paste(hour,":",min,sep="")) + } + timeOfDay <- Vectorize(timeOfDay,vectorize.args="num") + xtime <- timeOfDay(1:1440) + + len <- length(act_ind) + xseq <- seq(24/len,24,by=24/len) + d <- band*len + lws <- function(x,ff) lowess(x[c((len-d+1):len,1:len,1:d)],f=ff,iter=1)$y[(d+1):(d+len)] + vec <- lws(act_ori,band) + vec_ind <- lws(act_ind,band) + vec_all <- lws(act_all,band) + #vec <- act_ori + title_name <- paste("ID = ",id," #Day = ",id_Nday,sep="") + + if(lw==TRUE & ori==FALSE) { + rbokeh::figure(title=title_name,legend_location = "top_left", + width=600, height=300, + #xlab="Time of the Day", ylab="Activity Count", + xlim=c(0,24), ylim=c(0, max(c(act_max,act_ind,act_all))+50), + tools=c("pan","wheel_zoom","box_zoom","resize", + "box_select","reset","save","help")) %>% + #ly_points(xseq, act_ori, color="black", size=1, hover= "Time @xtime", alpha=0.5) %>% + ly_lines(xseq,act_ori,color="black",alpha=0.2,legend=FALSE) %>% + ly_points(xseq,vec,color="white",size=2,alpha=0, + hover=data.frame(Time=xtime,Activity=floor(vec)),legend=FALSE) %>% + ly_lines(xseq,vec,color="black",width=2,legend="activity lowess") %>% + ly_lines(xseq,vec_ind,width=2,color="blue",legend="individual mean lowess") %>% + ly_lines(xseq,vec_all,color="red",width=2,legend="global mean lowess") %>% + x_axis(label="Time of the Day", + desired_num_ticks = 10, num_minor_ticks= 2) %>% + y_axis(label="Activity Count") %>% + theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", + axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% + theme_title(text_font_size="12pt",text_font="Helvetica") %>% + theme_legend(label_text_font_size = "8pt") + } else if (lw==TRUE & ori==TRUE) { + rbokeh::figure(title=title_name,legend_location = "top_left", + width=600, height=300, + #xlab="Time of the Day", ylab="Activity Count", + xlim=c(0,24), ylim=c(0, max(c(act_max,act_ind,act_all))+50), + tools=c("pan","wheel_zoom","box_zoom","resize", + "box_select","reset","save","help")) %>% + #ly_points(xseq, act_ori, color="black", size=1, hover= "Time @xtime", alpha=0.5) %>% + ly_lines(xseq,act_ori,color="black",alpha=0.2,legend=FALSE) %>% + ly_points(xseq,vec,color="white",size=2,alpha=0, + hover=data.frame(Time=xtime,Activity=floor(vec)),legend=FALSE) %>% + ly_lines(xseq,act_ind,width=2,color="blue",alpha=0.2,legend="individual mean") %>% + ly_lines(xseq,act_all,color="red",alpha=0.2,width=2,legend="global mean") %>% + ly_lines(xseq,vec,color="black",width=2,legend="activity lowess") %>% + ly_lines(xseq,vec_ind,width=2,color="blue",legend="individual mean lowess") %>% + ly_lines(xseq,vec_all,color="yellow",width=2,legend="global mean lowess") %>% + x_axis(label="Time of the Day", + desired_num_ticks = 10, num_minor_ticks= 2) %>% + y_axis(label="Activity Count") %>% + theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", + axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% + theme_title(text_font_size="12pt",text_font="Helvetica") %>% + theme_legend(label_text_font_size = "8pt") + } else if(lw==FALSE & ori==TRUE) { + rbokeh::figure(title=title_name,legend_location = "top_left", + width=600, height=300, + #xlab="Time of the Day", ylab="Activity Count", + xlim=c(0,24), ylim=c(0, max(c(act_max,act_ind,act_all))+50), + tools=c("pan","wheel_zoom","box_zoom","resize", + "box_select","reset","save","help")) %>% + #ly_points(xseq, act_ori,color ="black", size=1, hover= "Time @xtime", alpha=0.5) %>% + ly_lines(xseq,act_ori,color="black",alpha=0.2,legend=FALSE) %>% + ly_points(xseq,vec,color="white",size=2,alpha=0, + hover=data.frame(Time=xtime,Activity=floor(vec)),legend=FALSE) %>% + ly_lines(xseq,act_ind,width=2,color="blue",alpha=0.2,legend="individual mean") %>% + ly_lines(xseq,act_all,color="red",alpha=0.2,width=2,legend="global mean") %>% + x_axis(label="Time of the Day", + desired_num_ticks = 10, num_minor_ticks= 2) %>% + y_axis(label="Activity Count") %>% + theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", + axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% + theme_title(text_font_size="12pt",text_font="Helvetica") %>% + theme_legend(label_text_font_size = "8pt") + } else { + stop("Error: at least plot one of the original and lowess!") + } +} diff --git a/R/bandSelect.R b/R/bandSelect.R index b784db6..704e9c2 100644 --- a/R/bandSelect.R +++ b/R/bandSelect.R @@ -1,168 +1,168 @@ -#' Penalized multi-band learning function -#' -#' In a group of individuals with physcial activity data, this function utilize Fast Fourier Transform (FFT) -#' and L-1/L-2 penalties to select significant harmonics/periodicities and describe the main activity pattern -#' (circadian rhythm) among the population. -#' -#' @param df the tbl_df data frame containing at least two variables: subject ID and activity. The function \code{form} can help prepare the data frame. -#' @param Nlength the length of observations necessary for each individual, note that it should be consistent among all -#' @param Nlambda \eqn{\lambda}'s take values from 0 to 2max(||X_k||^2), as 0 gives no penalty and the latter suppresses all \eqn{\theta}'s to 0. Therefore, we divide 2max(||X_k||^2) into Nlambda (default to be 100) \eqn{\lambda}'s to pick frequencies/harmonics/periodicities. -#' @param alpha the tuning parameter controlling the balance between L-1 and L-2 penalty. The default is 1, using complete Lasso/ L-1 penalty. -#' @param Ntop the number of frequencies/harmonics/periodicities picked for the population. The default is 5. -#' @param plot whether to plot: MSE against the number of nonzero \eqn{\theta}'s, and only the points at which the number of nonzero \eqn{\theta}'s changes (as \eqn{\lambda} changes) are be plotted. The default is TRUE. -#' @param cross whether to perform cross-validation. The default is FALSE. -#' @param Ncross the number of groups of data for cross-validation. If cross=TRUE, the data shall be divided into Ncross groups. -#' -#' @return if no cross-validation is conducted, return a list; if cross-validation, return a list of lists, with the last list consisting of all FFT results and cross-validation groups (showing the subject IDs leave-out /NOT used each time). -#' @return \item{topfreq}{vector of length \code{Ntop}: top frequencies selected.} -#' @return \item{mse}{vector of length \code{Nlambda}: mean squared error for each lambda (penalty). If no cross-validation, mse is calculated based on all available data; if cross-validation, mse is calculated based on the rest observations.} -#' @return \item{nonzero}{vector of length \code{Nlambda}: the number of nonzero \eqn{\theta}'s (frequencies) for each lambda (penalty).} -#' @return \item{deltazero}{vector of length \code{Nlambda}: the change in the number of nonzero \eqn{\theta}'s (frequencies) for each lambda (penalty).} -#' @return \item{lambda}{vector of length \code{Nlambda}: the value of lambda.} -#' @return \item{theta}{\code{Nfreq} by \code{Nlambda} matrix: estimated \eqn{\theta}'s (frequencies) at each lambda (penalty). Nfreq is the total number of frequencies given by FFT.} -#' @return \item{xscore}{\code{Nind} by \code{Nfreq} matrix: the original FFT scores for each individual. Nind is the number of individuals in the population, and Nfreq is the total number of frequencies given by FFT.} -#' @return \item{xprop}{\code{Nind} by \code{Nfreq} matrix: the original FFT results expressed as the proportion of variances explained by each frequency for each individual. Nind is the number of individuals in the population, and Nfreq is the total number of frequencies given by FFT.} -#' @return \item{freq}{vector of length \code{Nfreq}: list of frequencies in FFT results.} -#' -#' @examples -#' data(pa3) -#' re <- bandSelect(df=pa3,Nlength=1440*3,Nlambda=100,alpha=1,Ntop=5, -#' cross=FALSE,Ncross=NULL,plot=TRUE) -#' -#' @seealso \code{\link{form}} -#' -#' @references Li, X. , Kane, M. , Zhang, Y. , Sun, W. , Song, Y. , Dong, S. , Lin, Q. , Zhu, Q. , Jiang, F. & Zhao, H. (2019). Penalized Selection of Periodicities Characterizes the Consolidation of Sleep-Wake Circadian Rhythms During Early Childhood Development. Submitted. -#' -#' @export -#' - -bandSelect <- function(df,Nlength,Nlambda=100,alpha=1,Ntop=5,cross=FALSE,Ncross=NULL,plot=TRUE) { - ### preparation - # Note: requre harmonic_func.R - dflist <- lapply(split(df$activity,df$ID),unlist) - re <- list() - for(i in 1:length(dflist)) { - if(length(dflist[[i]])%' -#' @importFrom"dplyr" -#' mutate -#' -#' @param lis the list of activity data, with each element corresponding to the observation by one individual and the name of each element coresponding to the individual id. Specifically, each element is a \code{nob} by \code{nday} matrix, where each column is an observation by day. -#' @param maxday the maximal number of days per individual in the observation, used to check the data format. The default is 14. -#' @param id a vector of id names corresponding to the \code{lis} activity data. -#' -#' @keywords trelliscope -#' -#' @return The activity data frame with 3 columns: ID, IDday, and activity. - -#' @examples -#' data(lis3) -#' pa3 <- form(lis3) -#' -#' @seealso \code{\link{bandSelect}} -#' -#' @export - -form <- function(lis,maxday=14,id=NULL) { - checkformat <- do.call("c",lapply(lis,ncol)) - if(length(checkformat)!=length(lis)) { - stop("Contain empty matrix in the data list.") - } - if(max(checkformat)>maxday) { - stop(paste("Maxday ",maxday," reached: data list format may be wrong / consider change maxday.",sep="")) - } - - #### ID info - ID <- list() - if(!is.null(id)) { - if(length(id)!=length(lis)) { - stop("The length of the ID vector is not the same as the length of the data list.") - } - for (i in 1:length(id)) ID[[i]] <- rep(names(id[i]),ncol(lis[[i]])) - } else { - if(is.null(names(lis))) { - stop("Names of the data list are null: ID should be supplied.") - } - for (i in 1:length(lis)) ID[[i]] <- rep(names(lis[i]),ncol(lis[[i]])) - } - - #### ID_Nday info - ID_Nday <- NULL #required to avoid NOTE in CMD check - act <- data.frame(ID=as.numeric(unlist(ID)), - ID_Nday=unlist(lapply(ID,function(x) seq(1,length(x))))) - act <- dplyr::tbl_df(act) - act <- dplyr::group_by(act,ID,ID_Nday) - act <- tidyr::nest(act) - - #### activity - liscol <- do.call("cbind",lis) - liscol2 <- list() - for(i in 1:ncol(liscol)) liscol2[[i]] <- liscol[,i] - act <- dplyr::mutate(act,activity=liscol2);rm(liscol,liscol2) - act$data <- NULL - - return(act) -} - +#' Function to generate activity data frame for penalized multi-band learning +#' +#' This function generates the data frame necessary for further penalized multi-band learning. +#' +#' @importFrom"dplyr" +#' '%>%' +#' @importFrom"dplyr" +#' mutate +#' @importFrom"tibble" +#' add_column +#' +#' @param lis the list of activity data, with each element corresponding to the observation by one individual and the name of each element coresponding to the individual id. Specifically, each element is a \code{nob} by \code{nday} matrix, where each column is an observation by day. +#' @param maxday the maximal number of days per individual in the observation, used to check the data format. The default is 14. +#' @param id a vector of id names corresponding to the \code{lis} activity data. +#' +#' @keywords trelliscope +#' +#' @return The activity data frame with 3 columns: ID, IDday, and activity. + +#' @examples +#' data(lis3) +#' pa3 <- form(lis3) +#' +#' @seealso \code{\link{bandSelect}} +#' +#' @export + +form <- function(lis,maxday=14,id=NULL) { + checkformat <- do.call("c",lapply(lis,ncol)) + if(length(checkformat)!=length(lis)) { + stop("Contain empty matrix in the data list.") + } + if(max(checkformat)>maxday) { + stop(paste("Maxday ",maxday," reached: data list format may be wrong / consider change maxday.",sep="")) + } + + #### ID info + ID <- list() + if(!is.null(id)) { + if(length(id)!=length(lis)) { + stop("The length of the ID vector is not the same as the length of the data list.") + } + for (i in 1:length(id)) ID[[i]] <- rep(names(id[i]),ncol(lis[[i]])) + } else { + if(is.null(names(lis))) { + stop("Names of the data list are null: ID should be supplied.") + } + for (i in 1:length(lis)) ID[[i]] <- rep(names(lis[i]),ncol(lis[[i]])) + } + + #### ID_Nday info + ID_Nday <- NULL #required to avoid NOTE in CMD check + act <- data.frame(ID=as.numeric(unlist(ID)), + ID_Nday=unlist(lapply(ID,function(x) seq(1,length(x))))) + act <- dplyr::tbl_df(act) + act <- dplyr::group_by(act,ID,ID_Nday) + act <- tidyr::nest(act) + + #### activity + liscol <- do.call("cbind",lis) + liscol2 <- list() + for(i in 1:ncol(liscol)) liscol2[[i]] <- liscol[,i] + #act <- dplyr::mutate(act, activity = liscol2) note: due to the changed version of tidyr + act <- tibble::add_column(act,activity = liscol2) + rm(liscol,liscol2) + act$data <- NULL + + return(act) +} + diff --git a/R/gharmonic.R b/R/gharmonic.R index 0f53fdf..fac977e 100644 --- a/R/gharmonic.R +++ b/R/gharmonic.R @@ -1,102 +1,102 @@ -#' harmonic analysis test: g-value calculation -#' -#' This function calculates the g-value for the harmonic analysis test developed by R.A. Fisher (1929). -#' Harmonic analysis refers to Fast Fourier Transform (FFT) results. -#' Specifically, g is the proportion (squared modulus of one frequency divided by the sum of all squared moduli). -#' In order for g to be statistically significant in the harmonic analysis test, it needs to be at least g-value -#' at significance level \eqn{\alpha}. Please note that for the rth largest frequency, if any of the previous (r-1) -#' frequencies is not significant, then the rth largest frequency is also non-significant. -#' -#' @param n the total number of frequencies in FFT results. -#' @param r the modulus of the tested frequency is ranked as the rth largest among all frequencies. -#' @param p the FFT result of the tested frequency expressed as the squared modulus divided by the sum of the squared moduli by all frequencies (proportion: m_r^2/(m_1^2+...+m_n^2)). -#' @param tol the tolerance level during calculation. The default is 10^-7. -#' @param init the crude estimate for g-value if known. It is not called to calculate usual g-values. -#' @keywords harmonic test -#' -#' @return The g-value calculated by the harmonic test. - -#' @examples -#' gharmonic(n=100,r=1,p=0.05) -#' -#' @seealso \code{\link{pharmonic}} -#' -#' @references Fisher, R. A. (1929). Tests of significance in harmonic analysis. Proceedings of the Royal Society of London. Series A, 125(796), 54-59. -#' -#' @export - -gharmonic <- function(n,r,p,tol=10^-7,init=NULL) { - calc <- function(nn,rr,g) { - re <- 0 - k <- ifelse(floor(1/g)==1/g,1/g-1,floor(1/g)) - if(k0) { - temp <- lfactorial(nn)-lfactorial(rr-1)+(nn-1)*log(1-j*g)-log(j)-lfactorial(nn-j)-lfactorial(j-r) - if(exp(temp)==0) break - re <- re+(-1)^(j-rr)*exp(temp) - } else { - temp <- lfactorial(nn)-lfactorial(rr-1)+(nn-1)*log(j*g-1)-log(j)-lfactorial(nn-j)-lfactorial(j-r) - if(exp(temp)==0) break - re <- re+(-1)^(j-rr+nn-1)*exp(temp) - } - } - return(re) - } - calc <- Vectorize(calc,vectorize.args ="g") - #<70: 0.1-0.9 - #<983: 0.01-0.1 - if(is.null(init)) { - first <- seq(0.1,0.9,by=0.1) - firstre <- calc(nn=n,rr=r,g=first) - if(!any(firstre>p & firstre<1,na.rm = TRUE)) { - first <- seq(0.01,0.1,by=0.01) - firstre <- calc(nn=n,rr=r,g=first) - if(!any(firstre>p & firstre<1,na.rm = TRUE)) { - first <- seq(0.001,0.01,by=0.001) - firstre <- calc(nn=n,rr=r,g=first) - if(!any(firstre>p & firstre<1,na.rm = TRUE)) { - first <- seq(0.0001,0.001,by=0.0001) - firstre <- calc(nn=n,rr=r,g=first) - if(!any(firstre>p & firstre<1,na.rm = TRUE)) { - first <- seq(0.00001,0.0001,by=0.00001) - firstre <- calc(nn=n,rr=r,g=first) - if(!any(firstre>p & firstre<1,na.rm = TRUE)) { - stop("check calculation: specified range (0.00001-0.9) does not cover g") - } - } - } - } - } - index1 <- which(firstre-p>0 & firstre<1) - if(length(index1)>=1) { - index1 <- index1[length(index1)] - lr <- first[c(index1,index1+1)] - } else { - index1 <- which(firstre-p<0 & firstre>0) - index1 <- index1[1] - lr <- first[c(index1-1,index1)] - } - } else { - lr <- c(init-10^floor(log(init,10)),init+10^floor(log(init,10))) - } - - if(lr[1]>lr[2]) lr <- c(lr[2],lr[1]) - - lr <- c(lr[1],(lr[1]+lr[2])/2,lr[2]) - pnew <- p+0.05 - step <- 0 - while(abs(p-pnew)>tol & step < 1000) { - step <- step+1 - temp <- calc(nn=n,rr=r,g=lr) - pnew <- temp[2] - if(abs(p-pnew)>tol) { - if(temp[2]

0){ - lr <- c(lr[1],(lr[1]+lr[2])/2,lr[2]) - } else { - lr <- c(lr[2],(lr[2]+lr[3])/2,lr[3]) - } - } - } - return(lr[2]) -} +#' harmonic analysis test: g-value calculation +#' +#' This function calculates the g-value for the harmonic analysis test developed by R.A. Fisher (1929). +#' Harmonic analysis refers to Fast Fourier Transform (FFT) results. +#' Specifically, g is the proportion (squared modulus of one frequency divided by the sum of all squared moduli). +#' In order for g to be statistically significant in the harmonic analysis test, it needs to be at least g-value +#' at significance level \eqn{\alpha}. Please note that for the rth largest frequency, if any of the previous (r-1) +#' frequencies is not significant, then the rth largest frequency is also non-significant. +#' +#' @param n the total number of frequencies in FFT results. +#' @param r the modulus of the tested frequency is ranked as the rth largest among all frequencies. +#' @param p the FFT result of the tested frequency expressed as the squared modulus divided by the sum of the squared moduli by all frequencies (proportion: m_r^2/(m_1^2+...+m_n^2)). +#' @param tol the tolerance level during calculation. The default is 10^-7. +#' @param init the crude estimate for g-value if known. It is not called to calculate usual g-values. +#' @keywords harmonic test +#' +#' @return The g-value calculated by the harmonic test. + +#' @examples +#' gharmonic(n=100,r=1,p=0.05) +#' +#' @seealso \code{\link{pharmonic}} +#' +#' @references Fisher, R. A. (1929). Tests of significance in harmonic analysis. Proceedings of the Royal Society of London. Series A, 125(796), 54-59. +#' +#' @export + +gharmonic <- function(n,r,p,tol=10^-7,init=NULL) { + calc <- function(nn,rr,g) { + re <- 0 + k <- ifelse(floor(1/g)==1/g,1/g-1,floor(1/g)) + if(k0) { + temp <- lfactorial(nn)-lfactorial(rr-1)+(nn-1)*log(1-j*g)-log(j)-lfactorial(nn-j)-lfactorial(j-r) + if(exp(temp)==0) break + re <- re+(-1)^(j-rr)*exp(temp) + } else { + temp <- lfactorial(nn)-lfactorial(rr-1)+(nn-1)*log(j*g-1)-log(j)-lfactorial(nn-j)-lfactorial(j-r) + if(exp(temp)==0) break + re <- re+(-1)^(j-rr+nn-1)*exp(temp) + } + } + return(re) + } + calc <- Vectorize(calc,vectorize.args ="g") + #<70: 0.1-0.9 + #<983: 0.01-0.1 + if(is.null(init)) { + first <- seq(0.1,0.9,by=0.1) + firstre <- calc(nn=n,rr=r,g=first) + if(!any(firstre>p & firstre<1,na.rm = TRUE)) { + first <- seq(0.01,0.1,by=0.01) + firstre <- calc(nn=n,rr=r,g=first) + if(!any(firstre>p & firstre<1,na.rm = TRUE)) { + first <- seq(0.001,0.01,by=0.001) + firstre <- calc(nn=n,rr=r,g=first) + if(!any(firstre>p & firstre<1,na.rm = TRUE)) { + first <- seq(0.0001,0.001,by=0.0001) + firstre <- calc(nn=n,rr=r,g=first) + if(!any(firstre>p & firstre<1,na.rm = TRUE)) { + first <- seq(0.00001,0.0001,by=0.00001) + firstre <- calc(nn=n,rr=r,g=first) + if(!any(firstre>p & firstre<1,na.rm = TRUE)) { + stop("check calculation: specified range (0.00001-0.9) does not cover g") + } + } + } + } + } + index1 <- which(firstre-p>0 & firstre<1) + if(length(index1)>=1) { + index1 <- index1[length(index1)] + lr <- first[c(index1,index1+1)] + } else { + index1 <- which(firstre-p<0 & firstre>0) + index1 <- index1[1] + lr <- first[c(index1-1,index1)] + } + } else { + lr <- c(init-10^floor(log(init,10)),init+10^floor(log(init,10))) + } + + if(lr[1]>lr[2]) lr <- c(lr[2],lr[1]) + + lr <- c(lr[1],(lr[1]+lr[2])/2,lr[2]) + pnew <- p+0.05 + step <- 0 + while(abs(p-pnew)>tol & step < 1000) { + step <- step+1 + temp <- calc(nn=n,rr=r,g=lr) + pnew <- temp[2] + if(abs(p-pnew)>tol) { + if(temp[2]

0){ + lr <- c(lr[1],(lr[1]+lr[2])/2,lr[2]) + } else { + lr <- c(lr[2],(lr[2]+lr[3])/2,lr[3]) + } + } + } + return(lr[2]) +} diff --git a/R/ind_plot.R b/R/ind_plot.R index d80059e..cc34ed6 100644 --- a/R/ind_plot.R +++ b/R/ind_plot.R @@ -1,101 +1,101 @@ -#' trelliscope auxillary function -#' -#' plot individual observations. -#' -#' @param id individual ID. -#' @param act_ind individual mean observation. -#' @param act_all global mean observation. -#' @param band smoothing band paramter. -#' -#' @keywords internal -#' -#' @return individual to day observations. - -#### individual plot -ind_plot <- function(id,act_ind,act_all,band,ori=TRUE,lw=FALSE) { - timeOfDay <- function(num) { - hour <- floor(num/60) - min <- num - 60*hour - min <- paste(0,min,sep="") - min <- substr(min,nchar(min)-1,nchar(min)) - return(paste(hour,":",min,sep="")) - } - timeOfDay <- Vectorize(timeOfDay,vectorize.args = "num") - xtime <- timeOfDay(1:1440) - - len <- length(act_ind) - xseq <- seq(24/len,24,by=24/len) - d <- band*len - lws <- function(x,ff) lowess(x[c((len-d+1):len,1:len,1:d)],f=ff,iter=1)$y[(d+1):(d+len)] - vec_ind_s <- lws(act_ind,band) - vec_ind <- act_ind - vec_all_s <- lws(act_all,band) - delta <- act_ind-act_all - df <- data.frame(Time=xtime,Individual_Mean=act_ind,Population_Mean=act_all,Delta=delta) - title_name <- paste("ID = ",id, sep="") - - if(lw==FALSE & ori==TRUE) { - rbokeh::figure(title=title_name,legend_location = "top_left", - width=600, height=300, - #xlab="Time of the Day", ylab="Activity Count", - xlim=c(0,24), ylim=c(0, max(c(act_ind,act_all))+50), - tools=c("pan","wheel_zoom","box_zoom","resize", - "box_select","reset","save","help")) %>% - - ly_lines(xseq,vec_ind,width=1,color="blue",alpha=0.2,legend="ind mean") %>% - ly_points(xseq,vec_ind_s,color="white",size=2,alpha=0, - hover=df,legend=FALSE) %>% - ly_lines(xseq,act_all,color="red",width=2,legend="global mean") %>% - x_axis(label="Time of the Day", - desired_num_ticks = 10, num_minor_ticks= 2) %>% - y_axis(label="Activity Count") %>% - theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", - axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% - theme_title(text_font_size="12pt",text_font="Helvetica") %>% - theme_legend(label_text_font_size = "8pt") - } else if(lw==TRUE & ori==TRUE) { - rbokeh::figure(title=title_name,legend_location = "top_left", - width=600, height=300, - #xlab="Time of the Day", ylab="Activity Count", - xlim=c(0,24), ylim=c(0, max(c(act_ind,act_all))+50), - tools=c("pan","wheel_zoom","box_zoom","resize", - "box_select","reset","save","help")) %>% - - ly_lines(xseq,vec_ind,width=1,color="blue",alpha=0.2,legend="ind mean") %>% - ly_points(xseq,vec_ind_s,color="white",size=2,alpha=0, - hover=df,legend=FALSE) %>% - ly_lines(xseq,vec_ind_s,width=2,color="blue",legend="ind mean lowess") %>% - ly_lines(xseq,act_all,color="red",width=2,legend="global mean") %>% - ly_lines(xseq,vec_all_s,color="yellow",width=2,legend="global mean lowess") %>% - x_axis(label="Time of the Day", - desired_num_ticks = 10, num_minor_ticks= 2) %>% - y_axis(label="Activity Count") %>% - theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", - axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% - theme_title(text_font_size="12pt",text_font="Helvetica") %>% - theme_legend(label_text_font_size = "8pt") - } else if(lw==TRUE & ori==FALSE) { - rbokeh::figure(title=title_name,legend_location = "top_left", - width=600, height=300, - #xlab="Time of the Day", ylab="Activity Count", - xlim=c(0,24), ylim=c(0, max(c(act_ind,act_all))+50), - tools=c("pan","wheel_zoom","box_zoom","resize", - "box_select","reset","save","help")) %>% - - ly_lines(xseq,vec_ind,width=1,color="blue",alpha=0.2,legend="ind mean") %>% - ly_points(xseq,vec_ind_s,color="white",size=2,alpha=0, - hover=df,legend=FALSE) %>% - ly_lines(xseq,vec_ind_s,width=2,color="blue",legend="ind mean lowess") %>% - ly_lines(xseq,vec_all_s,color="yellow",width=2,legend="global mean lowess") %>% - x_axis(label="Time of the Day", - desired_num_ticks = 10, num_minor_ticks= 2) %>% - y_axis(label="Activity Count") %>% - theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", - axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% - theme_title(text_font_size="12pt",text_font="Helvetica") %>% - theme_legend(label_text_font_size = "8pt") - } else { - stop("Error: at least plot one of the original and lowess!") - } -} - +#' trelliscope auxillary function +#' +#' plot individual observations. +#' +#' @param id individual ID. +#' @param act_ind individual mean observation. +#' @param act_all global mean observation. +#' @param band smoothing band paramter. +#' +#' @keywords internal +#' +#' @return individual to day observations. + +#### individual plot +ind_plot <- function(id,act_ind,act_all,band,ori=TRUE,lw=FALSE) { + timeOfDay <- function(num) { + hour <- floor(num/60) + min <- num - 60*hour + min <- paste(0,min,sep="") + min <- substr(min,nchar(min)-1,nchar(min)) + return(paste(hour,":",min,sep="")) + } + timeOfDay <- Vectorize(timeOfDay,vectorize.args = "num") + xtime <- timeOfDay(1:1440) + + len <- length(act_ind) + xseq <- seq(24/len,24,by=24/len) + d <- band*len + lws <- function(x,ff) lowess(x[c((len-d+1):len,1:len,1:d)],f=ff,iter=1)$y[(d+1):(d+len)] + vec_ind_s <- lws(act_ind,band) + vec_ind <- act_ind + vec_all_s <- lws(act_all,band) + delta <- act_ind-act_all + df <- data.frame(Time=xtime,Individual_Mean=act_ind,Population_Mean=act_all,Delta=delta) + title_name <- paste("ID = ",id, sep="") + + if(lw==FALSE & ori==TRUE) { + rbokeh::figure(title=title_name,legend_location = "top_left", + width=600, height=300, + #xlab="Time of the Day", ylab="Activity Count", + xlim=c(0,24), ylim=c(0, max(c(act_ind,act_all))+50), + tools=c("pan","wheel_zoom","box_zoom","resize", + "box_select","reset","save","help")) %>% + + ly_lines(xseq,vec_ind,width=1,color="blue",alpha=0.2,legend="ind mean") %>% + ly_points(xseq,vec_ind_s,color="white",size=2,alpha=0, + hover=df,legend=FALSE) %>% + ly_lines(xseq,act_all,color="red",width=2,legend="global mean") %>% + x_axis(label="Time of the Day", + desired_num_ticks = 10, num_minor_ticks= 2) %>% + y_axis(label="Activity Count") %>% + theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", + axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% + theme_title(text_font_size="12pt",text_font="Helvetica") %>% + theme_legend(label_text_font_size = "8pt") + } else if(lw==TRUE & ori==TRUE) { + rbokeh::figure(title=title_name,legend_location = "top_left", + width=600, height=300, + #xlab="Time of the Day", ylab="Activity Count", + xlim=c(0,24), ylim=c(0, max(c(act_ind,act_all))+50), + tools=c("pan","wheel_zoom","box_zoom","resize", + "box_select","reset","save","help")) %>% + + ly_lines(xseq,vec_ind,width=1,color="blue",alpha=0.2,legend="ind mean") %>% + ly_points(xseq,vec_ind_s,color="white",size=2,alpha=0, + hover=df,legend=FALSE) %>% + ly_lines(xseq,vec_ind_s,width=2,color="blue",legend="ind mean lowess") %>% + ly_lines(xseq,act_all,color="red",width=2,legend="global mean") %>% + ly_lines(xseq,vec_all_s,color="yellow",width=2,legend="global mean lowess") %>% + x_axis(label="Time of the Day", + desired_num_ticks = 10, num_minor_ticks= 2) %>% + y_axis(label="Activity Count") %>% + theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", + axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% + theme_title(text_font_size="12pt",text_font="Helvetica") %>% + theme_legend(label_text_font_size = "8pt") + } else if(lw==TRUE & ori==FALSE) { + rbokeh::figure(title=title_name,legend_location = "top_left", + width=600, height=300, + #xlab="Time of the Day", ylab="Activity Count", + xlim=c(0,24), ylim=c(0, max(c(act_ind,act_all))+50), + tools=c("pan","wheel_zoom","box_zoom","resize", + "box_select","reset","save","help")) %>% + + ly_lines(xseq,vec_ind,width=1,color="blue",alpha=0.2,legend="ind mean") %>% + ly_points(xseq,vec_ind_s,color="white",size=2,alpha=0, + hover=df,legend=FALSE) %>% + ly_lines(xseq,vec_ind_s,width=2,color="blue",legend="ind mean lowess") %>% + ly_lines(xseq,vec_all_s,color="yellow",width=2,legend="global mean lowess") %>% + x_axis(label="Time of the Day", + desired_num_ticks = 10, num_minor_ticks= 2) %>% + y_axis(label="Activity Count") %>% + theme_axis(which=c("x","y"),axis_label_text_font = "Helvetica", + axis_label_text_font_size = "12pt", axis_label_text_font_style = "normal") %>% + theme_title(text_font_size="12pt",text_font="Helvetica") %>% + theme_legend(label_text_font_size = "8pt") + } else { + stop("Error: at least plot one of the original and lowess!") + } +} + diff --git a/R/ind_to_day.R b/R/ind_to_day.R index ce4af8c..35d3e42 100644 --- a/R/ind_to_day.R +++ b/R/ind_to_day.R @@ -1,35 +1,35 @@ -#' trelliscope auxillary function -#' -#' Take individual characteristics and map it to individual-day observations. For example, individual A/B is 40/45 -#' years old and has 4/2 day observations in dataset D. Therefore, mapping age to dataset D generates -#' 40, 40, 40, 40, 45, 45 for day 1,2,3,4/1,2 observation of A/B respectively. -#' -#' @param x covariates data for individuals to be merged. -#' @param df day observation dataset. -#' -#' @keywords internal -#' -#' @return individual to day observations. - -ind_to_day <- function(x,df) { - if(is.list(x)) { - temp <- unlist(lapply(split(df$ID,df$ID),length)) - temp2 <- do.call("cbind",x) - re <- list() - for(i in 1:length(temp)) { - re[[i]] <- temp2[,rep(i,temp[i])] - } - re <- do.call("cbind",re) - temp <- list() - for (i in 1:ncol(re)) temp[[i]] <- re[,i] - return(temp) - - } else if (is.vector(x)) { - temp <- unlist(lapply(split(df$ID,df$ID),length)) - temp <- cbind(temp,x) - temp <- unlist(apply(temp,1,function(x) rep(x[2],x[1]))) - names(temp) <- NULL - temp <- as.numeric(temp) - return(temp) - } -} +#' trelliscope auxillary function +#' +#' Take individual characteristics and map it to individual-day observations. For example, individual A/B is 40/45 +#' years old and has 4/2 day observations in dataset D. Therefore, mapping age to dataset D generates +#' 40, 40, 40, 40, 45, 45 for day 1,2,3,4/1,2 observation of A/B respectively. +#' +#' @param x covariates data for individuals to be merged. +#' @param df day observation dataset. +#' +#' @keywords internal +#' +#' @return individual to day observations. + +ind_to_day <- function(x,df) { + if(is.list(x)) { + temp <- unlist(lapply(split(df$ID,df$ID),length)) + temp2 <- do.call("cbind",x) + re <- list() + for(i in 1:length(temp)) { + re[[i]] <- temp2[,rep(i,temp[i])] + } + re <- do.call("cbind",re) + temp <- list() + for (i in 1:ncol(re)) temp[[i]] <- re[,i] + return(temp) + + } else if (is.vector(x)) { + temp <- unlist(lapply(split(df$ID,df$ID),length)) + temp <- cbind(temp,x) + temp <- unlist(apply(temp,1,function(x) rep(x[2],x[1]))) + names(temp) <- NULL + temp <- as.numeric(temp) + return(temp) + } +} diff --git a/R/lis3.R b/R/lis3.R index 949e0d2..d44881f 100644 --- a/R/lis3.R +++ b/R/lis3.R @@ -1,20 +1,20 @@ -#' Activity data for three individuals -#' -#' A synthetic data list of three elements, each consisting of an activity matrix for one individual, and -#' each column of an activity matrix is one-day activity observation. Therefore, an activity matrix is -#' \code{nob} by \code{nday}. It is a named list, the name of which is the individual ID. It is for illustration -#' purposes of the functions in the package \code{PML} only. -#' -#' @docType data -#' -#' @usage data(lis3) -#' -#' @format An object of class \code{list}. -#' -#' @keywords datasets -#' -#' @examples -#' data(lis3) -#' pa3 <- form(lis3) -#' +#' Activity data for three individuals +#' +#' A synthetic data list of three elements, each consisting of an activity matrix for one individual, and +#' each column of an activity matrix is one-day activity observation. Therefore, an activity matrix is +#' \code{nob} by \code{nday}. It is a named list, the name of which is the individual ID. It is for illustration +#' purposes of the functions in the package \code{PML} only. +#' +#' @docType data +#' +#' @usage data(lis3) +#' +#' @format An object of class \code{list}. +#' +#' @keywords datasets +#' +#' @examples +#' data(lis3) +#' pa3 <- form(lis3) +#' #' @seealso \code{\link{form}}, \code{\link{pa3}} \ No newline at end of file diff --git a/R/pa3.R b/R/pa3.R index 3206923..8e9be7d 100644 --- a/R/pa3.R +++ b/R/pa3.R @@ -1,20 +1,20 @@ -#' Activity tbl_df data for three individuals -#' -#' A synthetic data in tbl_df format. It has 13 observations for 3 individuals, and the variables are "ID", -#' "ID_Nday" (the ith day observation for an individual), and activity. The activity variable is an embedded -#' list with each element consisting of a vector of one-day observation. It is for illustration purposes of -#' the functions in the package \code{PML} only. -#' -#' @docType data -#' -#' @usage data(pa3) -#' -#' @format An object of classes \code{tbl_df}, \code{tbl}, \code{data.frame}. -#' -#' @keywords datasets -#' -#' @examples -#' data(pa3) -#' re <- bandSelect(df=pa3,Nlength=1440*3,Nlambda=100,alpha=1,Ntop=5,cross=FALSE,Ncross=NULL,plot=TRUE) -#' +#' Activity tbl_df data for three individuals +#' +#' A synthetic data in tbl_df format. It has 13 observations for 3 individuals, and the variables are "ID", +#' "ID_Nday" (the ith day observation for an individual), and activity. The activity variable is an embedded +#' list with each element consisting of a vector of one-day observation. It is for illustration purposes of +#' the functions in the package \code{PML} only. +#' +#' @docType data +#' +#' @usage data(pa3) +#' +#' @format An object of classes \code{tbl_df}, \code{tbl}, \code{data.frame}. +#' +#' @keywords datasets +#' +#' @examples +#' data(pa3) +#' re <- bandSelect(df=pa3,Nlength=1440*3,Nlambda=100,alpha=1,Ntop=5,cross=FALSE,Ncross=NULL,plot=TRUE) +#' #' @seealso \code{\link{bandSelect}} \ No newline at end of file diff --git a/R/pharmonic.R b/R/pharmonic.R index 2938fa4..cf1c364 100644 --- a/R/pharmonic.R +++ b/R/pharmonic.R @@ -1,45 +1,45 @@ -#' Harmonic analysis test: p-value calculation -#' -#' This function calculates the p-value for the harmonic analysis test developed by R.A. Fisher (1929). Harmonic analysis specifically refers to Fast Fourier Transform (FFT) results. -#' -#' @param n the total number of frequencies in FFT results -#' @param r the modulus of the tested frequency is ranked as the rth largest among all frequencies -#' @param g the FFT result of the tested frequency expressed as the squared modulus divided by the sum of the squared moduli by all frequencies (proportion: m_r^2/(m_1^2+...+m_n^2)). -#' -#' @keywords harmonic test -#' -#' @return The p-value calculated by the harmonic test. - -#' @examples -#' pharmonic(n=100,r=2,g=0.1) -#' -#' @seealso \code{\link{gharmonic}} -#' -#' @references Fisher, R. A. (1929). Tests of significance in harmonic analysis. Proceedings of the Royal Society of London. Series A, 125(796), 54-59. -#' -#' @export - -pharmonic <- Vectorize(function(n,r,g) { - re <- 0 - k <- ifelse(floor(1/g)==1/g,1/g-1,floor(1/g)) - term0 <- 0 - if(k0) { - temp <- lfactorial(n)-lfactorial(r-1)+(n-1)*log(1-j*g)-log(j)-lfactorial(n-j)-lfactorial(j-r) - if(exp(temp)==0) { - term0 <- term0+1 - break - } - re <- re+(-1)^(j-r)*exp(temp) - } else { - temp <- lfactorial(n)-lfactorial(r-1)+(n-1)*log(j*g-1)-log(j)-lfactorial(n-j)-lfactorial(j-r) - if(exp(temp)==0) { - term0 <- term0+1 - break - } - re <- re+(-1)^(j-r+n-1)*exp(temp) - } - } - return(re) - },vec="g") +#' Harmonic analysis test: p-value calculation +#' +#' This function calculates the p-value for the harmonic analysis test developed by R.A. Fisher (1929). Harmonic analysis specifically refers to Fast Fourier Transform (FFT) results. +#' +#' @param n the total number of frequencies in FFT results +#' @param r the modulus of the tested frequency is ranked as the rth largest among all frequencies +#' @param g the FFT result of the tested frequency expressed as the squared modulus divided by the sum of the squared moduli by all frequencies (proportion: m_r^2/(m_1^2+...+m_n^2)). +#' +#' @keywords harmonic test +#' +#' @return The p-value calculated by the harmonic test. + +#' @examples +#' pharmonic(n=100,r=2,g=0.1) +#' +#' @seealso \code{\link{gharmonic}} +#' +#' @references Fisher, R. A. (1929). Tests of significance in harmonic analysis. Proceedings of the Royal Society of London. Series A, 125(796), 54-59. +#' +#' @export + +pharmonic <- Vectorize(function(n,r,g) { + re <- 0 + k <- ifelse(floor(1/g)==1/g,1/g-1,floor(1/g)) + term0 <- 0 + if(k0) { + temp <- lfactorial(n)-lfactorial(r-1)+(n-1)*log(1-j*g)-log(j)-lfactorial(n-j)-lfactorial(j-r) + if(exp(temp)==0) { + term0 <- term0+1 + break + } + re <- re+(-1)^(j-r)*exp(temp) + } else { + temp <- lfactorial(n)-lfactorial(r-1)+(n-1)*log(j*g-1)-log(j)-lfactorial(n-j)-lfactorial(j-r) + if(exp(temp)==0) { + term0 <- term0+1 + break + } + re <- re+(-1)^(j-r+n-1)*exp(temp) + } + } + return(re) + },vec="g") diff --git a/R/test.harmonic.R b/R/test.harmonic.R index 19449f8..16b3ad9 100644 --- a/R/test.harmonic.R +++ b/R/test.harmonic.R @@ -1,79 +1,79 @@ -#' Harmonic analysis test for Fast Fourier Transform -#' -#' This function conducts harmonic test sequentially based on observations or Fast Fourier Transform (FFT) -#' results. -#' -#' @param ob Either the original observation or FFT results. See parameter \code{fft}. -#' @param p The p-value to be considered statistically significant. -#' @param fft If TRUE, \code{ob} is FFT results, with the first column frequencies and the second column signals in standardized proportions; if FALSE, \code{ob} is a vector of the original observation. The default is FALSE. -#' @param maxfreq To conduct test on at most \code{maxfreq} frequencies. The default is 10. -#' -#' @keywords harmonic test -#' -#' @return A list of two elements: -#' @return \item{sig}{The significant frequencies plus the first insignificant frequency.} -#' @return \item{fft}{The FFT results expressed in standardized proportions.} -#' @examples -#' data(pa3) -#' -#' #### test on individuals -#' ob <- do.call("c",pa3$activity[1:4]) -#' re <- test.harmonic(ob,p=0.05/(length(ob)-1)/2) -#' re$sig;head(re$fft) ## no harmonic is significant -#' ob2 <- do.call("c",pa3$activity[11:13]) -#' re2 <- test.harmonic(ob2,p=0.05/(length(ob2)-1)/2) -#' re2$sig;head(re2$fft) ## 3 significant harmonics -#' -#' #### test on the population average -#' re0 <- bandSelect(df=pa3,Nlength=1440*3,Nlambda=100,alpha=1,Ntop=3, -#' cross=FALSE,Ncross=NULL,plot=TRUE) -#' freq <- data.frame(Frequency=re0$freq,Proportion=colMeans(re0$xprop)) -#' re3 <- test.harmonic(freq,p=0.05/nrow(freq),fft=TRUE) -#' print(re3$sig,digits=3,row.names=FALSE) -#' -#' @seealso \code{\link{pharmonic}} -#' -#' @references Fisher, R. A. (1929). Tests of significance in harmonic analysis. Proceedings of the Royal Society of London. Series A, 125(796), 54-59. -#' -#' @export - -test.harmonic <- function(ob,p,fft=FALSE,maxfreq=10) { - if(fft==TRUE) { - fre <- ob[order(ob[,2],decreasing=TRUE),] - } else { - f <- fft(ob) - fre <- fft.harmonic(f) - fre <- fre[,c(3,4)] - } - names(fre) <- c("freq","prop(g)") - n <- nrow(fre) - re <- data.frame(freq=0,prop=0,p=0) - for(i in 1:maxfreq) { - r <- i - p2 <- pharmonic(n,r,fre[i,2]) - if(p2maxfreq) { - message(paste("P-value selected too large and too many frequencies significant ( > ,", maxfreq, ",max freq reached), consider Bonferonni correction?",sep="")) - break - } - } - } else { - if (i==1) { - re[1,] <- c(fre[i,c(1,2)],p2) - } else { - re <- rbind(re,re[1,]) - re[i,] <- c(fre[i,c(1,2)],p2) - } - break - } - } - pp <- rep(p,nrow(re)) - re <- cbind(re,pp) - colnames(re) <- c("frequency","prop (g)","p-value","p-threshold") - return(list(sig=re,fft=fre)) -} +#' Harmonic analysis test for Fast Fourier Transform +#' +#' This function conducts harmonic test sequentially based on observations or Fast Fourier Transform (FFT) +#' results. +#' +#' @param ob Either the original observation or FFT results. See parameter \code{fft}. +#' @param p The p-value to be considered statistically significant. +#' @param fft If TRUE, \code{ob} is FFT results, with the first column frequencies and the second column signals in standardized proportions; if FALSE, \code{ob} is a vector of the original observation. The default is FALSE. +#' @param maxfreq To conduct test on at most \code{maxfreq} frequencies. The default is 10. +#' +#' @keywords harmonic test +#' +#' @return A list of two elements: +#' @return \item{sig}{The significant frequencies plus the first insignificant frequency.} +#' @return \item{fft}{The FFT results expressed in standardized proportions.} +#' @examples +#' data(pa3) +#' +#' #### test on individuals +#' ob <- do.call("c",pa3$activity[1:4]) +#' re <- test.harmonic(ob,p=0.05/(length(ob)-1)/2) +#' re$sig;head(re$fft) ## no harmonic is significant +#' ob2 <- do.call("c",pa3$activity[11:13]) +#' re2 <- test.harmonic(ob2,p=0.05/(length(ob2)-1)/2) +#' re2$sig;head(re2$fft) ## 3 significant harmonics +#' +#' #### test on the population average +#' re0 <- bandSelect(df=pa3,Nlength=1440*3,Nlambda=100,alpha=1,Ntop=3, +#' cross=FALSE,Ncross=NULL,plot=TRUE) +#' freq <- data.frame(Frequency=re0$freq,Proportion=colMeans(re0$xprop)) +#' re3 <- test.harmonic(freq,p=0.05/nrow(freq),fft=TRUE) +#' print(re3$sig,digits=3,row.names=FALSE) +#' +#' @seealso \code{\link{pharmonic}} +#' +#' @references Fisher, R. A. (1929). Tests of significance in harmonic analysis. Proceedings of the Royal Society of London. Series A, 125(796), 54-59. +#' +#' @export + +test.harmonic <- function(ob,p,fft=FALSE,maxfreq=10) { + if(fft==TRUE) { + fre <- ob[order(ob[,2],decreasing=TRUE),] + } else { + f <- fft(ob) + fre <- fft.harmonic(f) + fre <- fre[,c(3,4)] + } + names(fre) <- c("freq","prop(g)") + n <- nrow(fre) + re <- data.frame(freq=0,prop=0,p=0) + for(i in 1:maxfreq) { + r <- i + p2 <- pharmonic(n,r,fre[i,2]) + if(p2maxfreq) { + message(paste("P-value selected too large and too many frequencies significant ( > ,", maxfreq, ",max freq reached), consider Bonferonni correction?",sep="")) + break + } + } + } else { + if (i==1) { + re[1,] <- c(fre[i,c(1,2)],p2) + } else { + re <- rbind(re,re[1,]) + re[i,] <- c(fre[i,c(1,2)],p2) + } + break + } + } + pp <- rep(p,nrow(re)) + re <- cbind(re,pp) + colnames(re) <- c("frequency","prop (g)","p-value","p-threshold") + return(list(sig=re,fft=fre)) +} diff --git a/R/tre.R b/R/tre.R index e68fe6f..968472f 100644 --- a/R/tre.R +++ b/R/tre.R @@ -1,179 +1,176 @@ -#' Trelliscope Visualization for Accelerometer Data -#' -#' This function generates the data frame necessary for trelliscope visualization. -#' -#' @importFrom"dplyr" -#' '%>%' -#' @importFrom"dplyr" -#' mutate -#' @importFrom"trelliscopejs" -#' trelliscope -#' @importFrom"trelliscopejs" -#' pmap_plot -#' @import"rbokeh" -#' -#' @param lis the list of activity data, with each element corresponding to the observation by one individual and the name of each element coresponding to the individual id. Specifically, each element is a \code{nob} by \code{nday} matrix, where each column is an observation by day. -#' @param id a vector of id names corresponding to the \code{lis} activity data. -#' @param varlis optional data frame to be merged to activity data, and the covariates are of interest for plotting to see activity differences. The first variables needs to be "ID". -#' @param smband smoothing parameter for plotting smoothed activity data. the default is 1/12 (see function \code{lowess}). -#' @param maxday maxday the maximal number of days per individual in the observation, used to check the data format. The default is 14. -#' @param plot.ind whether to plot individual mean activity plots. If not, plot day activity plots. The default is TRUE. -#' @param plot.ori whether to plot the original activity curves (tend to have large variations). The default is TRUE. -#' @param plot.sm whether to plot lowess of the activity curves. The default is TRUE. -#' @param plot.tre whether to generate trelliscope plots. If so, no data will be returned; if not, a data frame will be returned containing all information including trelliscope panels. To generate trelliscope based on the data, one needs to set all activity list columns to NULL. The default is FALSE. -#' @param plot.tre.path If plot.tre is TRUE, then plot.tre.path specifies the path to generate trelliscope files. The default is current working directory. -#' @keywords trelliscope -#' -#' @return The data frame including activity, filtering stats, optional covariates, and trelliscope panels. (No data frame will be returned if plot.tre is TRUE.) - -#' @examples -#' \dontrun{ -#' data(lis3) -#' data(var3) -#' -#' #### individual mean activity plot: return a dataset with trelliscope panels -#' tre.ind <- tre(lis3,varlis=var3) -#' tre.ind$activity_ind <- tre.ind$activity_all <- NULL -#' trelliscopejs::trelliscope(tre.ind,name = "Individual Mean Activity Plot", -#' nrow = 2, ncol = 2,path=tempdir()) -#' -#' #### day activity plot: directly generating trelliscope visualization -#' tre(lis3,plot.ind=FALSE,plot.ori=FALSE,plot.tre=TRUE,plot.tre.path=tempdir()) -#' } -#' -#' @seealso \code{\link{form}} -#' -#' @export - -tre <- function(lis,id=NULL,varlis=NULL,smband=1/12,maxday=14,plot.ind=TRUE,plot.ori=TRUE,plot.sm=TRUE,plot.tre=FALSE,plot.tre.path=NULL) { - checkformat <- do.call("c",lapply(lis,ncol)) - if(length(checkformat)!=length(lis)) { - stop("Contain empty matrix in the data list.") - } - if(max(checkformat)>maxday) { - stop(paste("Maxday ",maxday," reached: data list format may be wrong / consider change maxday.",sep="")) - } - - #### ID info - ID <- list() - if(!is.null(id)) { - if(length(id)!=length(lis)) { - stop("The length of the ID vector is not the same as the length of the data list.") - } - for (i in 1:length(id)) ID[[i]] <- rep(names(id[i]),ncol(lis[[i]])) - } else { - if(is.null(names(lis))) { - stop("Names of the data list are null: ID should be supplied.") - } - for (i in 1:length(lis)) ID[[i]] <- rep(names(lis[i]),ncol(lis[[i]])) - } - - #### ID_Nday info - ID_Nday <- activity <- activity_ind <- activity_max <- activity_all <- NULL #required to avoid NOTE in CMD check - act <- data.frame(ID=as.numeric(unlist(ID)), - ID_Nday=unlist(lapply(ID,function(x) seq(1,length(x))))) - act <- dplyr::tbl_df(act) - act <- dplyr::group_by(act,ID,ID_Nday) - act <- tidyr::nest(act) - - #### activity - liscol <- do.call("cbind",lis) - liscol2 <- list() - for(i in 1:ncol(liscol)) liscol2[[i]] <- liscol[,i] - act <- dplyr::mutate(act,activity=liscol2);rm(liscol) - act$data <- NULL - ## smoothed activity - act <- dplyr::mutate(act,activity_sm=lapply(act$activity,function(x) round(lowess(x[c((length(x)*(1-smband)+1):length(x),1:length(x),1:(length(x)*smband))],f=smband)$y[(length(x)*smband+1):(length(x)*(1+smband))],1))) - ## individual mean activity - act <- dplyr::mutate(act,activity_ind=lapply(ind_to_day(lapply(lis,rowMeans),act),function(x) round(x,1))) - ## global mean activity - mean_global <- round(rowMeans(do.call("cbind",act$activity_ind)),1) - temp <- list() - for(i in 1:nrow(act)) temp[[i]] <- mean_global - act <- dplyr::mutate(act,activity_all=temp) - rm(temp,mean_global) - - #### filter stats - ##count_mean - act$count_mean <- unlist(lapply(liscol2,mean)) - ##max number of consecutive zeros - act$zero_consecmax <- unlist(lapply(act$activity,function(x) { - re <- rle(x) - if(sum(re$values==0)==0) { - return(0) ##no zeros at all - } else { - return(max(re$lengths[re$values==0])) - } - })) - ##total number of zeros - act$zero_Nmax <- unlist(lapply(act$activity,function(x) sum(x==0))) - - ####plot individuals (average over days) or days - if(plot.ind==TRUE) { - act <- act[!duplicated(act$ID),] - act$activity <- act$activity_sm <- NULL - - ### merge with other datasets - if(!is.null(varlis)) { - deltadf <- data.frame(ID=unique(act$ID)) - deltadf <- merge(deltadf,varlis,by="ID",all.x=TRUE) - act <- cbind(act,deltadf[,-1]);rm(deltadf) - } - - ### generate plot: ind vs global - ptm <- Sys.time() - message("Generating trelliscope individual plots... It may take some time.") - ind <- act[!duplicated(act$ID),] - ind <- dplyr::mutate(act,panel = trelliscopejs::pmap_plot(list(ind$ID,ind$activity_ind, - ind$activity_all,smband,plot.ori,plot.sm), ind_plot)) - message(paste("Total time: ",round(difftime(Sys.time(),ptm,units="mins")[[1]],2)," mins",sep="")) - - ## trelliscope plot - if(plot.tre==TRUE) { - ind$activity_ind <- ind$activity_all <- NULL - if(is.null(plot.tre.path)) { - trelliscopejs::trelliscope(ind,name = "Individual Mean Activity Plot", nrow = 2, ncol = 2, - path=getwd()) - } else { - trelliscopejs::trelliscope(ind,name = "Individual Mean Activity Plot", nrow = 2, ncol = 2, - path=plot.tre.path) - } - } else { - return(ind) - } - - } else { - ## merge with other datasets - if(!is.null(varlis)) { - deltadf <- data.frame(ID=unique(act$ID)) - deltadf <- merge(deltadf,varlis,by="ID",all.x=TRUE) - deltadf2 <- apply(deltadf,2,function(y) ind_to_day(x=y,df=act)) - act <- cbind(act,deltadf2[,-1]);rm(deltadf,deltadf2) - } - ## auxillary information for plotting: y-axis activity max - act <- dplyr::mutate(act,activity_max = ind_to_day(unlist(lapply(split(act$activity_sm,act$ID),function(x) - max(unlist(lapply(x,max))))),act)) - - ## generate plot: day observation - ptm <- Sys.time() - message("Generating trelliscope activity day plots... It may take some time.") - act <- dplyr::mutate(act,panel = trelliscopejs::pmap_plot(list(id=ID,id_Nday=ID_Nday,act_ori=activity,act_ind=activity_ind, - act_all=activity_all,act_max=activity_max,band=smband,ori=plot.ori,lw=plot.sm), act_plot)) - message(paste("Total time: ",round(difftime(Sys.time(),ptm,units="mins")[[1]],3)," mins",sep="")) - #check memory: format(object.size(act),units="Mb",standard="legacy") - - ## trelliscope plot - if(plot.tre==TRUE) { - act$activity <- act$activity_sm <- act$activity_ind <- act$activity_all <- NULL - if(is.null(plot.tre.path)) { - trelliscopejs::trelliscope(act,name = "Daily Activity Plot", nrow = 2, ncol = 2, - path=getwd()) - } else { - trelliscopejs::trelliscope(act,name = "Day Activity Plot", nrow = 2, ncol = 2, - path=plot.tre.path) - } - } else { - return(act) - } - } -} +#' Trelliscope Visualization for Accelerometer Data +#' +#' This function generates the data frame necessary for trelliscope visualization. +#' +#' @importFrom"dplyr" +#' '%>%' +#' @importFrom"tibble" +#' add_column +#' @import"rbokeh" +#' +#' @param lis the list of activity data, with each element corresponding to the observation by one individual and the name of each element coresponding to the individual id. Specifically, each element is a \code{nob} by \code{nday} matrix, where each column is an observation by day. +#' @param id a vector of id names corresponding to the \code{lis} activity data. +#' @param varlis optional data frame to be merged to activity data, and the covariates are of interest for plotting to see activity differences. The first variables needs to be "ID". +#' @param smband smoothing parameter for plotting smoothed activity data. the default is 1/12 (see function \code{lowess}). +#' @param maxday maxday the maximal number of days per individual in the observation, used to check the data format. The default is 14. +#' @param plot.ind whether to plot individual mean activity plots. If not, plot day activity plots. The default is TRUE. +#' @param plot.ori whether to plot the original activity curves (tend to have large variations). The default is TRUE. +#' @param plot.sm whether to plot lowess of the activity curves. The default is TRUE. +#' @param plot.tre whether to generate trelliscope plots. If so, no data will be returned; if not, a data frame will be returned containing all information including trelliscope panels. To generate trelliscope based on the data, one needs to set all activity list columns to NULL. The default is FALSE. +#' @param plot.tre.path If plot.tre is TRUE, then plot.tre.path specifies the path to generate trelliscope files. The default is current working directory. +#' @keywords trelliscope +#' +#' @return The data frame including activity, filtering stats, optional covariates, and trelliscope panels. (No data frame will be returned if plot.tre is TRUE.) + +#' @examples +#' data(lis3) +#' data(var3) +#' +#' #### individual mean activity plot: return a dataset with trelliscope panels +#' tre.ind <- tre(lis3,varlis=var3) +#' tre.ind$activity_ind <- tre.ind$activity_all <- NULL +#' +#' @seealso \code{\link{form}} +#' +#' @export + +tre <- function(lis,id=NULL,varlis=NULL,smband=1/12,maxday=14,plot.ind=TRUE,plot.ori=TRUE,plot.sm=TRUE,plot.tre=FALSE,plot.tre.path=NULL) { + checkformat <- do.call("c",lapply(lis,ncol)) + if(length(checkformat)!=length(lis)) { + stop("Contain empty matrix in the data list.") + } + if(max(checkformat)>maxday) { + stop(paste("Maxday ",maxday," reached: data list format may be wrong / consider change maxday.",sep="")) + } + + #### ID info + ID <- list() + if(!is.null(id)) { + if(length(id)!=length(lis)) { + stop("The length of the ID vector is not the same as the length of the data list.") + } + for (i in 1:length(id)) ID[[i]] <- rep(names(id[i]),ncol(lis[[i]])) + } else { + if(is.null(names(lis))) { + stop("Names of the data list are null: ID should be supplied.") + } + for (i in 1:length(lis)) ID[[i]] <- rep(names(lis[i]),ncol(lis[[i]])) + } + + #### ID_Nday info + ID_Nday <- activity <- activity_ind <- activity_max <- activity_all <- NULL #required to avoid NOTE in CMD check + act <- data.frame(ID=as.numeric(unlist(ID)), + ID_Nday=unlist(lapply(ID,function(x) seq(1,length(x))))) + ### merge with other datasets + if(!is.null(varlis)) { + deltadf <- data.frame(ID=act$ID) + deltadf <- merge(deltadf,varlis,by="ID",all.x=TRUE) + act <- cbind(act,deltadf[,-1]);rm(deltadf) + } + ##format + act <- dplyr::tbl_df(act) + act <- dplyr::group_by(act,ID,ID_Nday) + act <- tidyr::nest(act,data=c()) + + #### activity + liscol <- do.call("cbind",lis) + liscol2 <- list() + for(i in 1:ncol(liscol)) liscol2[[i]] <- liscol[,i] + #act <- dplyr::mutate(act,activity=liscol2);rm(liscol) updated function due to changes in the tidyr package + act <- tibble::add_column(act,activity = liscol2);rm(liscol) + act$data <- NULL + ## smoothed activity + #act <- dplyr::mutate(act,activity_sm=lapply(act$activity,function(x) round(lowess(x[c((length(x)*(1-smband)+1):length(x),1:length(x),1:(length(x)*smband))],f=smband)$y[(length(x)*smband+1):(length(x)*(1+smband))],1))) updated function due to changes in the tidyr package + act <- tibble::add_column(act,activity_sm=lapply(act$activity,function(x) round(lowess(x[c((length(x)*(1-smband)+1):length(x),1:length(x),1:(length(x)*smband))],f=smband)$y[(length(x)*smband+1):(length(x)*(1+smband))],1))) + ## individual mean activity + #act <- dplyr::mutate(act,activity_ind=lapply(ind_to_day(lapply(lis,rowMeans),act),function(x) round(x,1))) updated function due to changes in the tidyr package + act <- tibble::add_column(act,activity_ind=lapply(ind_to_day(lapply(lis,rowMeans),act),function(x) round(x,1))) + ## global mean activity + mean_global <- round(rowMeans(do.call("cbind",act$activity_ind)),1) + temp <- list() + for(i in 1:nrow(act)) temp[[i]] <- mean_global + #act <- dplyr::mutate(act,activity_all=temp) updated function due to changes in the tidyr package + act <- tibble::add_column(act,activity_all=temp) + rm(temp,mean_global) + + #### filter stats + ##count_mean + act$count_mean <- unlist(lapply(liscol2,mean)) + ##max number of consecutive zeros + act$zero_consecmax <- unlist(lapply(act$activity,function(x) { + re <- rle(x) + if(sum(re$values==0)==0) { + return(0) ##no zeros at all + } else { + return(max(re$lengths[re$values==0])) + } + })) + ##total number of zeros + act$zero_Nmax <- unlist(lapply(act$activity,function(x) sum(x==0))) + + ####plot individuals (average over days) or days + if(plot.ind==TRUE) { + ### generate plot: ind vs global + #ptm <- Sys.time() + message("No trelliscope individual plots are generated due to the archived trelliscopejs package.") + ind <- act[!duplicated(act$ID),] + ind$activity <- ind$activity_sm <- NULL + #ind <- tibble::add_column(ind,panel = trelliscopejs::pmap_plot(list(ind$ID,ind$activity_ind, + # ind$activity_all,smband,plot.ori,plot.sm), ind_plot)) + #message(paste("Total time: ",round(difftime(Sys.time(),ptm,units="mins")[[1]],2)," mins",sep="")) + + ## trelliscope plot + #if(plot.tre==TRUE) { + # ind$activity_ind <- ind$activity_all <- NULL + # if(is.null(plot.tre.path)) { + # trelliscopejs::trelliscope(ind,name = "Individual Mean Activity Plot", nrow = 2, ncol = 2, + # path=getwd()) + # } else { + # trelliscopejs::trelliscope(ind,name = "Individual Mean Activity Plot", nrow = 2, ncol = 2, + # path=plot.tre.path) + # } + #} else { + # return(ind) + #} + + if(plot.tre==TRUE) { + print("plot.tre=TRUE but no plot is being generated due to the archived trelliscopejs package.") + } + return(ind) + + } else { + + ## auxillary information for plotting: y-axis activity max + act <- tibble::add_column(act,activity_max = ind_to_day(unlist(lapply(split(act$activity_sm,act$ID),function(x) + max(unlist(lapply(x,max))))),act)) + + ## generate plot: day observation + #ptm <- Sys.time() + message("No trelliscope individual plots are generated due to the archived trelliscopejs package.") + #act <- tibble::add_column(act,panel = trelliscopejs::pmap_plot(list(id=act$ID,id_Nday=act$ID_Nday, + # act_ori=act$activity,act_ind=act$activity_ind, + # act_all=act$activity_all,act_max=act$activity_max, + # band=smband,ori=plot.ori,lw=plot.sm), act_plot)) + #message(paste("Total time: ",round(difftime(Sys.time(),ptm,units="mins")[[1]],3)," mins",sep="")) + #check memory: format(object.size(act),units="Mb",standard="legacy") + + ## trelliscope plot + #if(plot.tre==TRUE) { + # act$activity <- act$activity_sm <- act$activity_ind <- act$activity_all <- NULL + # if(is.null(plot.tre.path)) { + # trelliscopejs::trelliscope(act,name = "Daily Activity Plot", nrow = 2, ncol = 2, + # path=getwd()) + # } else { + # trelliscopejs::trelliscope(act,name = "Day Activity Plot", nrow = 2, ncol = 2, + # path=plot.tre.path) + # } + #} else { + # return(act) + #} + + if(plot.tre==TRUE) { + print("plot.tre=TRUE but no plot is being generated due to the archived trelliscopejs package.") + } + return(act) + } +} diff --git a/R/var3.R b/R/var3.R index 567ecf6..773707a 100644 --- a/R/var3.R +++ b/R/var3.R @@ -1,17 +1,17 @@ -#' Covariate data for three individuals -#' -#' This is a synthetic data for 3 individuals, and the variables are "ID", age" and "gender". It is for -#' illustration purposes of the functions in the package \code{PML} only. -#' -#' @docType data -#' -#' @usage data(pa3) -#' -#' @format An object of class \code{data.frame}. -#' -#' @keywords datasets -#' -#' @examples -#' data(var3) -#' +#' Covariate data for three individuals +#' +#' This is a synthetic data for 3 individuals, and the variables are "ID", age" and "gender". It is for +#' illustration purposes of the functions in the package \code{PML} only. +#' +#' @docType data +#' +#' @usage data(pa3) +#' +#' @format An object of class \code{data.frame}. +#' +#' @keywords datasets +#' +#' @examples +#' data(var3) +#' #' @seealso \code{\link{tre}} \ No newline at end of file