Skip to content

Commit

Permalink
Merge pull request #2 from andjar/dev
Browse files Browse the repository at this point in the history 
Dev
  • Loading branch information
@andjar
andjar authored Mar 20, 2021
2 parents 5617a67 + ba69008 commit c55c3c3
Show file tree
Hide file tree
Showing 63 changed files with 1,331 additions and 543 deletions.
7 changes: 6 additions & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,18 +8,23 @@ Authors@R:
email = "anders.h.jarmund@ntnu.no",
comment = c(ORCID = "0000-0002-3923-1637"))
Description: Perform ASCA with various underlying models
URL: https://github.com/andjar/ALASCA
License: GPL-3
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.1.1
Depends:
data.table
Imports:
ggplot2,
ggpubr,
lme4,
lmerTest,
reshape2,
ggrepel
ggrepel,
Rfast,
data.table
Suggests:
knitr,
rmarkdown
Expand Down
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ S3method(summary,ALASCA)
export(ALASCA)
export(RMASCA)
export(SMASCA)
export(appendModel)
export(flipIt)
export(getCovars)
export(getLoadings)
Expand Down
6 changes: 6 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
# ALASCA 0.0.0.9

## Known bugs

* Reference level cannot be set when using multiple interactions
* `plotPart()` sometimes fails to provide more than a blank plot
295 changes: 195 additions & 100 deletions R/ALASCA.R
View file

Large diffs are not rendered by default.

163 changes: 47 additions & 116 deletions R/buildModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ buildModel <- function(object){
}
object <- runRegression(object)
if(object$method != "Rfast"){
# With Rfast, we've already got the coefficients
object <- getRegressionCoefficients(object)
}

