曾老师钩栗数据分析1229

d1-d3, h1-h3分别代表3年的地径和株高生长量。

1.种源间的差异（3个地点）

数据结构:

gouli3p20161228 <- read.csv("D:/R/zeng/gouli3p20161228.csv")
str(gouli3p20161228)

## 'data.frame':    554 obs. of  8 variables:
##  $ Provenance: Factor w/ 3 levels "AH","LQ","TS": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Family    : Factor w/ 18 levels "AH_1","AH_10",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ d1        : num  2.83 2.99 4.19 3.91 2.6 3.44 3.92 3.57 5.95 2.47 ...
##  $ h1        : num  5.5 5.4 15.2 17.1 8.7 13.1 10.3 9.8 14.5 10.1 ...
##  $ d2        : num  6.97 4.76 7.97 9.52 6.11 6.87 6.03 7.61 7.25 7.5 ...
##  $ h2        : num  50.6 23 48 71.2 46.5 54.2 41 67 52.3 38.4 ...
##  $ d3        : num  6.86 5.59 8.6 13.5 6.56 11.4 9.2 8.9 7.39 8.9 ...
##  $ h3        : num  53.1 52.1 74.2 85.9 51.3 59.8 59.2 79.2 68.5 48.9 ...

对数据进行探索分析、作图并进行方差分析。

1. 探索分析

library(psych)

# 变量名给X
X <- names(gouli3p20161228)
# 建一个空白矩阵
result.discribe <- matrix(data = NA, nrow = length(X) - 2, ncol = 13, byrow = T)
rownames(result.discribe) <- X[3:length(X)]

# 响应变量在第3列
for (i in 3:length(X)) {
    data.i <- gouli3p20161228[, X[i]]
    # 单个变量的统计结果
    temp.des <- describe(data.i)
    # 讲结果给空矩阵
    result.discribe[c(i - 2), ] <- as.matrix(temp.des[1, ])
}

# 列名给矩阵列名
colnames(result.discribe) <- colnames(temp.des)
options(digits = 3)
knitr::kable(result.discribe)

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
d1	1	554	3.51	0.974	3.52	3.49	0.971	1.33	7.13	5.80	0.300	0.383	0.041
h1	1	554	14.41	6.505	13.00	13.73	5.930	3.20	37.00	33.80	0.947	0.516	0.276
d2	1	554	7.12	1.750	7.09	7.16	1.824	2.17	11.73	9.56	-0.182	-0.321	0.074
h2	1	554	53.90	18.411	54.20	53.65	18.088	2.80	103.40	100.60	0.101	-0.191	0.782
d3	1	456	9.51	2.195	9.37	9.42	2.157	4.34	21.50	17.16	0.589	1.473	0.103
h3	1	456	73.46	22.121	71.65	72.49	21.572	23.50	161.00	137.50	0.477	0.427	1.036

2. 作图

library(multcompView)
# I need to group the treatments that are not different each other together.
generate_label_df <- function(TUKEY, variable) {
    # Extract labels and factor levels from Tukey post-hoc
    Tukey.levels <- TUKEY[[variable]][, 4]
    # multcompLetters这个函数确定了字母
    Tukey.labels <- data.frame(multcompLetters(Tukey.levels, threshold = 0.05)["Letters"])
    # I need to put the labels in the same order as in the boxplot :
    Tukey.labels$treatment = rownames(Tukey.labels)
    Tukey.labels = Tukey.labels[order(Tukey.labels$treatment), ]
    return(Tukey.labels)
}

