-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathAQUACOSM2022_Statistics_Table3.R
145 lines (108 loc) · 4.8 KB
/
AQUACOSM2022_Statistics_Table3.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#_______________________________________________________________________________
#HOUSEHOLD
rm(list=ls())
setwd("/Users/Anika/Library/Mobile Documents/com~apple~CloudDocs/PhD/AquaCosm/SupOpTrons_2022/4 NP Incubations/8 R Scripts")
#Load libraries,
library(readxl)
library(tidyverse)
library(MASS)
library(nlme)
library(reshape2)
library(sjPlot)
library(lattice)
library(mgcv)
library(lme4)
library(tidyverse)
library(psych)
library(ggplot2)
library(effects)
library(dplyr)
library(car)
library(performance)
library(modelsummary)
library(multcomp)
#Read data and prepare for usage for r and N:P
Data_Growth <- read_excel("Excel Files for Stats/STATISTICS_TABLE.xlsx")
Data_Growth$Temp <- as.factor(Data_Growth$Temp)
Data_Growth$Run <- as.factor(Data_Growth$Run)
Data_Growth$T_Treat <- as.factor(Data_Growth$T_Treat)
Data_Growth$T_ID <- as.factor(Data_Growth$T_ID)
Data_Growth <- mutate(Data_Growth, Treat_Unique = paste(T_ID, Run, sep = "_"))
Data_Growth <- mutate(Data_Growth, Balanced = "Balanced")
Data_Growth$Balanced[Data_Growth$Final_Ratio >= 42] <- "P-Limited"
Data_Growth$Balanced[Data_Growth$Final_Ratio <= 11] <- "N-Limited"
Data_Growth <- filter(Data_Growth, Treat_Unique != "Constant_12_1")
Data_Growth <- filter(Data_Growth, Treat_Unique != "Constant_18_1")
#Read data and prepare for usage for LRR
Data_LRR <- read_excel("Excel Files for Stats/SupOpTrons_LRR_NutAv.xlsx")
Data_LRR$Temp <- as.factor(Data_LRR$Temp)
Data_LRR <- mutate(Data_LRR, Mean_Ratio = (Final_Ratio_Run1+Final_Ratio_Run2)/2)
Data_LRR <- mutate(Data_LRR, Balanced = "Balanced")
Data_LRR$Balanced[Data_LRR$Mean_Ratio >= 42] <- "P-Limited"
Data_LRR$Balanced[Data_LRR$Mean_Ratio <= 11] <- "N-Limited"
#_______________________________________________________________________________
#### TABLE 3 GLM: GROWTH ~ TEMP * RATIO * NUTAV ####
#Here the data has been box-cox transformed due to a strong right-shift of the
#data. Afterwards, four outliers have been removed based on the assumption plots.
# 1) Explore your data
Data_Growth_H2 <- mutate(Data_Growth, Growth_Lin_3 = Growth_Lin+1) #Because box-cox transformation does not support negative values
# 2) Check for outliers
#Maybe transform the data since the distribution is right-shifted?
Data_Growth_H2 <- Data_Growth_H2[-171, ]
Data_Growth_H2 <- Data_Growth_H2[-153, ]
Data_Growth_H2 <- Data_Growth_H2[-171, ]
Data_Growth_H2 <- Data_Growth_H2[-240, ]
par(mfrow = c(1, 3)) #Three figures in one graph
boxplot((Data_Growth_H2$Growth_Lin_3)^3) #Does it look normally distributed?
hist((Data_Growth_H2$Growth_Lin_3)^3) #Does it look normally distributed?
dotchart((Data_Growth_H2$Growth_Lin_3)^3) #Points very far away are outliers
# 6) Building the model
#T*N Interaction with all levels
glm.growth.TNP <- glm((Growth_Lin_3)^3 ~ Temp * Balanced * Run, data = Data_Growth_H2, family = gaussian)
plot(glm.growth.TNP) #Check the assumptions of the model!
hist(resid(glm.growth.TNP)) #Does the data look normally distributed?
#summary(glm.growth.TNP)
summary(aov(glm.growth.TNP))
#This is to test which exponent works best!
bc = boxcox(glm.growth.TNP, lambda = seq(-4,4))
best.lam <- bc$x[which(bc$y==max(bc$y))]
best.lam #This gives you the best exponent!
#_______________________________________________________________________________
#### TABLE 3 GLM: LRR ~ TEMP * RATIO ####
# 2) Check for outliers
Data_LRR <- Data_LRR[-1, ]
Data_LRR <- Data_LRR[-70, ]
Data_LRR <- Data_LRR[-69, ]
Data_LRR <- Data_LRR[-1, ]
par(mfrow = c(1, 3)) #Three figures in one graph
boxplot(Data_LRR$LRR) #Does it look normally distributed?
hist(Data_LRR$LRR) #Does it look normally distributed?
dotchart(Data_LRR$LRR) #Points very far away are outliers
# 6) Building the model
#Temp + Run
#T*N Interaction with all levels
glm_LRR <- glm(LRR ~ (Temp * Balanced), family = gaussian, data = Data_LRR)
plot(glm_LRR) #Check the assumptions of the model!
hist(resid(glm_LRR)) #Does the data look normally distributed?
summary(aov(glm_LRR))
#_______________________________________________________________________________
#### TABLE 3 GLM: N:P ~ TEMP * RATIO * NUTAV ####
Data_NP_2 <- Data_Growth
# 2) Check for outliers
Data_NP_2 <- Data_NP_2[-126, ]
Data_NP_2 <- Data_NP_2[-14, ]
Data_NP_2 <- Data_NP_2[-76, ]
Data_NP_2 <- Data_NP_2[-75, ]
Data_NP_2 <- Data_NP_2[-147, ]
par(mfrow = c(1, 3)) #Three figures in one graph
boxplot(Data_NP_2$NP_Ratio) #Does it look normally distributed?
hist(Data_NP_2$NP_Ratio) #Does it look normally distributed?
dotchart(Data_NP_2$NP_Ratio) #Points very far away are outliers
# 6) Building the model
#Temp + Run
#T*N Interaction with all levels
glm_3 <- glm(log(NP_Ratio) ~ (Temp * Balanced * Run), family = gaussian, data = Data_NP_2)
plot(glm_3) #Check the assumptions of the model!
hist(resid(glm_3)) #Does the data look normally distributed?
summary(aov(glm_3))
#_______________________________________________________________________________