lmer_results$HIV_p[r] <- stat_sum_r$coefficients[3,5]
lmer_results$MSM_p[r] <- stat_sum_r$coefficients[4,5]
lmer_results$Tm_p[r] <- stat_sum_r$coefficients[5,5]
if(as.numeric(lmer_results$Diet_p[r])<0.05) {
if(as.numeric(lmer_results$Diet_Coef[r])>0) {
lmer_results$Diet_dir[r] <- "Up"
} else {
lmer_results$Diet_dir[r] <- "Down"
}
} else {
if(as.numeric(lmer_results$Diet_Coef[r])>0) {
lmer_results$Diet_dir[r] <- "Up (Insignificant)"
} else {
lmer_results$Diet_dir[r] <- "Down (Insignificant)"
}
}
if(as.numeric(lmer_results$HIV_p[r])<0.05) {
if(as.numeric(lmer_results$HIV_Coef[r])>0) {
lmer_results$HIV_dir[r] <- "Up"
} else {
lmer_results$HIV_dir[r] <- "Down"
}
} else {
if(as.numeric(lmer_results$HIV_Coef[r])>0) {
lmer_results$HIV_dir[r] <- "Up (Insignificant)"
} else {
lmer_results$HIV_dir[r] <- "Down (Insignificant)"
}
}
if(as.numeric(lmer_results$MSM_p[r])<0.05) {
if(as.numeric(lmer_results$MSM_Coef[r])>0) {
lmer_results$MSM_dir[r] <- "Up"
} else {
lmer_results$MSM_dir[r] <- "Down"
}
} else {
if(as.numeric(lmer_results$MSM_Coef[r])>0) {
lmer_results$MSM_dir[r] <- "Up (Insignificant)"
} else {
lmer_results$MSM_dir[r] <- "Down (Insignificant)"
}
}
if(as.numeric(lmer_results$Tm_p[r])<0.05) {
if(as.numeric(lmer_results$Tm_Coef[r])>0) {
lmer_results$Tm_dir[r] <- "Up"
} else {
lmer_results$Tm_dir[r] <- "Down"
}
} else {
if(as.numeric(lmer_results$Tm_Coef[r])>0) {
lmer_results$Tm_dir[r] <- "Up (Insignificant)"
} else {
lmer_results$Tm_dir[r] <- "Down (Insignificant)"
}
}
}
vars_lmer <- c("Diet_dir","HIV_dir","MSM_dir","Tm_dir")
lmer_results_long <- reshape(lmer_results, varying=vars_lmer, direction="long", v.names="Direction", times=vars_lmer)
lmer_results_long$time <- gsub("_dir", "", lmer_results_long$time)
lmer_results_long$var_names <- gsub("_Norm", "", lmer_results_long$var_names)
cd4_lmer_results <- lmer_results_long[which(substr(lmer_results_long$var_names,1,4)=="CD4+"),]
cd4_lm_plot <- ggplot(cd4_lmer_results, aes(x = time, y = var_names, fill=Direction)) +
scale_fill_manual(values=c("darkblue","lightblue","darkred","darksalmon"))+
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 1) +
labs(title = "Linear Mixed Effects Models", x = "Clinical Variables", y = "Normalized Measure")+
theme_classic() +
theme(legend.position = "right",
legend.title = element_text(size = 15),
legend.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.title.y = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.text.y = element_text(size = 15),
strip.text.x = element_text(size=15),
strip.text.y = element_text(size=15),
ggh4x.facet.nestline = element_blank(),
strip.background.x = element_rect(
color="black", fill="gray", linetype="solid"
),
strip.background.y = element_rect(
color="black", fill="gray", linetype="solid"
))
cd8_lmer_results <- lmer_results_long[which(substr(lmer_results_long$var_names,1,4)=="CD8+"),]
cd8_lm_plot <- ggplot(cd8_lmer_results, aes(x = time, y = var_names, fill=Direction)) +
scale_fill_manual(values=c("darkblue","lightblue","darkred","darksalmon"))+
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 1) +
labs(title = "Linear Mixed Effects Models", x = "Clinical Variables", y = "Normalized Measure")+
theme_classic() +
theme(legend.position = "right",
legend.title = element_text(size = 15),
legend.text = element_text(size = 15),
axis.text.x = element_text(size = 15),
axis.title.y = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.text.y = element_text(size = 15),
strip.text.x = element_text(size=15),
strip.text.y = element_text(size=15),
ggh4x.facet.nestline = element_blank(),
strip.background.x = element_rect(
color="black", fill="gray", linetype="solid"
),
strip.background.y = element_rect(
color="black", fill="gray", linetype="solid"
))
out_vars <- c("HDLNorm","LDLNorm","TriglyceridesNorm","HOMAIRNorm","IL6Norm","CRPNorm","LBPNorm")
lm_out_data <- unique(expand.grid(un_gps,sig_ag_proc_cells,out_vars))
colnames(lm_out_data) <- c("Group","Cell","Outcome")
lm_out_data["Coef"] <- NA
lm_out_data["pval"] <- NA
lm_out_data["Sig"] <- NA
for (k in 1:nrow(lm_out_data)) {
gp_k <- as.character(lm_out_data$Group[k])
cl_k <- as.character(lm_out_data$Cell[k])
out_k <- as.character(lm_out_data$Outcome[k])
data_k <- bd_all_proc_imp[which(bd_all_proc_imp$Group==gp_k),which(colnames(bd_all_proc_imp) %in% c("StudyID", out_k, cl_k))]
colnames(data_k)[which(colnames(data_k)==out_k)] <- "Outcome"
colnames(data_k)[which(colnames(data_k)==cl_k)] <- "Cell"
data_k$Cell <- as.numeric(data_k$Cell)
data_k$Outcome <- as.numeric(data_k$Outcome)
lmer_k <- lmer(Outcome ~ Cell + (1 | StudyID), data = data_k)
stat_mixed_k <- mixed(lmer_k, data=data_k)
stat_sum_k <- summary(stat_mixed_k)
lm_out_data$Coef[k] <- stat_sum_k$coefficients[2,1]
lm_out_data$pval[k] <- stat_sum_k$coefficients[2,5]
if (lm_out_data$Coef[k]>0) {
sig_k <- "Positive"
} else {
sig_k <- "Negative"
}
if (lm_out_data$pval[k]<0.05) {
sig_k <- paste0(sig_k, " (Significant)")
}
lm_out_data$Sig[k] <- sig_k
}
lm_out_data$Cell <- gsub("_Norm","", lm_out_data$Cell)
lm_out_data$Cell <- gsub("(% of CD3- CD45+)","", lm_out_data$Cell)
lm_out_data$Outcome <- gsub("Norm","",lm_out_data$Outcome)
lm_out_plot <- ggplot(lm_out_data, aes(x = Group, y = Cell, fill=Sig)) +
scale_fill_manual(values=c("lightblue","darkblue","darksalmon","darkred"))+
geom_dotplot(binaxis = "y", stackdir = "center", dotsize = 1) +
labs(title = "Linear Mixed Effects Models", x = "Outcome", y = "Cell Population")+
theme_classic() +
theme(legend.position = "right",
legend.title = element_text(size = 15),
legend.text = element_text(size = 15),
axis.text.x = element_text(size = 12, angle=-45, hjust=0),
axis.title.y = element_text(size = 15),
axis.title.x = element_text(size = 15),
axis.text.y = element_text(size = 15),
strip.text.x = element_text(size=15),
strip.text.y = element_text(size=15),
ggh4x.facet.nestline = element_blank(),
strip.background.x = element_rect(
color="black", fill="gray", linetype="solid"
),
strip.background.y = element_rect(
color="black", fill="gray", linetype="solid"
))+
facet_nested(cols=vars(Outcome),
scales="free_y")
ann_data <- bd_all_proc_imp[,1:17]
pcoa_data <- bd_all_proc_imp[,18:ncol(bd_all_proc_imp)]
#pcoa_data_raw <- bd_all_raw[,18:ncol(bd_all_proc_imp)]
for (j in 1:nrow(pcoa_data_raw)) {
for (k in 1:ncol(pcoa_data_raw)) {
if (is.na(pcoa_data_raw[j,k])) {
pcoa_data_raw[j,k] <- 0
}
}
}
data=meta_data_bd
meta_data_bd$HDLNorm <- as.numeric(meta_data_bd$HDLNorm)
meta_data_bd$LDLNorm <- as.numeric(meta_data_bd$LDLNorm)
sub_vars=c("Diet","HIV","MSM")
dep_vars=c("HDLNorm","LDLNorm")
gp_var="Group"
normality=TRUE
adjustment="none"
#First we pull the length of the subset variables
sub_var_num <- length(sub_vars)
#We generate a blank list that we subsequently add to
sub_var_list <- list()
#We go through each subset variables and add to create a list of variables with select values
for (j in 1:sub_var_num) {
#We pull the individual subset variables
sub_name_j <- sub_vars[j]
#We then pull their values
sub_vals_j <- as.vector(unique(data[,which(colnames(data)==sub_name_j)]))
#We then add them to the list
sub_var_list[[sub_name_j]] <- sub_vals_j
}
#Now we need to pull the group values to generate pairwise comparisons
gp_vals <- as.vector(unique(data[,which(colnames(data)==gp_var)]))
#We generate pairwise combinations of 2
tm_combs <- combn(gp_vals, 2, simplify = FALSE)
#We collapse them into one string
tm_strings <- sapply(tm_combs, function(x) paste(x, collapse = "-xxx-"))
#We now add them to the sub_var_list to allow us to generate a table
sub_var_list[["combo"]] <- tm_strings
#Now we generate a table with all the subset variables those values
subset_table <- expand.grid(sub_var_list)
#Now we broaden to include that with every dependent variable
#First we pull the number of dep_vars
dep_var_num <- length(dep_vars)
#Now we generate a blank table that we add to with each dependent variable
stats_sum_table <- c()
#We go through each dependent variable to add it to the stats_sum_table
for (d in 1:dep_var_num) {
#For each variable we adjoin a column in the front with the dependent variable name
stats_sum_table_d <- cbind(rep(dep_vars[d],nrow(subset_table)),subset_table)
#We rename that first adjoining column Measure
colnames(stats_sum_table_d)[1] <- "Measure"
#We now adjoin it to the stats_sum_table
stats_sum_table <- rbind(stats_sum_table,stats_sum_table_d)
}
#We now add new fields for additional specified fields
#First we split the combination into 2 groups
stats_sum_table <- separate(stats_sum_table, combo, into = c("group1", "group2"), sep = "-xxx-")
#Now we add blank fields for the stats
#One for the p-value
stats_sum_table["p-value"] <- NA
#One for the adjusted p-value
stats_sum_table["padj"] <- NA
#Now we go through ever combination of subsets of data
View(stats_sum_table)
subset_table
sub_var_list
test <- as.data.frame(sub_var_list)
test <- as.data.frame(expand.grid(sub_var_list))
View(test)
subs <- unique(stats_sum_table[,1:(ncol(stats_sum_table)-4)])
View(subs)
test <- "which(colnames(data)=='HDLNorm')"
get(test)
r=1
#We pull the measure
meas_r <- as.character(subs$Measure[r])
#We pull the measure columns from the data set
meas_col <- which(colnames(data)==meas_r)
#Now we subset the data by repeatedly filtering based on the subset filters
#We start with the full data set
data_r <- data
#We then go through the columns and subsequently filter ou the values
for (c in 2:ncol(subs)) {
#We pull the column name
colname_c <- as.character(colnames(subs)[c])
#We then identify that column in the data set
col_c <- which(colnames(data) <- colname_c)
#We pull the column value
sub_val_c <- as.character(subs[r,c])
data_r <- data_r[which(data_r[,col_c]==sub_val_c)]
}
c
colnames(c)
colname_c
col_c
data=meta_data_bd
meta_data_bd$HDLNorm <- as.numeric(meta_data_bd$HDLNorm)
meta_data_bd$LDLNorm <- as.numeric(meta_data_bd$LDLNorm)
sub_vars=c("Diet","HIV","MSM")
dep_vars=c("HDLNorm","LDLNorm")
gp_var="Group"
normality=TRUE
adjustment="none"
#First we pull the length of the subset variables
sub_var_num <- length(sub_vars)
#We generate a blank list that we subsequently add to
sub_var_list <- list()
#We go through each subset variables and add to create a list of variables with select values
for (j in 1:sub_var_num) {
#We pull the individual subset variables
sub_name_j <- sub_vars[j]
#We then pull their values
sub_vals_j <- as.vector(unique(data[,which(colnames(data)==sub_name_j)]))
#We then add them to the list
sub_var_list[[sub_name_j]] <- sub_vals_j
}
#Now we need to pull the group values to generate pairwise comparisons
gp_vals <- as.vector(unique(data[,which(colnames(data)==gp_var)]))
#We generate pairwise combinations of 2
tm_combs <- combn(gp_vals, 2, simplify = FALSE)
#We collapse them into one string
tm_strings <- sapply(tm_combs, function(x) paste(x, collapse = "-xxx-"))
#We now add them to the sub_var_list to allow us to generate a table
sub_var_list[["combo"]] <- tm_strings
#Now we generate a table with all the subset variables those values
subset_table <- expand.grid(sub_var_list)
#Now we broaden to include that with every dependent variable
#First we pull the number of dep_vars
dep_var_num <- length(dep_vars)
#Now we generate a blank table that we add to with each dependent variable
stats_sum_table <- c()
#We go through each dependent variable to add it to the stats_sum_table
for (d in 1:dep_var_num) {
#For each variable we adjoin a column in the front with the dependent variable name
stats_sum_table_d <- cbind(rep(dep_vars[d],nrow(subset_table)),subset_table)
#We rename that first adjoining column Measure
colnames(stats_sum_table_d)[1] <- "Measure"
#We now adjoin it to the stats_sum_table
stats_sum_table <- rbind(stats_sum_table,stats_sum_table_d)
}
#We now add new fields for additional specified fields
#First we split the combination into 2 groups
stats_sum_table <- separate(stats_sum_table, combo, into = c("group1", "group2"), sep = "-xxx-")
#Now we add blank fields for the stats
#One for the p-value
stats_sum_table["p-value"] <- NA
#One for the adjusted p-value
stats_sum_table["padj"] <- NA
#Now we go through ever combination of subsets of data and measure and perform the statistical comparisons
#Because the 4th to last columns though the last column specify the pariwise groups and the p values we
subs <- unique(stats_sum_table[,1:(ncol(stats_sum_table)-4)])
#Now we go thorugh each measure (first column in subs) and data subset filters (the remaining columns) and pull that
#data do statistics and add statistics to the stats table
r=1
#We pull the measure
meas_r <- as.character(subs$Measure[r])
#We pull the measure columns from the data set
meas_col <- which(colnames(data)==meas_r)
#Now we subset the data by repeatedly filtering based on the subset filters
#We start with the full data set
data_r <- data
#We then go through the columns and subsequently filter ou the values
for (c in 2:ncol(subs)) {
#We pull the column name
colname_c <- as.character(colnames(subs)[c])
#We then identify that column in the data set
col_c <- which(colnames(data) = colname_c)
#We pull the measure
meas_r <- as.character(subs$Measure[r])
#We pull the measure columns from the data set
meas_col <- which(colnames(data)==meas_r)
#Now we subset the data by repeatedly filtering based on the subset filters
#We start with the full data set
data_r <- data
#We then go through the columns and subsequently filter ou the values
for (c in 2:ncol(subs)) {
#We pull the column name
colname_c <- as.character(colnames(subs)[c])
#We then identify that column in the data set
col_c <- which(colnames(data) == colname_c)
#We pull the column value
sub_val_c <- as.character(subs[r,c])
data_r <- data_r[which(data_r[,col_c]==sub_val_c)]
}
c
colname_c
col_c
sub_val_c <- as.character(subs[r,c])
sub_val_c
data_r <- data_r[which(data_r[,col_c]==sub_val_c)]
#We pull the column name
colname_c <- as.character(colnames(subs)[c])
#We then identify that column in the data set
col_c <- which(colnames(data) == colname_c)
colnames(data_r)[col_c] <- "Filter"
#We pull the column value
sub_val_c <- as.character(subs[r,c])
data=meta_data_bd
meta_data_bd$HDLNorm <- as.numeric(meta_data_bd$HDLNorm)
meta_data_bd$LDLNorm <- as.numeric(meta_data_bd$LDLNorm)
sub_vars=c("Diet","HIV","MSM")
dep_vars=c("HDLNorm","LDLNorm")
gp_var="Group"
normality=TRUE
adjustment="none"
#First we pull the length of the subset variables
sub_var_num <- length(sub_vars)
#We generate a blank list that we subsequently add to
sub_var_list <- list()
#We go through each subset variables and add to create a list of variables with select values
for (j in 1:sub_var_num) {
#We pull the individual subset variables
sub_name_j <- sub_vars[j]
#We then pull their values
sub_vals_j <- as.vector(unique(data[,which(colnames(data)==sub_name_j)]))
#We then add them to the list
sub_var_list[[sub_name_j]] <- sub_vals_j
}
#Now we need to pull the group values to generate pairwise comparisons
gp_vals <- as.vector(unique(data[,which(colnames(data)==gp_var)]))
#We generate pairwise combinations of 2
tm_combs <- combn(gp_vals, 2, simplify = FALSE)
#We collapse them into one string
tm_strings <- sapply(tm_combs, function(x) paste(x, collapse = "-xxx-"))
#We now add them to the sub_var_list to allow us to generate a table
sub_var_list[["combo"]] <- tm_strings
#Now we generate a table with all the subset variables those values
subset_table <- expand.grid(sub_var_list)
#Now we broaden to include that with every dependent variable
#First we pull the number of dep_vars
dep_var_num <- length(dep_vars)
#Now we generate a blank table that we add to with each dependent variable
stats_sum_table <- c()
#We go through each dependent variable to add it to the stats_sum_table
for (d in 1:dep_var_num) {
#For each variable we adjoin a column in the front with the dependent variable name
stats_sum_table_d <- cbind(rep(dep_vars[d],nrow(subset_table)),subset_table)
#We rename that first adjoining column Measure
colnames(stats_sum_table_d)[1] <- "Measure"
#We now adjoin it to the stats_sum_table
stats_sum_table <- rbind(stats_sum_table,stats_sum_table_d)
}
#We now add new fields for additional specified fields
#First we split the combination into 2 groups
stats_sum_table <- separate(stats_sum_table, combo, into = c("group1", "group2"), sep = "-xxx-")
#Now we add blank fields for the stats
#One for the p-value
stats_sum_table["p-value"] <- NA
#One for the adjusted p-value
stats_sum_table["padj"] <- NA
#Now we go through ever combination of subsets of data and measure and perform the statistical comparisons
#Because the 4th to last columns though the last column specify the pariwise groups and the p values we
subs <- unique(stats_sum_table[,1:(ncol(stats_sum_table)-4)])
#Now we go thorugh each measure (first column in subs) and data subset filters (the remaining columns) and pull that
#data do statistics and add statistics to the stats table
for (r in 1:nrow(subs)) {
#We pull the measure
meas_r <- as.character(subs$Measure[r])
#We pull the measure columns from the data set
meas_col <- which(colnames(data)==meas_r)
#Now we subset the data by repeatedly filtering based on the subset filters
#We start with the full data set
data_r <- data
#We then go through the columns and subsequently filter ou the values
for (c in 2:ncol(subs)) {
#We pull the column name
colname_c <- as.character(colnames(subs)[c])
#We then identify that column in the data set
col_c <- which(colnames(data) == colname_c)
colnames(data_r)[col_c] <- "Filter"
#We pull the column value
sub_val_c <- as.character(subs[r,c])
data_r <- data_r[which(data_r$Filter==sub_val_c)]
colnames(data_r)[col_c] <- colname_c
}
#Now we focus on the columns of the measures and the groups from our stastic
data_r_filt <- data_r[,which(colnames(data_r) %in% c(meas_r, gp_var))]
)
data_r_filt <- data_r[,which(colnames(data_r) %in% c(meas_r, gp_var))]
data_r_filt
data_r_filt
#We pull the measure
meas_r <- as.character(subs$Measure[r])
#We pull the measure columns from the data set
meas_col <- which(colnames(data)==meas_r)
#Now we subset the data by repeatedly filtering based on the subset filters
#We start with the full data set
data_r <- data
#We then go through the columns and subsequently filter ou the values
for (c in 2:ncol(subs)) {
#We pull the column name
colname_c <- as.character(colnames(subs)[c])
#We then identify that column in the data set
col_c <- which(colnames(data) == colname_c)
colnames(data_r)[col_c] <- "Filter"
#We pull the column value
sub_val_c <- as.character(subs[r,c])
data_r <- data_r[which(data_r$Filter==sub_val_c),]
colnames(data_r)[col_c] <- colname_c
}
#Now we focus on the columns of the measures and the groups from our stastic
data_r_filt <- data_r[,which(colnames(data_r) %in% c(meas_r, gp_var))]
View(data_r_filt)
nrow(data_r_filt)
test <- data[which(data$Diet=="Western" & data$Group=="HIV(+)MSM"),]
nrow(test)
#Now we focus on the columns of the measures and the groups for the statistics
data_r_filt <- data_r[,which(colnames(data_r) %in% c(meas_r, gp_var))]
colnames(data_r_filt)[which(colnames(data_r_filt)==meas_r)] <- "Amount"
colnames(group_r_filt)[which(colnames(data_r_filt)==meas_r)] <- "Group"
#Now we focus on the columns of the measures and the groups for the statistics
data_r_filt <- data_r[,which(colnames(data_r) %in% c(meas_r, gp_var))]
colnames(data_r_filt)[which(colnames(data_r_filt)==meas_r)] <- "Amount"
colnames(data_r_filt)[which(colnames(data_r_filt)==meas_r)] <- "Group"
stat_r <- anova_test(data=data_r_filt, Amount ~ Group)
data_r_filt <- data
#We change the column names to apply to all measures and group names for the stats
colnames(data_r_filt)[which(colnames(data_r_filt)==meas_r)] <- "Amount"
colnames(data_r_filt)[which(colnames(data_r_filt)==meas_r)] <- "Group"
#We ensure that the amount column is numeric for statistics
data_r_filt$Amount <- as.numeric(data_r_filt$Amount)
#If the data is normal we do an ANOVA
stat_r <- anova_test(data=data_r_filt, Amount ~ Group)
version(pkg = "tidyr")
?anova_test
version(pkg = "rstatix")
?p.adjust
?aov
?TukeyHSD
model_r <- aov(data=data_r_filt, Amount ~ Group)
model_r
model_r$coefficients
sum_r <- summary(model_r)
sum_r
sum_r[1,5]
model_r$coefficients
model_r$terms
#The ANOVA model is put in
model_r <- aov(data=data_r_filt, Amount ~ Group)
#We generate a summary
sum_r <- as.data.frame(summary(model_r))
model_r <- aov(data=data_r_filt, Amount ~ Group)
#We generate a summary
sum_r <- summary(model_r)
sum_r
coef(sum_r)
coef(model-r)
coef(model_r)
sum_r
summary(model_r)[[1]][["Pr(>F)"]][1]
?TukeyHSD