Expand Down Expand Up @@ -44,117 +45,55 @@ buildModel <- function(object){
#' @return An ALASCA object
runRegression <- function(object){

regr.model <- list()
fdf <- data.frame()
cyts <- unique(object$df$variable)
cc <- 1
ccc <- 1

if(object$method == "Rfast"){
object$RegressionCoefficients <- Reduce(rbind,lapply(unique(object$df$variable), function(x){
start.time <- Sys.time()
df <- subset(object$df, variable == x)
ust <- unique(df$studyno)
for(i in 1:length(ust)){
df$studyno2[df$studyno == ust[i]] <- i
#start.time <- Sys.time()
df <- object$df[variable == x]
modmat <- model.matrix(object$newformula, data = df)
if(object$forceEqualBaseline){
modmat <- modmat[,!grepl(paste0("time",levels(object$df$time)[1]), colnames(modmat))]
}
modmat <- model.matrix(object$formula, data = df)
mod <- data.frame(
estimate = Rfast::rint.reg(y = df$estimate, x = as.matrix(modmat[,2:ncol(modmat)]), id = df$studyno2, ranef = FALSE)$be,
estimate = Rfast::rint.reg(y = df[,value],
x = modmat[,2:ncol(modmat)],
id = as.numeric(factor(df[,ID])),
ranef = FALSE)$be,
pvalue = 1,
covar = as.character(x),
variable = colnames(modmat)
)
end.time <- Sys.time()
cat("\n\n",end.time - start.time,"\n")
#end.time <- Sys.time()
#cat("\n\n",end.time - start.time,"\n")
mod
}))

return(object)
}else if(object$method == "LM"){
object$regr.model <- lapply(unique(object$df$variable), function(x){
modmat <- model.matrix(object$formula, data = object$df[variable == x])
modmat <- modmat[,-1]
if(object$forceEqualBaseline){
modmat <- modmat[,!grepl(paste0("time",unique(object$df$time)[1]), colnames(modmat))]
}
environment(object$newformula) <- environment()
regr.model <- lm(object$newformula, data = object$df[variable == x])
attr(regr.model, "name") <- x
regr.model
})
}else{
makeX <- function(object, regr.model){
X <- model.matrix(regr.model)
baselineLabel <- paste0("time", unique(object$df$time)[1])
if(object$doDebug){
cat(".... Baseline found to be: ",baselineLabel,"\n")
cat(".... Column names of X: ",colnames(X),"\n")
}
X <- X[, substr(colnames(X),1,nchar(baselineLabel)) != baselineLabel]
newFormula <- as.character(object$formula)
newFormulaPred <- strsplit(as.character(newFormula[3]), "\\+")[[1]]
newFormulaPred <- newFormulaPred[Reduce(cbind,lapply(newFormulaPred, function (x) grepl("\\(",x)))]
newFormula <- paste(newFormula[2],"~","X +",newFormulaPred, collapse = " ")
object$newFormula <- formula(newFormula)
if(object$doDebug){
cat(".... Old formula (given as ~ `outcome` `predictors`): ",as.character(object$formula),"\n")
cat(".... New formula (given as ~ `outcome` `predictors`): ",as.character(object$newFormula),"\n")
cat(".... Column names of X: ",colnames(X),"\n")
cat(".... Size of X (rowsxcols): ",nrow(X),"x",ncol(X),"\n")
}
object$X <- X
return(object)
}

if(object$forceEqualBaseline){
if(object$doDebug){
cat(".... Starts removing interaction between groups and first time point\n")
}
# Remove interaction between group and first time point
if(!object$missingMeasurements){
if(object$method == "LM"){
regr.model <- lm(object$formula, data = subset(object$df, variable == unique(object$df$variable)[1]))
}else{
regr.model <- lmerTest::lmer(object$formula, data = subset(object$df, variable == unique(object$df$variable)[1]))
}
object <- makeX(object, regr.model)
}
}

if(object$doDebug){
cat(".. Beginning regression\n")
}
object$regr.model <- lapply(unique(object$df$variable), function(i){
if(object$doDebug){
cat(".... Variable: ",i," (",levels(object$df$variable)[i],")\n")
cat(".... Number of rows in subset: ",nrow(subset(object$df, variable == i)),"\n")
}
object$regr.model <- lapply(unique(object$df$variable), function(x){
modmat <- model.matrix(object$formula, data = object$df[variable == x])
modmat <- modmat[,-1]
if(object$forceEqualBaseline){
if(!object$missingMeasurements){
if(object$method == "LM"){
regr.model <- lm(object$newFormula, data = subset(object$df, variable == i))
}else{
regr.model <- lmerTest::lmer(object$newFormula, data = subset(object$df, variable == i))
}
}else{

if(object$method == "LM"){
regr.model <- lm(object$formula, data = subset(object$df, variable == i))
}else{
regr.model <- lmerTest::lmer(object$formula, data = subset(object$df, variable == i))
}

object <- makeX(object, regr.model)

if(object$method == "LM"){
regr.model <- lm(object$newFormula, data = subset(object$df, variable == i))
}else{
regr.model <- lmerTest::lmer(object$newFormula, data = subset(object$df, variable == i))
}
}
}else{
if(object$method == "LM"){
regr.model <- lm(object$formula, data = subset(object$df, variable == i))
}else{
regr.model <- lmerTest::lmer(object$formula, data = subset(object$df, variable == i))
}
modmat <- modmat[,!grepl(paste0("time",levels(object$df$time)[1]), colnames(modmat), fixed = TRUE)]
}
attr(regr.model, "name") <- i
#modmat <- modmat[,ncol(modmat):1]
environment(object$newformula) <- environment()
regr.model <- lmerTest::lmer(object$newformula, data = object$df[variable == x])
attr(regr.model, "name") <- x
regr.model
}
)

names(object$regr.model) <- unique(object$df$variable)
})
}

names(object$regr.model) <- unique(object$df$variable)
return(object)
}

Expand All @@ -178,9 +117,7 @@ getRegressionCoefficients <- function(object){
rownames(a) <- NULL
a
}))
if(object$forceEqualBaseline){
fdf$variable[substr(fdf$variable,1,1) == "X"] <- substr(fdf$variable[substr(fdf$variable,1,1) == "X"],2,nchar(fdf$variable[substr(fdf$variable,1,1) == "X"]))
}
fdf$variable <- gsub("modmat","",fdf$variable, fixed = TRUE)

