如何汇总统计2WAPost HOC？
How to summary statistics 2WA Post HOC?
structure(list(AGI = c(“ATCG01240”, “ATCG01310”, “ATMG00070”), aox2_0h__1 = c(15.79105291, 14.82652303, 14.70630068), aox2_0h__2 = c(16.06494674, 14.50610036, 14.52189807), aox2_0h__3 = c(14.64596287, 14.73266459, 13.07143141), aox2_0h__4 = c(15.71713641, 15.15430026, 16.32190068 ), aox2_12h__1 = c(14.99030606, 15.08046949, 15.8317372), aox2_12h__2 = c(15.15569857, 14.98996474, 14.64862254), aox2_12h__3 = c(15.12144791, 14.90111092, 14.59618842), aox2_12h__4 = c(14.25648197, 15.09832061, 14.64442686), aox2_24h__1 = c(15.23997241, 14.80968391, 14.22573239 ), aox2_24h__2 = c(15.57551513, 14.94861669, 15.18808897), aox2_24h__3 = c(15.04928714, 14.83758685, 13.06948037), aox2_24h__4 = c(14.79035385, 14.93873234, 14.70402827), aox5_0h__1 = c(15.8245918, 14.9351844, 14.67678306), aox5_0h__2 = c(15.75108628, 14.85867002, 14.45704948 ), aox5_0h__3 = c(14.36545859, 14.79296855, 14.82177912), aox5_0h__4 = c(14.80626019, 13.43330964, 16.33482718), aox5_12h__1 = c(14.66327372, 15.22571466, 16.17761867), aox5_12h__2 = c(14.58089039, 14.98545497, 14.4331578), aox5_12h__3 = c(14.58091828, 14.86139511, 15.83898617 ), aox5_12h__4 = c(14.48097297, 15.1420725, 13.39369381), aox5_24h__1 = c(15.41855602, 14.9890092, 13.92629626), aox5_24h__2 = c(15.78386057, 15.19372889, 14.63254456), aox5_24h__3 = c(15.55321382, 14.82013321, 15.74324956), aox5_24h__4 = c(14.53085803, 15.12196994, 14.81028556 ), WT_0h__1 = c(14.0535031, 12.45484834, 14.89102226), WT_0h__2 = c(13.64720361, 15.07144643, 14.99836235), WT_0h__3 = c(14.28295759, 13.75283646, 14.98220861), WT_0h__4 = c(14.79637443, 15.1108037, 15.21711524 ), WT_12h__1 = c(15.05711898, 13.33689777, 14.81064042), WT_12h__2 = c(14.83846779, 13.62497318, 14.76356308), WT_12h__3 = c(14.77215863, 14.72814995, 13.0835214), WT_12h__4 = c(14.70685445, 14.98527337, 16.12727292), WT_24h__1 = c(15.43813077, 14.56918572, 14.92146565 ), WT_24h__2 = c(16.05986898, 14.70583866, 15.64566505), WT_24h__3 = c(14.87721853, 13.22461859, 16.34119942), WT_24h__4 = c(14.92822133, 14.74382383, 12.79146694)), class = “data.frame”, row.names = c(NA, -3L))
我必须总结每个基因（即“ATCG01240”、“ATCG01310”、“ATMG00070”）时间点的数据； 4 次重复的平均值、标准误差并对每个基因、每个时间点进行多重比较（Posthoc 2way ANOVA；即 WT-aox2、WT-aox5、aox2-aox5）。
试一试。由于有关数据帧是列表的错误消息，我最终无法获得方差分析的一种方式。看看您可以用它做什么 - 或者其他人可以提供帮助。
df <- structure(list(AGI = c("ATCG01240", "ATCG01310", "ATMG00070"), aox2_0h__1 = c(15.79105291, 14.82652303, 14.70630068), aox2_0h__2 = c(16.06494674, 14.50610036, 14.52189807), aox2_0h__3 = c(14.64596287, 14.73266459, 13.07143141), aox2_0h__4 = c(15.71713641, 15.15430026, 16.32190068 ), aox2_12h__1 = c(14.99030606, 15.08046949, 15.8317372), aox2_12h__2 = c(15.15569857, 14.98996474, 14.64862254), aox2_12h__3 = c(15.12144791, 14.90111092, 14.59618842), aox2_12h__4 = c(14.25648197, 15.09832061, 14.64442686), aox2_24h__1 = c(15.23997241, 14.80968391, 14.22573239 ), aox2_24h__2 = c(15.57551513, 14.94861669, 15.18808897), aox2_24h__3 = c(15.04928714, 14.83758685, 13.06948037), aox2_24h__4 = c(14.79035385, 14.93873234, 14.70402827), aox5_0h__1 = c(15.8245918, 14.9351844, 14.67678306), aox5_0h__2 = c(15.75108628, 14.85867002, 14.45704948 ), aox5_0h__3 = c(14.36545859, 14.79296855, 14.82177912), aox5_0h__4 = c(14.80626019, 13.43330964, 16.33482718), aox5_12h__1 = c(14.66327372, 15.22571466, 16.17761867), aox5_12h__2 = c(14.58089039, 14.98545497, 14.4331578), aox5_12h__3 = c(14.58091828, 14.86139511, 15.83898617 ), aox5_12h__4 = c(14.48097297, 15.1420725, 13.39369381), aox5_24h__1 = c(15.41855602, 14.9890092, 13.92629626), aox5_24h__2 = c(15.78386057, 15.19372889, 14.63254456), aox5_24h__3 = c(15.55321382, 14.82013321, 15.74324956), aox5_24h__4 = c(14.53085803, 15.12196994, 14.81028556 ), WT_0h__1 = c(14.0535031, 12.45484834, 14.89102226), WT_0h__2 = c(13.64720361, 15.07144643, 14.99836235), WT_0h__3 = c(14.28295759, 13.75283646, 14.98220861), WT_0h__4 = c(14.79637443, 15.1108037, 15.21711524 ), WT_12h__1 = c(15.05711898, 13.33689777, 14.81064042), WT_12h__2 = c(14.83846779, 13.62497318, 14.76356308), WT_12h__3 = c(14.77215863, 14.72814995, 13.0835214), WT_12h__4 = c(14.70685445, 14.98527337, 16.12727292), WT_24h__1 = c(15.43813077, 14.56918572, 14.92146565 ), WT_24h__2 = c(16.05986898, 14.70583866, 15.64566505), WT_24h__3 = c(14.87721853, 13.22461859, 16.34119942), WT_24h__4 = c(14.92822133, 14.74382383, 12.79146694)), class = "data.frame", row.names = c(NA, -3L))
#https://www.r-statistics.com/tag/transpose/
mdf<-melt(df)
names(mdf)<-c("gene","TRT_time","value") #name the new columns
#need to parse TRT_time into separate TRT, time, and sample number columns.
#
foo <- data.frame(do.call('rbind', strsplit(as.character(mdf$TRT_time),'_',fixed=TRUE)))
mdf[2] <- NULL
typeof(foo)
typeof(mdf)
#need to merge these two lists
mdf2<-cbind(foo, mdf)
names(mdf2)<-c("treatment","time","blank", "sampleID","gene","value") #name the new columns
#remove the blank column
mdf2[3] <- NULL
agg.mean <- aggregate(value ~ gene + treatment + time,data=mdf2,mean)
names(agg.mean) <- c("gene","treatment","time", "mean") #name the new columns
#create a function for standard error of the mean
st.err <- function(x) {
sd(x)/sqrt(length(x))
}
st.err <- function(x, na.rm=FALSE) {
if(na.rm==TRUE) x <- na.omit(x)
sd(x)/sqrt(length(x))
}
agg.SEM <- aggregate(value ~ gene + treatment + time,data=mdf2,st.err)
names(agg.SEM) <- c("gene","treatment","time", "SEM") #name the new columns
agg.meanSEM <- cbind (agg.mean, agg.SEM)#prove they match by row
agg.meanSEM <- cbind (agg.mean, agg.SEM$SEM)
#sort the data in mdf2 in preparation for 2 way ANOVA repeated measures
names(mdf2)
mdf3<-mdf2[with(mdf2, order(-sampleID, gene, treatment, time)), ]
#https://rcompanion.org/rcompanion/d_08.html
#Two-way Anova
# An R Companion for the Handbook of Biological Statistics Salvatore S. Mangiafico
install.packages("Rmisc")
library(Rmisc)
sum = summarySE(mdf2, measurevar="value", groupvars=c("gene", "treatment","time"))
sum
#In case you are interested.
boxplot(value ~ gene:treatment:time,
data = mdf2,
xlab = "Gene X Treatment X Time",
ylab = "Gene Expression Level")
#Fit the linear model and conduct ANOVA
#original from source above
# model = lm(value ~ treatment + time + treatment:time,
# data=mdf2)
# Anova(model, type="III")
#I could not get the type added into the split commands; so I do not know which type is being run.
install.packages("car")
library(car)
#using split to 'aggregate' by gene for 2 way Anova.
#
s <- split(mdf2, mdf2$gene)
sink("sink-examp.txt")#this will capture results in a text file to your working directory
do.call(rbind, lapply(s, function(x) Anova(lm(value ~ treatment + time + treatment:time, x))[ ]))
sink()
###############
#Now to arrange for a one way ANOVA to run the TukeyHSD.
df_posthoc<-t(df)
df_posthoc<-data.frame(df_posthoc)
typeof(df_posthoc)
#Convert first row to headers.
#
names(df_posthoc) <- as.matrix(df_posthoc[1, ])
df_posthoc <- df_posthoc[-1, ]
df_posthoc[] <- lapply(df_posthoc, function(x) type.convert(as.character(x)))
df_posthoc
names(df_posthoc)
#
df_posthoc <- cbind(rownames(df_posthoc), data.frame(AD, row.names=NULL))
colnames(df_posthoc)[1] <- "TRT_time" #name the first column
#need to parse TRT_time into separate TRT, time, and sample number columns.
#
foo2 <- data.frame(do.call('rbind', strsplit(as.character(AD$TRT_time),'_',fixed=TRUE)))
foo2[3] <- NULL
colnames(foo2)[1] <- "treatment" #name the first column
colnames(foo2)[2] <- "time" #name the first column
colnames(foo2)[3] <- "sampleID" #name the first column
df_posthoc[1] <- NULL
typeof(foo2)
typeof(df_posthoc)
#need to merge these two lists
df_posthoc<-cbind(foo2, df_posthoc)
typeof(df_posthoc)
AD<-data.frame(df_posthoc)
typeof(df_posthoc)
mode(df_posthoc)
data.matrix(df_posthoc)
typeof(df_posthoc)
mode(df_posthoc)
df_posthoc <- mapply(df_posthoc, FUN=as.numeric)
df_posthoc<-data.frame(df_posthoc)
#Useful site. https://www.statmethods.net/stats/anova.html
1:ncol(df_posthoc)
AVz<- rep(NA,ncol(df_posthoc))# Creates a table, AV, with the same number of columns as in file.
for (i in 2:ncol(df_posthoc)){
AVz<-summary(aov(df_posthoc~treatment, data=df_posthoc))
tk<-TukeyHSD((aov(df_posthoc~treatment, data=df_posthoc)))
capture.output(AVz,file="One-WayANOVASummarycolumn.doc")
capture.output(tk,file="TukeyHSDResultscolumn.doc")
}
structure(list(AGI = c(“ATCG01240”, “ATCG01310”, “ATMG00070”), aox2_0h__1 = c(15.79105291, 14.82652303, 14.70630068), aox2_0h__2 = c(16.06494674, 14.50610036, 14.52189807), aox2_0h__3 = c(14.64596287, 14.73266459, 13.07143141), aox2_0h__4 = c(15.71713641, 15.15430026, 16.32190068 ), aox2_12h__1 = c(14.99030606, 15.08046949, 15.8317372), aox2_12h__2 = c(15.15569857, 14.98996474, 14.64862254), aox2_12h__3 = c(15.12144791, 14.90111092, 14.59618842), aox2_12h__4 = c(14.25648197, 15.09832061, 14.64442686), aox2_24h__1 = c(15.23997241, 14.80968391, 14.22573239 ), aox2_24h__2 = c(15.57551513, 14.94861669, 15.18808897), aox2_24h__3 = c(15.04928714, 14.83758685, 13.06948037), aox2_24h__4 = c(14.79035385, 14.93873234, 14.70402827), aox5_0h__1 = c(15.8245918, 14.9351844, 14.67678306), aox5_0h__2 = c(15.75108628, 14.85867002, 14.45704948 ), aox5_0h__3 = c(14.36545859, 14.79296855, 14.82177912), aox5_0h__4 = c(14.80626019, 13.43330964, 16.33482718), aox5_12h__1 = c(14.66327372, 15.22571466, 16.17761867), aox5_12h__2 = c(14.58089039, 14.98545497, 14.4331578), aox5_12h__3 = c(14.58091828, 14.86139511, 15.83898617 ), aox5_12h__4 = c(14.48097297, 15.1420725, 13.39369381), aox5_24h__1 = c(15.41855602, 14.9890092, 13.92629626), aox5_24h__2 = c(15.78386057, 15.19372889, 14.63254456), aox5_24h__3 = c(15.55321382, 14.82013321, 15.74324956), aox5_24h__4 = c(14.53085803, 15.12196994, 14.81028556 ), WT_0h__1 = c(14.0535031, 12.45484834, 14.89102226), WT_0h__2 = c(13.64720361, 15.07144643, 14.99836235), WT_0h__3 = c(14.28295759, 13.75283646, 14.98220861), WT_0h__4 = c(14.79637443, 15.1108037, 15.21711524 ), WT_12h__1 = c(15.05711898, 13.33689777, 14.81064042), WT_12h__2 = c(14.83846779, 13.62497318, 14.76356308), WT_12h__3 = c(14.77215863, 14.72814995, 13.0835214), WT_12h__4 = c(14.70685445, 14.98527337, 16.12727292), WT_24h__1 = c(15.43813077, 14.56918572, 14.92146565 ), WT_24h__2 = c(16.05986898, 14.70583866, 15.64566505), WT_24h__3 = c(14.87721853, 13.22461859, 16.34119942), WT_24h__4 = c(14.92822133, 14.74382383, 12.79146694)), class = “data.frame”, row.names = c(NA, -3L))
我必须总结每个基因（即“ATCG01240”、“ATCG01310”、“ATMG00070”）时间点的数据； 4 次重复的平均值、标准误差并对每个基因、每个时间点进行多重比较（Posthoc 2way ANOVA；即 WT-aox2、WT-aox5、aox2-aox5）。
试一试。由于有关数据帧是列表的错误消息，我最终无法获得方差分析的一种方式。看看您可以用它做什么 - 或者其他人可以提供帮助。
df <- structure(list(AGI = c("ATCG01240", "ATCG01310", "ATMG00070"), aox2_0h__1 = c(15.79105291, 14.82652303, 14.70630068), aox2_0h__2 = c(16.06494674, 14.50610036, 14.52189807), aox2_0h__3 = c(14.64596287, 14.73266459, 13.07143141), aox2_0h__4 = c(15.71713641, 15.15430026, 16.32190068 ), aox2_12h__1 = c(14.99030606, 15.08046949, 15.8317372), aox2_12h__2 = c(15.15569857, 14.98996474, 14.64862254), aox2_12h__3 = c(15.12144791, 14.90111092, 14.59618842), aox2_12h__4 = c(14.25648197, 15.09832061, 14.64442686), aox2_24h__1 = c(15.23997241, 14.80968391, 14.22573239 ), aox2_24h__2 = c(15.57551513, 14.94861669, 15.18808897), aox2_24h__3 = c(15.04928714, 14.83758685, 13.06948037), aox2_24h__4 = c(14.79035385, 14.93873234, 14.70402827), aox5_0h__1 = c(15.8245918, 14.9351844, 14.67678306), aox5_0h__2 = c(15.75108628, 14.85867002, 14.45704948 ), aox5_0h__3 = c(14.36545859, 14.79296855, 14.82177912), aox5_0h__4 = c(14.80626019, 13.43330964, 16.33482718), aox5_12h__1 = c(14.66327372, 15.22571466, 16.17761867), aox5_12h__2 = c(14.58089039, 14.98545497, 14.4331578), aox5_12h__3 = c(14.58091828, 14.86139511, 15.83898617 ), aox5_12h__4 = c(14.48097297, 15.1420725, 13.39369381), aox5_24h__1 = c(15.41855602, 14.9890092, 13.92629626), aox5_24h__2 = c(15.78386057, 15.19372889, 14.63254456), aox5_24h__3 = c(15.55321382, 14.82013321, 15.74324956), aox5_24h__4 = c(14.53085803, 15.12196994, 14.81028556 ), WT_0h__1 = c(14.0535031, 12.45484834, 14.89102226), WT_0h__2 = c(13.64720361, 15.07144643, 14.99836235), WT_0h__3 = c(14.28295759, 13.75283646, 14.98220861), WT_0h__4 = c(14.79637443, 15.1108037, 15.21711524 ), WT_12h__1 = c(15.05711898, 13.33689777, 14.81064042), WT_12h__2 = c(14.83846779, 13.62497318, 14.76356308), WT_12h__3 = c(14.77215863, 14.72814995, 13.0835214), WT_12h__4 = c(14.70685445, 14.98527337, 16.12727292), WT_24h__1 = c(15.43813077, 14.56918572, 14.92146565 ), WT_24h__2 = c(16.05986898, 14.70583866, 15.64566505), WT_24h__3 = c(14.87721853, 13.22461859, 16.34119942), WT_24h__4 = c(14.92822133, 14.74382383, 12.79146694)), class = "data.frame", row.names = c(NA, -3L))
#https://www.r-statistics.com/tag/transpose/
mdf<-melt(df)
names(mdf)<-c("gene","TRT_time","value") #name the new columns
#need to parse TRT_time into separate TRT, time, and sample number columns.
#
foo <- data.frame(do.call('rbind', strsplit(as.character(mdf$TRT_time),'_',fixed=TRUE)))
mdf[2] <- NULL
typeof(foo)
typeof(mdf)
#need to merge these two lists
mdf2<-cbind(foo, mdf)
names(mdf2)<-c("treatment","time","blank", "sampleID","gene","value") #name the new columns
#remove the blank column
mdf2[3] <- NULL
agg.mean <- aggregate(value ~ gene + treatment + time,data=mdf2,mean)
names(agg.mean) <- c("gene","treatment","time", "mean") #name the new columns
#create a function for standard error of the mean
st.err <- function(x) {
sd(x)/sqrt(length(x))
}
st.err <- function(x, na.rm=FALSE) {
if(na.rm==TRUE) x <- na.omit(x)
sd(x)/sqrt(length(x))
}
agg.SEM <- aggregate(value ~ gene + treatment + time,data=mdf2,st.err)
names(agg.SEM) <- c("gene","treatment","time", "SEM") #name the new columns
agg.meanSEM <- cbind (agg.mean, agg.SEM)#prove they match by row
agg.meanSEM <- cbind (agg.mean, agg.SEM$SEM)
#sort the data in mdf2 in preparation for 2 way ANOVA repeated measures
names(mdf2)
mdf3<-mdf2[with(mdf2, order(-sampleID, gene, treatment, time)), ]
#https://rcompanion.org/rcompanion/d_08.html
#Two-way Anova
# An R Companion for the Handbook of Biological Statistics Salvatore S. Mangiafico
install.packages("Rmisc")
library(Rmisc)
sum = summarySE(mdf2, measurevar="value", groupvars=c("gene", "treatment","time"))
sum
#In case you are interested.
boxplot(value ~ gene:treatment:time,
data = mdf2,
xlab = "Gene X Treatment X Time",
ylab = "Gene Expression Level")
#Fit the linear model and conduct ANOVA
#original from source above
# model = lm(value ~ treatment + time + treatment:time,
# data=mdf2)
# Anova(model, type="III")
#I could not get the type added into the split commands; so I do not know which type is being run.
install.packages("car")
library(car)
#using split to 'aggregate' by gene for 2 way Anova.
#
s <- split(mdf2, mdf2$gene)
sink("sink-examp.txt")#this will capture results in a text file to your working directory
do.call(rbind, lapply(s, function(x) Anova(lm(value ~ treatment + time + treatment:time, x))[ ]))
sink()
###############
#Now to arrange for a one way ANOVA to run the TukeyHSD.
df_posthoc<-t(df)
df_posthoc<-data.frame(df_posthoc)
typeof(df_posthoc)
#Convert first row to headers.
#
names(df_posthoc) <- as.matrix(df_posthoc[1, ])
df_posthoc <- df_posthoc[-1, ]
df_posthoc[] <- lapply(df_posthoc, function(x) type.convert(as.character(x)))
df_posthoc
names(df_posthoc)
#
df_posthoc <- cbind(rownames(df_posthoc), data.frame(AD, row.names=NULL))
colnames(df_posthoc)[1] <- "TRT_time" #name the first column
#need to parse TRT_time into separate TRT, time, and sample number columns.
#
foo2 <- data.frame(do.call('rbind', strsplit(as.character(AD$TRT_time),'_',fixed=TRUE)))
foo2[3] <- NULL
colnames(foo2)[1] <- "treatment" #name the first column
colnames(foo2)[2] <- "time" #name the first column
colnames(foo2)[3] <- "sampleID" #name the first column
df_posthoc[1] <- NULL
typeof(foo2)
typeof(df_posthoc)
#need to merge these two lists
df_posthoc<-cbind(foo2, df_posthoc)
typeof(df_posthoc)
AD<-data.frame(df_posthoc)
typeof(df_posthoc)
mode(df_posthoc)
data.matrix(df_posthoc)
typeof(df_posthoc)
mode(df_posthoc)
df_posthoc <- mapply(df_posthoc, FUN=as.numeric)
df_posthoc<-data.frame(df_posthoc)
#Useful site. https://www.statmethods.net/stats/anova.html
1:ncol(df_posthoc)
AVz<- rep(NA,ncol(df_posthoc))# Creates a table, AV, with the same number of columns as in file.
for (i in 2:ncol(df_posthoc)){
AVz<-summary(aov(df_posthoc~treatment, data=df_posthoc))
tk<-TukeyHSD((aov(df_posthoc~treatment, data=df_posthoc)))
capture.output(AVz,file="One-WayANOVASummarycolumn.doc")
capture.output(tk,file="TukeyHSDResultscolumn.doc")
}