---
title: "demoscript_new"
author: "Jan-Peter george"
date: "2026-05-28"
output: html_document
---

This is a short tutorial, which walks the user through the most crucial steps during data analysis. These steps include:

-Selection of subsets

-Handling dates

-Creating new columns

-Finding and eliminating duplicates

\- Filtering data and outlier removal

-Plotting data

-Building simple transfer models

First, load some R packages which are necessary to carry out the work...

```{r}
require("data.table")
require("ggplot2")
require("dplyr")
require("tidyverse")
require("lubridate")
require("stringr")
require("stringx")
require("breedR")
require("rnaturalearth")
require("car")
require("EnvStats")
```

Then import the Rdata sources for Trial sites, provenances, and tree observations. Make sure that you use your local paths here...

```{r}
tree_obs<-readRDS("your local path")
provs<-readRDS("your local path")
sites<-readRDS("your local path")
```

Lets select Quercus petraea as an example species and Height as a focal trait. This subset gives us 54,865 observations

```{r}
QP<-tree_obs[grepl("QUPE",OptFORESTS.trial.code)]
str(QP)
###transform into data.table
QP<-as.data.table(QP)
####Selecting trait of interest
unique(QP$Trait)
QP_height <- QP[Trait == "Height"]
str(QP_height)
```

Then, join the tree observations with trial site and provenance metadata to have it all in one data table

```{r}
QP_height_prov<-QP_height %>% left_join(provs,by=c("OptFORESTS.trial.code","OptFORESTS.provenance.code"))
QP_height_prov_site<- QP_height_prov %>% left_join(sites, by=("OptFORESTS.trial.code"))
head(QP_height_prov_site)
```

Lets bring the establishment year of the trial into a more handy data format...

```{r}
unique(QP_height_prov_site$Establishment.year)
QP_height_prov_site[, establishment.year.new := str_replace_all(Establishment.year,c("2007"="2007","2008"="2008","1992"="1992","1978"="1978","1990-1995" = "1990","1990-1996"="1990")) ]  
QP_height_prov_site[, establishment.year.new :=ymd(QP_height_prov_site$establishment.year.new, truncated = 2L)]
unique(QP_height_prov_site$establishment.year.new)
```

With the wrangled dates we can create a new column, which we call "age_at_assessment" and plot the different ages in a histogram

```{r}
QP_height_prov_site[,age_at_assessment:=as.numeric(Assessment.date-establishment.year.new)/365]
hist(QP_height_prov_site$age_at_assessment, xlab= "Age")
```

We select all ages between 20 & 30

```{r}
QP_height_prov_site<-QP_height_prov_site[age_at_assessment>=20 & age_at_assessment<=30]
```

Which observations are repeated in the subset?

```{r}
summary<-QP_height_prov_site[,.(Dates=list(unique(Assessment.date))),by=OptFORESTS.trial.code]

summary
```

Eliminating repeated observations

```{r}
QP_final<-setDT(QP_height_prov_site)[order(Assessment.date),.SD[.N], by = .(OptFORESTS.tree.id)]
QP_final[,.(Dates=list(unique(Assessment.date))),by=OptFORESTS.trial.code]
QP_final[OptFORESTS.trial.code=="QUPE-DK-01"][1]
QP_final[OptFORESTS.trial.code=="QUPE-DK-01"][4257]
```

Remove double-trees with ID pattern

```{r}
QP_final_filt<-QP_final %>%
  filter(!grepl("_DK", OptFORESTS.tree.id))
QP_final_filt[,.(Dates=list(unique(Assessment.date))),by=OptFORESTS.trial.code]
```

Check if all duplicates have been removed

```{r}
which(duplicated(QP_final_filt$OptFORESTS.tree.id))
```

Plot the data in a histogram to see how the data is distributed

```{r}
hist(QP_final_filt$Value,
     xlab = "Height",
     main = "Height",
     breaks = sqrt(length(QP_final_filt$Value))) # set number of bins
```

Another histogram to see what data types have been used. Obviously, not all data have been measured, but some data have been imputed from dbh-height growth curves...

```{r}
ggplot(QP_final_filt, aes(x = Value, fill = Data.type)) + 
  geom_histogram()
```

Applying a Hampel filter:

The Hampel Filter identifies outliers based on the median absolute deviation (MAD), a measure less affected by outliers in the data than the standard deviation

```{r}
lower_bound <- median(na.omit(QP_final_filt$Value)) - 3 * mad(na.omit(QP_final_filt$Value), constant = 1)
lower_bound

upper_bound <- median(na.omit(QP_final_filt$Value)) + 3 * mad(na.omit(QP_final_filt$Value), constant = 1)
upper_bound
outlier_ind <- which(QP_final_filt$Value < lower_bound | QP_final_filt$Value > upper_bound)
outlier_ind
```

Test if there are significant outliers:

```{r}
qqPlot(QP_final_filt$Value)
outlier.test <- rosnerTest(QP_final_filt$Value, k=3)
outlier.test
outlier.test$all.stats
```

Create an overview map to see which provenances have been tested at which trial site. Obviously, the romanian provs have been tested only in the romanian trial site.

```{r}
world <- ne_countries(scale = 50, returnclass = 'sf')

world %>% 
  ggplot() +
  geom_sf() +
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65))+
  geom_point(data = QP_final_filt, size=2,position = position_jitterdodge(dodge.width = 1, jitter.width = 0.1),aes(x = Longitude.x, y = Latitude.x,col=OptFORESTS.trial.code))+
  geom_point(data = QP_final_filt,size=2, aes(x = Longitude.y, y = Latitude.y),color="black")
```

Data Analysis part: Building a a simple mixed model with unstructured random effects

```{r}
str(QP_final_filt)
unique(QP_final_filt$Genetic.background)
```

Building a remlf90 model with provenance and block as random effect and intercept as fixed effect

```{r}
model_1 <- remlf90(
  fixed = Value ~ 1,
  random=~ OptFORESTS.provenance.code + Block.no ,
  data = QP_final_filt)
```

Look at the model summary to see if genetics play a significant role:

```{r}
summary(model_1)
```

Add the trial site as an additional random effect:

```{r}
model_2 <- remlf90(
  fixed = Value ~ 1,
  random=~ OptFORESTS.provenance.code + OptFORESTS.trial.code + Block.no ,
  data = QP_final_filt
)
summary(model_2)
```

Account for different ages in the model:

```{r}
model_3 <- remlf90(
  fixed = Value ~ 1,
  random=~ OptFORESTS.provenance.code + Block.no + age_at_assessment ,
  data = QP_final_filt
)
summary(model_3)
```

Are there signatures of local adaptation in the data?

First, create a new variable called "distance_from_home":

```{r}
QP_final_filt[,distance_from_home:=Latitude.x-Latitude.y]
```

Create a basic plot with Height as a function of transfer distance:

```{r}
ggplot(QP_final_filt,aes(x=distance_from_home,y=Value)) + 
  geom_boxplot(aes(group = OptFORESTS.provenance.code),size=0.3)+
  geom_dotplot(aes(group = OptFORESTS.provenance.code),
    binaxis='y', stackdir='center',
    color = "black", fill = "#999999",dotsize = 0.1,
    position = position_dodge(0.8)
  )+
  geom_jitter(
    position = position_jitter(0.2),
    size = 1.2
  )+
  stat_smooth(aes(y = Value),method = "lm", formula = y ~ x + I(x^2), size = 1)
```

Wrap the plot to include age as a co-variate:

```{r}
ggplot(QP_final_filt,aes(x=distance_from_home,y=Value)) + 
  geom_boxplot(aes(group = OptFORESTS.provenance.code),size=0.3)+
  geom_dotplot(aes(group = OptFORESTS.provenance.code),
    binaxis='y', stackdir='center',
    color = "black", fill = "#999999",dotsize = 0.1,
    position = position_dodge(0.8)
  )+
  geom_jitter(
    position = position_jitter(0.2),
    size = 1.2
  )+
  stat_smooth(aes(y = Value),method = "lm", formula = y ~ x + I(x^2), size = 1)+
  facet_wrap(vars(age_at_assessment))
```

Include age as a random factor in the model:

```{r}
model_4 <- remlf90(
  fixed = Value ~ 1,
  random=~ distance_from_home + Block.no + age_at_assessment ,
  data = QP_final_filt
)
summary(model_4)
```

Part 2: Survival data

Load the data, make sure to use your local path!

```{r}
surv_2.0<-readRDS("your local path")
str(surv_2.0)
unique(surv_2.0$data.type)
```

Select Laric decidua as an example...

```{r}
LD<-surv_2.0[grepl("LADE",OptFORESTS.trial.code)]
```

Import provenance and trial site metadata and select a subset:

```{r}
LD_prov<-LD %>% left_join(provs,by=c("OptFORESTS.trial.code","OptFORESTS.provenance.code"))
LD_prov_site<- LD_prov %>% left_join(sites, by=("OptFORESTS.trial.code"))
LD_prov_site<-LD_prov_site[!is.na(amr.corr)]
LD_prov_site<-LD_prov_site[amr.corr>0 & amr.corr<0.9]
```

Lets have a look how the data is distributed

```{r}
ggplot(LD_prov_site, aes(x = amr.corr, fill = data.type)) + 
  geom_histogram()
```

Building another reml model with provenance, block, and trial site as random effects

```{r}
model_5 <- remlf90(
  fixed =amr.raw ~ 1,
  random=~ OptFORESTS.provenance.code + Block.no + OptFORESTS.trial.code,
  data = LD_prov_site
)
summary(model_5)
```