colnames(fdf) <- c("estimate", "pvalue", "covar", "variable")

Expand Down Expand Up @@ -220,7 +157,7 @@ separateLMECoefficients <- function(object){
object$RegressionCoefficients$comp <- "TIME"
if(object$separateTimeAndGroup){
object$RegressionCoefficients$comp[!(object$RegressionCoefficients$variable == "(Intercept)" |
(substr(object$RegressionCoefficients$variable, 1, 4) == "time" & !grepl(":",object$RegressionCoefficients$variable)))
(substr(object$RegressionCoefficients$variable, 1, 4) == "time" & !grepl(":",object$RegressionCoefficients$variable, fixed = "TRUE")))
] <- "GROUP"
}
return(object)
Expand All @@ -236,22 +173,11 @@ getEffectMatrix <- function(object){
if(!object$minimizeObject){
cat("Calculating effect matrix\n")
}
parts <- subset(object$df, variable == object$df$variable[1])
if(object$method == "LM"){
Dmatrix <- model.matrix(object$regr.model[[1]],"X")
}else if(object$method == "LMM"){
Dmatrix <- lme4::getME(object$regr.model[[1]],"X")
}else if(object$method == "Rfast"){

df <- subset(object$df, variable == object$df$variable[1])
parts <- object$df[variable == object$df$variable[1]]
#parts <- object$df[!duplicated(cbind(object$df$ID, object$df$time))]
Dmatrix <- model.matrix(object$formula, data = object$df[variable == object$df$variable[1]])
#Dmatrix <- Dmatrix[,ncol(Dmatrix):1]

Dmatrix <- model.matrix(object$formula, data = df)

}

if(object$forceEqualBaseline){
colnames(Dmatrix)[substr(colnames(Dmatrix),1,1) == "X"] <- substr(colnames(Dmatrix)[substr(colnames(Dmatrix),1,1) == "X"],2,nchar(colnames(Dmatrix)[substr(colnames(Dmatrix),1,1) == "X"]))
}
if(object$separateTimeAndGroup){
BmatrixTime <- object$RegressionCoefficients[object$RegressionCoefficients$comp == "TIME",c("covar","estimate","variable")]
BmatrixTime <- reshape2::dcast(BmatrixTime, formula = variable~covar, value.var = "estimate")
Expand Down Expand Up @@ -300,7 +226,12 @@ getEffectMatrix <- function(object){
}

object$parts$time <- parts$time
object$parts$group <- parts$group
if(object$keepTerms != ""){
keepTerms <- c("group", object$keepTerms)
object$parts$group <- apply( parts[ , ..keepTerms ] , 1 , paste , collapse = " - " )
}else{
object$parts$group <- parts$group
}
if(!object$minimizeObject){
cat("Finished calculating effect matrix!\n")
}
Expand Down
54 changes: 37 additions & 17 deletions R/plot.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,11 +178,14 @@ getLoadingPlot <- function(object, component = 1, effect = "time", decreasingLoa
if(object$validate){
g <- ggplot2::ggplot(loadings, ggplot2::aes(x = covars, y = loading, ymin = low, ymax = high)) + ggplot2::geom_pointrange(size = pointSize)
}else{
g <- ggplot2::ggplot(loadings, ggplot2::aes(x = covars, y = loading)) + ggplot2::geom_point()
if(any(colnames(loadings) == "model")){
g <- ggplot2::ggplot(loadings, ggplot2::aes(x = covars, y = loading, shape = model)) + ggplot2::geom_point()
}else{
g <- ggplot2::ggplot(loadings, ggplot2::aes(x = covars, y = loading)) + ggplot2::geom_point()
}
}

g <- g + myTheme +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, vjust = 1, hjust=1)) +
g <- g + myTheme + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, vjust = 1, hjust=1), legend.position = c(0.8, 0.8)) +
ggplot2::labs(x = "Variable",
y = paste0("PC",component, " (", round(100*ifelse(effect == "time", object$ALASCA$loading$explained$time[component],object$ALASCA$loading$explained$group[component]),2),"%)"))
if(!any(is.na(highlight))){
Expand Down Expand Up @@ -303,7 +306,7 @@ plotParts <- function(object,
valueColumn <- valueColumn
}
}else if(is(object, "ALASCA")){
df <- object$df
df <- object$dfRaw
valueColumn <- as.character(object$formula)[2]
if(any(participantColumn == FALSE)){
if(any(object$participantColumn == FALSE)){
Expand All @@ -325,7 +328,7 @@ plotParts <- function(object,
return(g)
}
if(any(is.na(variable))){
variabel <- unique(df$variable)
variable <- unique(df$variable)
}
g <- lapply(variable, function(xi){
plotFunction(df, timeColumn, valueColumn, participantColumn, xi, addSmooth, myTheme = myTheme)
Expand Down Expand Up @@ -417,7 +420,7 @@ plotVal <- function(object, component = 1, myTheme = ggplot2::theme_bw()){
)
})
)
dff$plotGroup <- paste0(dff$model,dff$group)
dff$plotGroup <- paste0(dff$model,"-",dff$group)
gsg <- ggplot2::ggplot(dff, ggplot2::aes(x = time, y = score, group = plotGroup, color = group)) +
ggplot2::geom_point(alpha = 0.2) + ggplot2::geom_line(alpha = 0.2) +
ggplot2::geom_point(data = subset(getScores(object)$group, PC == component), group = NA, alpha = 1, color = "black") +
Expand Down Expand Up @@ -460,7 +463,7 @@ plotVal <- function(object, component = 1, myTheme = ggplot2::theme_bw()){
)
})
)
dff$plotGroup <- paste0(dff$model,dff$group)
dff$plotGroup <- paste0(dff$model,"-",dff$group)
gs <- ggplot2::ggplot(dff, ggplot2::aes(x = time, y = score, group = plotGroup, color = group)) +
ggplot2::geom_point(alpha = 0.2) + ggplot2::geom_line(alpha = 0.2) +
ggplot2::geom_point(data = subset(getScores(object)$time, PC == component), group = NA, alpha = 1, color = "black") +
Expand Down Expand Up @@ -497,27 +500,43 @@ plotVal <- function(object, component = 1, myTheme = ggplot2::theme_bw()){
#' @return A ggplot2 objects\.
#'
#' @export
plotCovar <- function(object, covar = NA, tlab = NA, return_data = FALSE, myTheme = ggplot2::theme_bw()){
plotCovar <- function(object, covar = NA, tlab = NA, return_data = FALSE, myTheme = ggplot2::theme_bw(), pvalue = "shape"){
df <- getCovars(object)
if(any(is.na(covar))){
covar <- object$covars
covar <- unique(df$variable)
}
if(any(is.na(tlab))){
tlab <- covar
}
df <- getCovars(object)
for(i in 1:length(covar)){
df$tlab[df$variable == covar[i]] <- tlab[i]
}
df$pvalue_label <- ifelse(df$pvalue >= 0.05, "Not significant", ifelse(df$pvalue < 0.001, "< 0.001", ifelse(df$pvalue < 0.01, "< 0.01", "< 0.05")))
df$pvalue_sign <- ifelse(df$pvalue >= 0.05, "", ifelse(df$pvalue < 0.001, "***", ifelse(df$pvalue < 0.01, "**", "*")))
df$covar <- factor(df$covar, levels = unique(df$covar[order(df$estimate)]))
if(return_data){
return(df)
}else{
g <- ggplot2::ggplot(df, ggplot2::aes(x = estimate, y = covar, shape = pvalue_label)) +
ggplot2::geom_point() + ggplot2::geom_vline(xintercept = 0) +
ggplot2::scale_shape_manual(values = c("Not significant"=3, "< 0.05"=15, "< 0.01"=16, "< 0.001"=17, "Baseline"=5)) +
ggplot2::facet_wrap(~tlab) + ggplot2::labs(x = "Coefficient", y = "", shape = "P value") +
myTheme + ggplot2::theme(legend.position = "bottom", legend.box="vertical", legend.margin=ggplot2::margin())
if(pvalue == "shape"){
g <- ggplot2::ggplot(df, ggplot2::aes(x = estimate, y = covar, shape = pvalue_label)) +
ggplot2::geom_point() + ggplot2::geom_vline(xintercept = 0) +
ggplot2::scale_shape_manual(values = c("Not significant"=3, "< 0.05"=15, "< 0.01"=16, "< 0.001"=17, "Baseline"=5)) +
ggplot2::facet_wrap(~tlab) + ggplot2::labs(x = "Coefficient", y = "", shape = "P value") +
myTheme + ggplot2::theme(legend.position = "bottom", legend.box="vertical", legend.margin=ggplot2::margin())
}else if(pvalue == "star" | pvalue == "asterisk" | pvalue == "stars"){
g <- ggplot2::ggplot(df, ggplot2::aes(x = estimate, y = covar, label = pvalue_sign)) +
ggplot2::geom_point() + ggplot2::geom_vline(xintercept = 0) +
ggplot2::geom_text(hjust = 0.5, vjust = 0) +
ggplot2::facet_wrap(~tlab) + ggplot2::labs(x = "Coefficient", y = "", shape = "P value") +
myTheme
}else{
g <- ggplot2::ggplot(df, ggplot2::aes(x = estimate, y = covar)) +
ggplot2::geom_point() + ggplot2::geom_vline(xintercept = 0) +
ggplot2::facet_wrap(~tlab) + ggplot2::labs(x = "Coefficient", y = "", shape = "P value") +
myTheme + ggplot2::theme(legend.position = "bottom", legend.box="vertical", legend.margin=ggplot2::margin())
}


return(g)
}
}
Expand All @@ -536,7 +555,7 @@ plotCovar <- function(object, covar = NA, tlab = NA, return_data = FALSE, myThem
#' @export
plotProjection <- function(object, comp = c(1,2), return_data = FALSE, myTheme = ggplot2::theme_bw()){
df <- object$df
df$ID <- df[, colnames(df) == object$participantColumn]
df$ID <- df[, ID]
loadings_Time <- subset(getLoadings(object)$time, PC %in% comp)
loadings_Time <- reshape2::dcast(data = loadings_Time, covars ~ paste0("PC",PC), value.var = "loading")
df_time <- merge(df, loadings_Time, by.x = "variable", by.y = "covars")
Expand Down Expand Up @@ -587,7 +606,8 @@ plotProjection <- function(object, comp = c(1,2), return_data = FALSE, myTheme =
}))
df_time <- df_time[!duplicated(df_time),]
colnames(df_time) <- c("part", paste0("PC", comp[1]), paste0("PC", comp[2]), "time", "group")
g <- ggplot2::ggplot(df_time, ggplot2::aes_string(x = paste0("PC", comp[1]), y = paste0("PC", comp[2]), group = "part", color = "group")) + ggplot2::geom_point() + ggplot2::geom_line(alpha = 0.7, arrow = ggplot2::arrow(type="closed", length = ggplot2::unit(0.20,"cm"))) + myTheme
g <- ggplot2::ggplot(df_time, ggplot2::aes_string(x = paste0("PC", comp[1]), y = paste0("PC", comp[2]), group = "part", color = "group")) +
ggplot2::geom_point() + ggplot2::geom_line(alpha = 0.7, arrow = ggplot2::arrow(type="closed", length = ggplot2::unit(0.20,"cm"))) + myTheme
if(return_data){
return(df_time)
}else{
Expand Down
Loading

0 comments on commit c55c3c3

Please sign in to comment.