-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRelated_trees_quantile_regression.Rmd
187 lines (162 loc) · 6.75 KB
/
Related_trees_quantile_regression.Rmd
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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
---
title: "Related Trees"
---
Code used to create Figure 2
```{r}
library(readxl)
library(lme4)
library(quantreg)
library(ggplot2)
```
Import data
```{r}
related <- read.csv("~/Documents/Ash/related_trees.csv", sep=",")
pheno <- read.csv("~/Documents/Ash/mp_mastersheet.10.03.2022_tidy_gebv.csv")
```
```{r}
filtered_phenotype <- pheno[which(pheno$Type %in% c("Adult", "Juvenile") & !(is.na(pheno$GEBV))),]
replicates <- c("S31R", "S54R","S67R", "S76R", "S82R", "S223R", "S332Q", "S410Q", "S415R", "S571Q",
"S655R", "S868R", "S885R", "S895R", "S901Q", "S925Q", "S928Q", "S939Q", "S954R", "S987Q",
"SA", "SB", "SE", "SH", "SJ", "S394C")
filtered_phenotype <- filtered_phenotype[which(!(filtered_phenotype$Sample %in% replicates)),]
```
Count the number of offspring(sire or dam) per parent tree
```{r}
num_offspring <- table(c(related$dam, related$sire))
num_offspring_tab <- data.frame(num_offspring)
colnames(num_offspring_tab) <- c("Parent", "Num_Offspring")
```
```{r}
parents <- filtered_phenotype[which(filtered_phenotype$Label %in% names(num_offspring)),]
parent_health <- data.frame(parents$Label, parents$DBH, parents$PercentScore_2019, parents$PercentScore_2021, parents$GEBV)
colnames(parent_health) <- c("Parent", "DBH", "Percent_2019", "Percent_2021", "GEBV")
```
Plot trees individually and weight by the confidence of the parentage assignment
Weight by LLR/max(LLR)
```{r}
dam_health <- merge(related, parent_health, by.x = 'dam', by.y = 'Parent')
dam_health$LLR <- dam_health$LLRdam
dam_health$Parent_ID <- dam_health$dam
sire_health <- merge(related, parent_health, by.x = 'sire', by.y = 'Parent')
sire_health$LLR <- sire_health$LLRsire
sire_health$Parent_ID <- sire_health$sire
invd_parent_health <- rbind(dam_health, sire_health)
invd_parent_health <- merge(invd_parent_health, filtered_phenotype, by.x = 'id', by.y = 'Label')
head(invd_parent_health)
```
Quantile fit of juvenile score against adult GEBV
```{r}
# Fit quantile regression for the median (50th percentile)
model_50 <- rq(Score_2019 ~ GEBV.x, tau = 0.5, data = invd_parent_health)
# Fit for 25th and 75th percentiles
model_25 <- rq(Score_2019 ~ GEBV.x, tau = 0.25, data = invd_parent_health)
model_75 <- rq(Score_2019 ~ GEBV.x, tau = 0.75, data = invd_parent_health)
summary(model_25)
summary(model_50)
summary(model_75)
# Create a new data frame with GEBV values for prediction
new_data <- data.frame(GEBV.x = seq(min(invd_parent_health$GEBV.x), max(invd_parent_health$GEBV.x), length.out = 100))
# Predict Tree health for the different quantiles
new_data$Tree_health_25 <- predict(model_25, newdata = new_data)
new_data$Tree_health_50 <- predict(model_50, newdata = new_data)
new_data$Tree_health_75 <- predict(model_75, newdata = new_data)
parent_plot <- ggplot(invd_parent_health, aes(x = GEBV.x, y = Score_2019)) +
geom_point(color = "blue",
size = 2,
alpha = 0.6) + # Scatter plot of data points
geom_line(
data = new_data,
aes(x = GEBV.x, y = Tree_health_25),
color = "green",
linetype = "dashed",
size = 1
) + # 25th percentile
geom_line(
data = new_data,
aes(x = GEBV.x, y = Tree_health_50),
color = "red",
size = 1.2
) + # 50th percentile (median)
geom_line(
data = new_data,
aes(x = GEBV.x, y = Tree_health_75),
color = "orange",
linetype = "dashed",
size = 1
) + # 75th percentile
labs(x = "Parent GEBV", y = "2019 Health Score", title = "Quantile Regression of Juvenile Score vs Parent GEBV") +
theme_minimal()
parent_plot
```
juveniles
```{r}
juveniles <- filtered_phenotype[filtered_phenotype$Type=="Juvenile",]
adults <- filtered_phenotype[filtered_phenotype$Type=="Adult",]
```
```{r}
# Fit quantile regression for the median (50th percentile)
model_50 <- rq(Score_2019 ~ GEBV, tau = 0.5, data = juveniles)
# Fit for 25th and 75th percentiles
model_25 <- rq(Score_2019 ~ GEBV, tau = 0.25, data = juveniles)
model_75 <- rq(Score_2019 ~ GEBV, tau = 0.75, data = juveniles)
summary(model_25)
summary(model_50)
summary(model_75)
# Create a new data frame with GEBV values for prediction
new_data <- data.frame(GEBV = seq(min(juveniles$GEBV), max(juveniles$GEBV), length.out = 500))
# Predict Tree health for the different quantiles
new_data$Tree_health_25 <- predict(model_25, newdata = new_data)
new_data$Tree_health_50 <- predict(model_50, newdata = new_data)
new_data$Tree_health_75 <- predict(model_75, newdata = new_data)
juv_plot <- ggplot(juveniles, aes(x = GEBV, y = Score_2019)) +
geom_point(color = "blue",
size = 2,
alpha = 0.6) + # Scatter plot of data points
geom_line(
data = new_data,
aes(x = GEBV, y = Tree_health_25),
color = "green",
linetype = "dashed",
size = 1
) + # 25th percentile
geom_line(
data = new_data,
aes(x = GEBV, y = Tree_health_50),
color = "red",
size = 1.2
) + # 50th percentile (median)
geom_line(
data = new_data,
aes(x = GEBV, y = Tree_health_75),
color = "orange",
linetype = "dashed",
size = 1
) + # 75th percentile
labs(x = "GEBV", y = "2019 Health Score", title = "Quantile Regression of Juvenile Score vs GEBV") +
theme_minimal()
juv_plot
```
```{r}
# Fit quantile regression for the median (50th percentile)
model_50 <- rq(PercentScore_2019 ~ GEBV, tau = 0.5, data = adults)
# Fit for 25th and 75th percentiles
model_25 <- rq(PercentScore_2019 ~ GEBV, tau = 0.25, data = adults)
model_75 <- rq(PercentScore_2019 ~ GEBV, tau = 0.75, data = adults)
summary(model_25)
summary(model_50)
summary(model_75)
# Create a new data frame with GEBV values for prediction
new_data <- data.frame(GEBV = seq(min(adults$GEBV), max(adults$GEBV), length.out = 500))
# Predict Tree health for the different quantiles
new_data$Tree_health_25 <- predict(model_25, newdata = new_data)
new_data$Tree_health_50 <- predict(model_50, newdata = new_data)
new_data$Tree_health_75 <- predict(model_75, newdata = new_data)
adult_plot <- ggplot(adults, aes(x = GEBV, y = PercentScore_2019)) +
geom_point(color = "blue", size = 2, alpha = 0.6) + # Scatter plot of data points
geom_line(data = new_data, aes(x = GEBV, y = Tree_health_25), color = "green", linetype = "dashed", size = 1) + # 25th percentile
geom_line(data = new_data, aes(x = GEBV, y = Tree_health_50), color = "red", size = 1.2) + # 50th percentile (median)
geom_line(data = new_data, aes(x = GEBV, y = Tree_health_75), color = "orange", linetype = "dashed", size = 1) + # 75th percentile
labs(x = "GEBV", y = "2019 Canopy Cover", title = "Quantile Regression of Adult Score vs GEBV") +
theme_minimal()
adult_plot
```