par(mfrow=c(3,2),family="serif",mar = c(3, 4, .5, .5))
X<-names(gouli3p20161228) 
    for (i in 3:length(X)){
      data.i<-gouli3p20161228[, X[i]]
      model=lm(data.i~gouli3p20161228$Provenance)
      ANOVA=aov(model)
      TUKEY <- TukeyHSD(x=ANOVA, 'gouli3p20161228$Provenance', conf.level=0.95)
      # 保存p值
      # a<-TUKEY$`gouli3p20161228$Provenance`[,c(1,4)]
      # colnames(a)<-c(paste(X[j],"diff",sep = "_"),"p adj")
      # p.file<-paste0(filename,"_pvalue.csv")
      # write.table(a,file=p.file,sep=",",quote=F,append=T)
      LABELS=generate_label_df(TUKEY, "gouli3p20161228$Provenance")
      a=boxplot(data.i~gouli3p20161228$Provenance,
                ylim=c(min(data.i,na.rm=T) ,
                       1.01*max(data.i,na.rm=T)) , 
                #xlab="Crosses",
                ylab=X[i] , main=""
                # 箱体边界的颜色
                # border=c("grey","white")
                )
      abline(h=mean(data.i,na.rm=T),lty=2)
      means <- aggregate(data.i, list(gouli3p20161228$Provenance), mean, na.rm=T)
      points(means$x, col = "black", pch = 18)
      # I want to write the letter over each box. Over is how high I want to write it.
      over=0.04*max( a$stats[nrow(a$stats),] )
      #Add the labels
      text( c(1:nlevels(gouli3p20161228$Provenance)) , 
            a$stats[nrow(a$stats),]+over , LABELS[,1])
      
      # 普通图
      # boxplot(data.i~Provenance,
      #         xlab="Provenance", ylab=X[i],
      #         data = gouli3p20161228)
    }

3. 方差分析

summary(aov(d1~Provenance,data = gouli3p20161228))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## Provenance    2      6    2.97    3.16  0.043 *
## Residuals   551    518    0.94                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(h1~Provenance,data = gouli3p20161228))

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Provenance    2   8703    4352     163 <2e-16 ***
## Residuals   551  14697      27                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(d2~Provenance,data = gouli3p20161228))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## Provenance    2     22   11.20    3.69  0.025 *
## Residuals   551   1670    3.03                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(h2~Provenance,data = gouli3p20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    2   9571    4786    14.8 5.4e-07 ***
## Residuals   551 177880     323                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(d3~Provenance,data = gouli3p20161228))

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Provenance    2    226   112.9      26  2e-11 ***
## Residuals   453   1967     4.3                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 98 observations deleted due to missingness

summary(aov(h3~Provenance,data = gouli3p20161228))

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Provenance    2  21614   10807    24.4  9e-11 ***
## Residuals   453 201027     444                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 98 observations deleted due to missingness

2 种源与系的差异(2个地点)

数据结构:

gouli20161228 <- read.csv("D:/R/zeng/gouli20161228.csv")
str(gouli20161228)

## 'data.frame':    528 obs. of  8 variables:
##  $ Provenance: Factor w/ 2 levels "AH","LQ": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Family    : Factor w/ 16 levels "AH_1","AH_10",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ d1        : num  2.83 2.99 4.19 3.91 2.6 3.44 3.92 3.57 5.95 2.47 ...
##  $ h1        : num  5.5 5.4 15.2 17.1 8.7 13.1 10.3 9.8 14.5 10.1 ...
##  $ d2        : num  6.97 4.76 7.97 9.52 6.11 6.87 6.03 7.61 7.25 7.5 ...
##  $ h2        : num  50.6 23 48 71.2 46.5 54.2 41 67 52.3 38.4 ...
##  $ d3        : num  6.86 5.59 8.6 13.5 6.56 11.4 9.2 8.9 7.39 8.9 ...
##  $ h3        : num  53.1 52.1 74.2 85.9 51.3 59.8 59.2 79.2 68.5 48.9 ...

对数据进行探索分析、作图并进行方差分析。

1. 探索分析

library(psych)

# 变量名给X
X<-names(gouli20161228)
# 建一个空白矩阵
result.discribe<-matrix(data = NA, nrow = length(X)-2, ncol = 13, byrow = T)
rownames(result.discribe)<-X[3:length(X)]

