生物统计学作业三
Liang
2023-05-08
T 检验
df <- read.csv("Data1.csv")
# 叶片光合速率
ans <- t.test(A ~ Group, data = df)
# draw
dfm <- as.data.frame(list(Group = c("单叶", "复叶"), mean = c(mean(df$A[df$Group == 1]), mean(df$A[df$Group == 2])), se = c(se(df$A[df$Group == 1]), se(df$A[df$Group == 2]))))
p <- ggplot(dfm, aes(Group, mean)) +
geom_col(aes(fill = Group), width = 0.3) +
theme_classic() +
ylab("叶片光合速率") +
xlab("") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.1))+
annotate(geom = "text", x = 1, y = 14, label = "a") +
annotate(geom = "text", x = 2, y = 23, label = "b") +
annotate(geom = "text", x = 1, y = 23, label = paste("p-value =", round(ans$p.value, 3)), size = 5)
ans##
## Welch Two Sample t-test
##
## data: A by Group
## t = -2.9671, df = 3.7333, p-value = 0.045
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -12.0068915 -0.2281085
## sample estimates:
## mean in group 1 mean in group 2
## 10.4825 16.6000p
图一: 不同叶形树木的叶片光合速率比较。以单叶和复叶树的叶片光合速率没有差异为假设做t检验，拒绝原假设的概率为0.045，在95%的置信水平下可拒绝原假设。我们的证据支持叶片光合速率在单叶和复叶树之间有区别，复叶树的光合速率更高。误差线示正负标准误。
# 气孔导度
ans <- t.test(gs ~ Group, data = df)
# draw
dfm <- as.data.frame(list(Group = c("单叶", "复叶"), mean = c(mean(df$gs[df$Group == 1]), mean(df$gs[df$Group == 2])), se = c(se(df$gs[df$Group == 1]), se(df$gs[df$Group == 2]))))
p <- ggplot(dfm, aes(Group, mean)) +
geom_col(aes(fill = Group), width = 0.3) +
theme_classic() +
ylab("气孔导度") +
xlab("") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.1))+
annotate(geom = "text", x = 1, y = 0.2, label = "a") +
annotate(geom = "text", x = 2, y = 0.6, label = "b") +
annotate(geom = "text", x = 1, y = 0.75, label = paste("p-value =", round(ans$p.value, 3)), size = 5)
ans##
## Welch Two Sample t-test
##
## data: gs by Group
## t = -4.0419, df = 3.0299, p-value = 0.02675
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -0.4903254 -0.0596746
## sample estimates:
## mean in group 1 mean in group 2
## 0.1275 0.4025p
图二: 不同叶形树木的叶片气孔导读比较。以单叶和复叶树的叶片气孔导度没有差异为假设做t检验，拒绝原假设的概率为0.027，在95%的置信水平下可拒绝原假设。我们的证据支持叶片气孔导度在单叶和复叶树之间有区别，复叶树的气孔导度更高。误差线示正负标准误。
# 水分利用效率
ans <- t.test(WUEi ~ Group, data = df)
ans##
## Welch Two Sample t-test
##
## data: WUEi by Group
## t = 4.964, df = 3.6839, p-value = 0.009489
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## 18.0925 67.8075
## sample estimates:
## mean in group 1 mean in group 2
## 85.15 42.20# draw
dfm <- as.data.frame(list(Group = c("单叶", "复叶"), mean = c(mean(df$WUEi[df$Group == 1]), mean(df$WUEi[df$Group == 2])), se = c(se(df$WUEi[df$Group == 1]), se(df$WUEi[df$Group == 2]))))
p <- ggplot(dfm, aes(Group, mean)) +
geom_col(aes(fill = Group), width = 0.3) +
theme_classic() +
ylab("水分利用效率") +
xlab("") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.1))+
annotate(geom = "text", x = 1, y = 110, label = "a") +
annotate(geom = "text", x = 2, y = 52, label = "b") +
annotate(geom = "text", x = 2, y = 100, label = paste("p-value =", round(ans$p.value, 3)), size = 5)
p
图三: 不同叶形树木的叶片水分利用效率比较。以单叶和复叶树的叶片水分利用效率没有差异为假设做t检验，拒绝原假设的概率为0.009，在95%的置信水平下可拒绝原假设。我们的证据支持叶片水分利用效率在单叶和复叶树之间有区别，单叶树的水分利用效率更高。误差线示正负标准误。
单因素方差分析
df <- read.csv("Data11.csv")
ans <- summary(aov(Result ~ Treat, data = df))
# draw
theLevel <- levels(factor(df$Treat))
dfm <- as.data.frame(list(Treat = theLevel, mean = sapply(theLevel, function(x){mean(df$Result[df$Treat == x])}), se = sapply(theLevel, function(x){se(df$Result[df$Treat == x])})))
p <- ggplot(dfm, aes(Treat, mean)) +
geom_col(aes(fill = Treat), width = 0.3) +
theme_classic() +
ylab("疗效") +
xlab("") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.1))+
annotate(geom = "text", x = 1, y = 450, label = paste("p-value =", round(ans[[1]]$`Pr(>F)`[1], 3)), size = 5)
ans## Df Sum Sq Mean Sq F value Pr(>F)
## Treat 3 137206 45735 4.782 0.0114 *
## Residuals 20 191288 9564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p
图四: 四个治疗组别的疗效，误差线示正负标准误。采用固定效应模型，以四个治疗组之间疗效不存在差异为原假设，组别为因素作方差分析。拒绝原假设的概率为0.011，在95%的置信水平下，支持至少有两组之间的疗效是不一样的。
简单线性回归
先分别采用混合效应模型两两分析单叶（SL）和复叶树（CL）的叶片光合速率（A）、气孔导度（gs）与整枝导水率（Kshoot）的关系，由叶型提供随机截距。
df <- read.csv("Data1.csv")
ans <- lmer(Kshoot ~ A + (1 | Leaf_type), data = df)
random <- ranef(ans)$Leaf_type$`(Intercept)`
ans <- summary(ans)
ans## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Kshoot ~ A + (1 | Leaf_type)
## Data: df
##
## REML criterion at convergence: -8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.62060 -0.38744 -0.00878 0.57294 1.19631
##
## Random effects:
## Groups Name Variance Std.Dev.
## Leaf_type (Intercept) 0.010602 0.10297
## Residual 0.003645 0.06038
## Number of obs: 8, groups: Leaf_type, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.329918 0.132274 3.850043 -2.494 0.069610 .
## A 0.053349 0.008001 5.891683 6.667 0.000594 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## A -0.819ggplot(df, aes(A, Kshoot)) +
geom_point(aes(color = Leaf_type, size = 10))+
geom_abline(intercept = ans[[10]][1,1], slope = ans[[10]][2,1], linewidth = 1)+
geom_abline(intercept = ans[[10]][1,1] + random[1], slope = ans[[10]][2,1], color = "red", linewidth = 1)+
geom_abline(intercept = ans[[10]][1,1] + random[2], slope = ans[[10]][2,1], color = "#00bfc4", linewidth = 1)+
theme_classic()+
annotate(geom = "text", x = 12, y = 0.75, label = paste("p-value =", round(ans[[10]][2,5], 3)), size = 5)+
xlab("叶片光合速率") + ylab("整枝导水率")+guides(size="none")
图五: 叶片光合速率和整枝导水率之间的正相关关系。由叶片光合速率提供的斜率为零的概率为0.001，支持叶片光合速率和整枝导水率存在强正相关关系。黑线是对所有点的回归线，红色和蓝色的线则是在组内作回归。
ans <- lmer(Kshoot ~ gs + (1 | Leaf_type), data = df)## boundary (singular) fit: see help('isSingular')random <- ranef(ans)$Leaf_type$`(Intercept)`
ans <- summary(ans)
ans## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Kshoot ~ gs + (1 | Leaf_type)
## Data: df
##
## REML criterion at convergence: -13
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4209 -0.5412 0.1244 0.5904 1.3341
##
## Random effects:
## Groups Name Variance Std.Dev.
## Leaf_type (Intercept) 0.000000 0.00000
## Residual 0.006167 0.07853
## Number of obs: 8, groups: Leaf_type, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.04114 0.05353 6.00000 -0.769 0.471
## gs 1.63636 0.17269 6.00000 9.476 7.87e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## gs -0.855
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')ggplot(df, aes(gs, Kshoot)) +
geom_point(aes(color = Leaf_type, size = 10))+
geom_abline(intercept = ans[[10]][1,1], slope = ans[[10]][2,1], linewidth = 1)+
theme_classic()+
annotate(geom = "text", x = 0.25, y = 0.75, label = paste("p-value =", round(ans[[10]][2,5], 3)), size = 5)+
xlab("气孔导度") + ylab("整枝导水率")+guides(size="none")
图六: 叶片气孔导度和整枝导水率之间的正相关关系。由叶片气孔导度提供的斜率为零的概率接近零，支持叶片气孔导度和整枝导水率存在强正相关关系。此处由叶型提供的随机截距为零。
ans <- lmer(gs ~ A + (1 | Leaf_type), data = df)
random <- ranef(ans)$Leaf_type$`(Intercept)`
ans <- summary(ans)
ans## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: gs ~ A + (1 | Leaf_type)
## Data: df
##
## REML criterion at convergence: -9.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2388 -0.8904 0.3619 0.5953 1.0584
##
## Random effects:
## Groups Name Variance Std.Dev.
## Leaf_type (Intercept) 0.003558 0.05965
## Residual 0.003305 0.05749
## Number of obs: 8, groups: Leaf_type, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.160270 0.107232 3.635976 -1.495 0.21622
## A 0.031405 0.007124 5.520204 4.408 0.00552 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## A -0.900ggplot(df, aes(A, gs)) +
geom_point(aes(color = Leaf_type, size = 10))+
geom_abline(intercept = ans[[10]][1,1], slope = ans[[10]][2,1], linewidth = 1)+
geom_abline(intercept = ans[[10]][1,1] + random[1], slope = ans[[10]][2,1], color = "red", linewidth = 1)+
geom_abline(intercept = ans[[10]][1,1] + random[2], slope = ans[[10]][2,1], color = "#00bfc4", linewidth = 1)+
theme_classic()+
annotate(geom = "text", x = 12, y = 0.5, label = paste("p-value =", round(ans[[10]][2,5], 3)), size = 5)+
xlab("叶片光合速率") + ylab("气孔导度")+guides(size="none")
图七: 叶片光合速率和气孔导度之间的正相关关系。由叶片光合速率提供的斜率为零的概率为0.006，支持叶片光合速率和整枝导水率存在强正相关关系。黑线是对所有点的回归线，红色和蓝色的线则是在组内作回归。
进一步分析
我们发现这三个变量之间两两有正相关关系，我们猜测这三个变量之间存在极强的共线性，并且可能是由一个外部变量共同驱动的。我们用混合效应模型同时检验叶片光合速率和气孔导度和它们的交互效应对整枝导水率的影响，以叶型提供随机截距。
ans <- lmer(Kshoot ~ A + gs + A * gs + (1 | Leaf_type), data = df)
ans <- summary(ans)
ans## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Kshoot ~ A + gs + A * gs + (1 | Leaf_type)
## Data: df
##
## REML criterion at convergence: -7.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.03624 -0.45655 0.06407 0.38198 1.14291
##
## Random effects:
## Groups Name Variance Std.Dev.
## Leaf_type (Intercept) 0.115909 0.34045
## Residual 0.002383 0.04882
## Number of obs: 8, groups: Leaf_type, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.34094 0.84774 1.02430 -1.582 0.355
## A 0.08649 0.03871 1.55934 2.234 0.190
## gs 5.52485 3.67791 1.18559 1.502 0.345
## A:gs -0.21547 0.15749 1.21826 -1.368 0.370
##
## Correlation of Fixed Effects:
## (Intr) A gs
## A -0.933
## gs -0.941 0.913
## A:gs 0.951 -0.946 -0.995我们发现，在混合效应模型下，主要的效应并未归结为任何一个变量提供的，这可能是我们的变量的共线性极强，模型无法分辨这两个变量的作用。故我们的最终结论为这三个变量之间存在极强的线性正相关关系，我们可以用这其中某个变量去预测其他变量。但其之间是否存在因果关系，以及其是否由某个外部的变量共同驱动需要更多的研究。