-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathAQUACOSM2022_Statistics_Table2.R
144 lines (108 loc) · 4.61 KB
/
AQUACOSM2022_Statistics_Table2.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
144
#_______________________________________________________________________________
#
# AQUACOSM-Plus: Statistics for Speed
# 21 February 2022
# Author: Anika Happe
#
#_______________________________________________________________________________
#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)
library(rstatix) #for tukey_hsd()
#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 <- mutate(Data_Growth, growth_rate_2 = Growth_Lin+1)
Data_Growth <- filter(Data_Growth, Run == 1)
Data_Growth <- filter(Data_Growth, T_ID != "Constant_6")
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$Balanced <- as.factor(Data_Growth$Balanced)
Data_NP <- filter(Data_Growth, Rep == 1)
#Read data and prepare for usage for LRR
#(based on the calculated LRR in script 1) -> Take data from there
Data_LRR <- read_excel("Excel Files for Stats/SupOpTrons_LRR_Speed_2.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 2 GLM: GROWTH ~ TEMP * RATIO * SPEED ####
# 1) Explore your data
#Data_NP <- filter(Data_Growth, Rep == 1)
# 2) Check for outliers
Data_Growth <- Data_Growth[-11, ]
par(mfrow = c(1, 3)) #Three figures in one graph
boxplot(Data_Growth$Growth_Lin) #Does it look normally distributed?
hist(Data_Growth$Growth_Lin) #Does it look normally distributed?
dotchart(Data_Growth$Growth_Lin) #Points very far away are outliers
# 6) Building the model
#Temp + Run
#T*N Interaction with all levels
glm_r_speed <- glm(Growth_Lin ~ Temp * Balanced * T_Treat, family = gaussian, data = Data_Growth)
plot(glm_r_speed) #Check the assumptions of the model!
hist(resid(glm_r_speed)) #Does the data look normally distributed?
summary(glm_r_speed)
summary(aov(glm_r_speed))
#tukey_results <- tukey_hsd(glm_r_speed)
#print(tukey_results, n = Inf)
#_______________________________________________________________________________
#### TABLE 2 GLM: LRR ~ TEMP * RATIO ####
# 2) Check for outliers
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_speed <- glm(LRR ~ Temp * Balanced, family = gaussian, data = Data_LRR)
plot(glm_lrr_speed) #Check the assumptions of the model!
hist(resid(glm_lrr_speed)) #Does the data look normally distributed?
summary(aov(glm_lrr_speed))
#_______________________________________________________________________________
#### TABLE 2 GLM: N:P ~ TEMP * RATIO * SPEED ####
# 2) Check for outliers
Data_NP <- Data_NP[-1, ]
par(mfrow = c(1, 3)) #Three figures in one graph
boxplot(Data_NP$NP_Ratio) #Does it look normally distributed?
hist(Data_NP$NP_Ratio) #Does it look normally distributed?
dotchart(Data_NP$NP_Ratio) #Points very far away are outliers
# 6) Building the model
#Temp + Run
#T*N Interaction with all levels
glm_NP_speed <- glm(NP_Ratio ~ Temp * Balanced * T_Treat, family = gaussian, data = Data_NP)
plot(glm_NP_speed) #Check the assumptions of the model!
hist(resid(glm_NP_speed)) #Does the data look normally distributed?
summary(aov(glm_NP_speed))
#tukey_hsd(glm_NP_speed)
#_______________________________________________________________________________