# 响应变量在第3列
    for (i in 3:length(X)){
      data.i<-gouli20161228[,  X[i]]
      # 单个变量的统计结果
      temp.des<-describe(data.i)
      # 讲结果给空矩阵
      result.discribe[c(i-2),]<-as.matrix(temp.des[1,])
    }

# 列名给矩阵列名
colnames(result.discribe)<-colnames(temp.des)
options(digits = 3)
knitr::kable(result.discribe)

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
d1	1	528	3.49	0.974	3.49	3.47	0.964	1.33	7.13	5.80	0.362	0.480	0.042
h1	1	528	14.23	6.483	12.90	13.50	5.782	3.20	37.00	33.80	1.025	0.730	0.282
d2	1	528	7.10	1.739	7.07	7.14	1.801	2.17	11.73	9.56	-0.162	-0.284	0.076
h2	1	528	54.00	18.400	54.30	53.73	18.014	2.80	103.40	100.60	0.111	-0.146	0.801
d3	1	433	9.46	2.093	9.36	9.40	2.165	4.65	16.40	11.75	0.294	-0.125	0.101
h3	1	433	73.70	21.868	71.70	72.68	20.756	24.40	161.00	136.60	0.516	0.508	1.051

2. 作图

par(mfrow=c(3,2),family="serif",mar = c(4, 4, .5, .5))
X<-names(gouli20161228) 
    for (i in 3:length(X)){
      data.i<-gouli20161228[, X[i]]
      boxplot(data.i~Family,
              xlab="Family", ylab=X[i],
              col = c(rep("grey",10),rep("white",6)),
              data = gouli20161228,
              las=2)
    }

3. 方差分析

