因為要去掉數(shù)據(jù)中的自交四瘫，所以之前畫的boxplot需要重新做，代碼如下：
setwd("D:/Sync/R/8801/fullsib/160620")
Myfile<-c("8801fs_05_14.csv")
df<-read.csv(Myfile)
str(df)
# 'data.frame': 2226 obs. of 20 variables:
# $ TreeId : int 3141 3142 3143 3144 3145 3146 3147 3148 3149 3150 ...
# $ Mum : Factor w/ 23 levels "CF29","R1","R1_OP",..: 18 18 18 18 18 18 18 18 18 18 ...
# $ Dad : Factor w/ 13 levels "","CF29","R1",..: 8 8 8 8 8 8 8 8 8 8 ...
# $ Com : Factor w/ 66 levels "CF29RC402","CF29Y11",..: 52 52 52 52 52 52 52 52 52 52 ...
# $ Tree : int NA 123 NA 107 NA 120 104 122 118 121 ...
# $ Fam : int 84 84 84 84 84 84 84 84 84 84 ...
# $ Rep : int 1 1 1 1 1 1 1 1 1 1 ...
# $ Row : int 41 41 41 42 42 42 43 43 44 44 ...
# $ Col : int 52 53 55 51 54 55 52 54 51 53 ...
# $ Recipro: Factor w/ 18 levels "",".","RC13RC82",..: 12 12 12 12 12 12 12 12 12 12 ...
# $ Plot : Factor w/ 250 levels "","1CF29Y11",..: 42 42 42 42 42 42 42 42 42 42 ...
# $ h05 : num 15.5 15 15.5 16.5 17.5 15.8 17.2 17.2 15.5 16.1 ...
# $ d05 : num 14.5 16 13 13.3 13.4 14.7 18.3 17.3 14.5 14.8 ...
# $ v05 : num 0.326 0.384 0.262 0.292 0.314 ...
# $ h14 : num NA 19.9 NA 20.7 NA 20.8 23.2 23.1 19.2 21.6 ...
# $ d14 : num NA 20.2 NA 15.9 NA 19.1 23.1 24 17.1 17.4 ...
# $ v14 : num NA 0.812 NA 0.523 NA ...
# $ delta.h: num NA 4.9 NA 4.2 NA 5 6 5.9 3.7 5.5 ...
# $ delta.d: num NA 4.2 NA 2.6 NA 4.4 4.8 6.7 2.6 2.6 ...
# $ delta.v: num NA 0.428 NA 0.231 NA ...
names(df)
# [1] "TreeId" "Mum" "Dad" "Com" "Tree" "Fam" "Rep" "Row"
# [9] "Col" "Recipro" "Plot" "h05" "d05" "v05" "h14" "d14"
# [17] "v14" "delta.h" "delta.d" "delta.v"
names(df)[12:length(names(df))]<-c("H05","DBH05","V05","H14","DBH14","V14","HI", "DI","VI")
df$Mum<-as.character(df$Mum)
df$Dad<-as.character(df$Dad)
# df$Mum<-factor(df$Mum)
# df$Dad<-factor(df$Dad)
df.rm.self<-df[df$Mum!=df$Dad,]
levels(df.rm.self$Com)
str(df.rm.self)
df.rm.self[df.rm.self$Com=="Y9Y9",]
df.rm.self$Com<-as.character(df.rm.self$Com)
df.rm.self$Com<-factor(df.rm.self$Com)
# par(las=2)
# plot(df.rm.self$H05~df.rm.self$Com)
source("D:/Sync/R/rcodes/d.box.order.r")
names(df.rm.self)[4]<-c("com")
#pdf
pdf("grow.all.crosses.order.pdf", width = 7, height = 4,family="Times",pointsize = 8)
par(mar = c(6, 4, .5, .5),las=2)
X<-names(df.rm.self)
for (j in 12:length(X)){
data.i<-df.rm.self[, X[j]]
d.box.order(data.i, df.rm.self, X, j)
}
dev.off()
# win.metafile
win.metafile("grow.all.crosses.H05.wmf", width = 7, height = 4,family="serif",pointsize = 8)
par(mar = c(6, 4, .5, .5),las=2)
X<-names(df.rm.self)
data.i<-df.rm.self[, X[12]]
d.box.order(data.i, df.rm.self, X, 12)
dev.off()
com<-unique(df$Com)
com<-as.character(com)
comcom<-read.csv("comcom.growth.csv")
# comcom<-read.csv("comcom.growth.rmself.csv")
comcom<-as.matrix(comcom)
names(df)[4]<-c("com")
pdf("grow.performance.rmself.pdf",paper="a4", width = 0, height = 9,family="Times",pointsize = 10)
par(mfrow=c(3,3),mar = c(6, 4, .5, .5),las=2)
for (i in 1:nrow(comcom)){
df.com<-df[df$com %in% com[comcom[i,]],]
df.com$com<-factor(df.com$com)
X<-names(df.com)
for (j in 12:length(X)){
data.i<-df.com[, X[j]]
model=lm(data.i ~ df.com$com )
ANOVA=aov(model)
TUKEY <- TukeyHSD(x=ANOVA, 'df.com$com', conf.level=0.99)
# a<-TUKEY$`df.com$com`[,c(1,4)]
# colnames(a)<-c(paste(X[j],"diff",sep = "_"),"p adj")
# write.table(a,file="growth_p_value.csv",sep=",",quote=F,append=T)
LABELS=generate_label_df(TUKEY , "df.com$com")
a=boxplot(data.i~df.com$com,
ylim=c(min(data.i,na.rm=T) ,
1.01*max(data.i,na.rm=T)) ,
#xlab="Crosses",
ylab=X[j] , main="")
abline(h=mean(data.i,na.rm=T),lty=2)
# I want to write the letter over each box. Over is how high I want to write it.
over=0.02*max( a$stats[nrow(a$stats),] )
#Add the labels
text( c(1:nlevels(df.com$com)) ,
a$stats[nrow(a$stats),]+over , LABELS[,1])
}}
dev.off()
# win.metafile
win.metafile("grow.performance.rmself.wmf", width = 8.3, height = 11.7,family="serif",pointsize = 10)
par(mfrow=c(3,3),mar = c(6, 4, .5, .5),las=2)
df.com<-df[df$com %in% com[comcom[1,]],]
df.com$com<-factor(df.com$com)
X<-names(df.com)
for (j in 12:length(X)){
data.i<-df.com[, X[j]]
model=lm(data.i ~ df.com$com )
ANOVA=aov(model)
TUKEY <- TukeyHSD(x=ANOVA, 'df.com$com', conf.level=0.99)
# a<-TUKEY$`df.com$com`[,c(1,4)]
# colnames(a)<-c(paste(X[j],"diff",sep = "_"),"p adj")
# write.table(a,file="growth_p_value.csv",sep=",",quote=F,append=T)
LABELS=generate_label_df(TUKEY , "df.com$com")
a=boxplot(data.i~df.com$com,
ylim=c(min(data.i,na.rm=T) ,
1.01*max(data.i,na.rm=T)) ,
#xlab="Crosses",
ylab=X[j] , main="")
abline(h=mean(data.i,na.rm=T),lty=2)
# I want to write the letter over each box. Over is how high I want to write it.
over=0.02*max( a$stats[nrow(a$stats),] )
#Add the labels
text( c(1:nlevels(df.com$com)) ,
a$stats[nrow(a$stats),]+over , LABELS[,1])
}
dev.off()
save(df,df.rm.self,file = "8801.未替換家系名及去掉自交.RData")
Sys.time()
# [1] "2016-10-20 22:14:09 CST"
