-
Notifications
You must be signed in to change notification settings - Fork 0
/
MIxed Modeling
130 lines (113 loc) · 4.74 KB
/
MIxed Modeling
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
## Mixed Model Statistical procedures
## Developed by Max Gebhart
setwd("C:/Users/maxge/Downloads/Grad School Papers/Data-Analysis")
library(ggplot2)
library(dplyr)
library(tidyr)
library(gtools)
library(dismo)
library(raster)
library(sp)
library(lme4)
library(Matrix)
library(pspearman)
# First we will run a simple spearman correlation test to find what environmental variables
## we can expect within the mix model
x <- read.csv("Field Starry Data.csv", header=T)
head(x)
attach(x)
y<- cor.test(Bulbil, Temp, method="spearman", exact=F)
## This will make a blank matrix for the p values that are generated
pvalues <- matrix(nrow = 5, ncol=5, dimnames = list(c("Temp", "DO", "pH", "% Transmittance", "Submersed"), c("Density", "Bulbil", "Above", "Rhizoid", "Total")))
rvalues <- matrix(nrow = 5, ncol=5, dimnames = list(c("Temp", "DO", "pH", "% Transmittance", "Submersed"), c("Density", "Bulbil", "Above", "Rhizoid", "Total")))
{
##Density Run
{
a <- cor.test(Density, Temp, method="spearman", exact=F)
b <- cor.test(Density, DO, method="spearman", exact=F)
c <- cor.test(Density, pH, method="spearman", exact=F)
d <- cor.test(Density, X..Transmittance, method="spearman", exact=F)
e <-cor.test(Density, Submersed.Light, method="spearman", exact=F)
rvalues[1] <- a$estimate
rvalues[2] <- b$estimate
rvalues[3] <- c$estimate
rvalues[4] <- d$estimate
rvalues[5] <- e$estimate
}
## Bulbil Run
{
a <- cor.test(Bulbil, Temp, method="spearman", exact=F)
b <- cor.test(Bulbil, DO, method="spearman", exact=F)
c <- cor.test(Bulbil, pH, method="spearman", exact=F)
d <- cor.test(Bulbil, X..Transmittance, method="spearman", exact=F)
e <-cor.test(Bulbil, Submersed.Light, method="spearman", exact=F)
rvalues[6] <- a$estimate
rvalues[7] <- b$estimate
rvalues[8] <- c$estimate
rvalues[9] <- d$estimate
rvalues[10] <- e$estimate
}
## Above run
{
a <- cor.test(Above, Temp, method="spearman", exact=F)
b <- cor.test(Above, DO, method="spearman", exact=F)
c <- cor.test(Above, pH, method="spearman", exact=F)
d <- cor.test(Above, X..Transmittance, method="spearman", exact=F)
e <-cor.test(Above, Submersed.Light, method="spearman", exact=F)
rvalues[11] <- a$estimate
rvalues[12] <- b$estimate
rvalues[13] <- c$estimate
rvalues[14] <- d$estimate
rvalues[15] <- e$estimate
}
## Rhizoid run
{
a <- cor.test(Rhizoid, Temp, method="spearman", exact=F)
b <- cor.test(Rhizoid, DO, method="spearman", exact=F)
c <- cor.test(Rhizoid, pH, method="spearman", exact=F)
d <- cor.test(Rhizoid, X..Transmittance, method="spearman", exact=F)
e <-cor.test(Rhizoid, Submersed.Light, method="spearman", exact=F)
rvalues[16] <- a$estimate
rvalues[17] <- b$estimate
rvalues[18] <- c$estimate
rvalues[19] <- d$estimate
rvalues[20] <- e$estimate
}
## Total run
{
a <- cor.test(Total, Temp, method="spearman", exact=F)
b <- cor.test(Total, DO, method="spearman", exact=F)
c <- cor.test(Total, pH, method="spearman", exact=F)
d <- cor.test(Total, X..Transmittance, method="spearman", exact=F)
e <-cor.test(Total, Submersed.Light, method="spearman", exact=F)
rvalues[21] <- a$estimate
rvalues[22] <- b$estimate
rvalues[23] <- c$estimate
rvalues[24] <- d$estimate
rvalues[25] <- e$estimate
}
write.csv(rvalues, file = "Starry (R value) Correlations.csv")
}
#lme4 has built in data from an Arabidopsis study so we can use this as test data
## test data is important for establishing that the procedures are working
## this is also a small crash course into the syntax for this package
{
## this is to get a simple linear model from the data looking at the effects
## between two variables and the summary will give all values
Arabidopsis <- Arabidopsis
ar <- lm(gen ~ total.fruits, data = Arabidopsis)
summary(ar)
# this will be to get a mixed model using the syntax with rack as a random effect
## and total.fruits as a fixed effect you can add more to both following the next lines
ara <- lmer(gen ~ total.fruits + (1 | rack), data = Arabidopsis)
## adding further variables using this (1 | x) will add further to random effects
## x represents your random effects variable
ara <- lmer(gen ~ total.fruits + (1 | rack) + (1 | nutrient), data = Arabidopsis)
## by doing the syntax x ~ y + z will add more to your fixed effects list
## x, y, and z all being variables you want to understand
ara <- lmer(gen ~ total.fruits + nutrient + (1| rack), data = Arabidopsis)
## unfortunately you'll need to keep your variables similar as numbers
## I suggest using binary values for qualitative variables if possible
## the model created by the lmer function cannot be saved as a csv so you might need to save certain parts
## for a csv file if needed
}