summary(aov(d1~Provenance*Family,data = gouli20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    1      1    0.99    1.19    0.28    
## Family       14     75    5.34    6.44 5.1e-12 ***
## Residuals   512    424    0.83                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(h1~Provenance*Family,data = gouli20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    1   8321    8321  350.89 < 2e-16 ***
## Family       14   1683     120    5.07 6.1e-09 ***
## Residuals   512  12142      24                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(d2~Provenance*Family,data = gouli20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    1     18   18.49    7.06  0.0081 ** 
## Family       14    234   16.70    6.37 7.2e-12 ***
## Residuals   512   1342    2.62                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(h2~Provenance*Family,data = gouli20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    1   9456    9456   35.80 4.1e-09 ***
## Family       14  33729    2409    9.12 < 2e-16 ***
## Residuals   512 135239     264                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aov(d3~Provenance*Family,data = gouli20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    1    204   203.8   55.43 5.6e-13 ***
## Family       14    156    11.2    3.03 0.00018 ***
## Residuals   417   1533     3.7                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 95 observations deleted due to missingness

summary(aov(h3~Provenance*Family,data = gouli20161228))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Provenance    1  21106   21106   55.71 4.9e-13 ***
## Family       14  27497    1964    5.18 4.7e-09 ***
## Residuals   417 157977     379                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 95 observations deleted due to missingness

3.泰顺不同小群落的差异

数据结构:

taishun20161228 <- read.csv("D:/R/zeng/taishun20161228.csv")
str(taishun20161228)

## 'data.frame':    26 obs. of  8 variables:
##  $ Provenance: Factor w/ 1 level "TS": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Population: Factor w/ 2 levels "HLS","YS": 2 2 2 2 2 2 2 2 2 2 ...
##  $ d1        : num  4.85 4.64 4.59 4.36 3.19 3.64 4.44 3.32 1.89 4.41 ...
##  $ h1        : num  23 17 18 19 15 15 25 18 9 19 ...
##  $ d2        : num  9.36 7.99 10.04 4.92 8.48 ...
##  $ h2        : num  54.7 63.1 63.4 26.7 72.8 36.2 76.7 74.5 31 69.2 ...
##  $ d3        : num  10.3 NA 10.6 13.2 14.5 ...
##  $ h3        : num  65.3 NA 72.2 43.4 89.2 ...

对数据进行探索分析、作图并进行方差分析。

1. 探索分析

library(psych)

# 变量名给X
X<-names(taishun20161228)
# 建一个空白矩阵
result.discribe<-matrix(data = NA, nrow = length(X)-2, ncol = 13, byrow = T)
rownames(result.discribe)<-X[3:length(X)]

# 响应变量在第3列
    for (i in 3:length(X)){
      data.i<-taishun20161228[,  X[i]]
      # 单个变量的统计结果
      temp.des<-describe(data.i)
      # 讲结果给空矩阵
      result.discribe[c(i-2),]<-as.matrix(temp.des[1,])
    }

# 列名给矩阵列名
colnames(result.discribe)<-colnames(temp.des)
options(digits = 3)
knitr::kable(result.discribe)

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
d1	1	26	3.94	0.876	4.25	4.04	0.726	1.64	5.16	3.52	-1.049	0.445	0.172
h1	1	26	18.15	5.904	18.50	18.45	7.413	5.00	26.00	21.00	-0.370	-0.928	1.158
d2	1	26	7.50	1.949	7.85	7.60	2.031	3.53	10.40	6.87	-0.566	-0.763	0.382
h2	1	26	51.84	18.882	52.00	52.08	25.130	19.50	82.00	62.50	-0.088	-1.329	3.703
d3	1	23	10.46	3.552	9.67	10.20	2.461	4.34	21.50	17.16	1.115	1.849	0.741
h3	1	23	68.88	26.589	68.50	68.37	28.318	23.50	123.50	100.00	0.197	-0.904	5.544

2. 作图

par(mfrow=c(3,2),family="serif",mar = c(3, 4, .5, .5))
X<-names(taishun20161228) 
    for (i in 3:length(X)){
      data.i<-taishun20161228[, X[i]]
      boxplot(data.i~Population,
              xlab="Population", ylab=X[i],
              data = taishun20161228)
    }

3. 方差分析

summary(aov(d1~Population,data = taishun20161228))

##             Df Sum Sq Mean Sq F value Pr(>F)
## Population   1   0.88   0.881    1.15   0.29
## Residuals   24  18.32   0.763

summary(aov(h1~Population,data = taishun20161228))

##             Df Sum Sq Mean Sq F value Pr(>F)
## Population   1     65    65.0    1.93   0.18
## Residuals   24    806    33.6

summary(aov(d2~Population,data = taishun20161228))

##             Df Sum Sq Mean Sq F value Pr(>F)
## Population   1    0.9    0.88    0.22   0.64
## Residuals   24   94.1    3.92

summary(aov(h2~Population,data = taishun20161228))

##             Df Sum Sq Mean Sq F value Pr(>F)
## Population   1    620     620    1.79   0.19
## Residuals   24   8293     346

summary(aov(d3~Population,data = taishun20161228))

##             Df Sum Sq Mean Sq F value Pr(>F)
## Population   1   15.4    15.4    1.24   0.28
## Residuals   21  262.1    12.5               
## 3 observations deleted due to missingness

summary(aov(h3~Population,data = taishun20161228))

##             Df Sum Sq Mean Sq F value Pr(>F)
## Population   1    614     614    0.86   0.36
## Residuals   21  14939     711               
## 3 observations deleted due to missingness

保存工作对象

save(gouli20161228,gouli3p20161228,taishun20161228,file = "zeng1229.rdata")

曾老师钩栗数据分析1229

Leiming Dong

2016年12月29日

1.种源间的差异（3个地点）

1. 探索分析

2. 作图

3. 方差分析

2 种源与系的差异(2个地点)

1. 探索分析

2. 作图

3. 方差分析

3.泰顺不同小群落的差异

1. 探索分析

2. 作图

3. 方差分析