Get stuff

library(tidyverse)
library(lme4)

Estimate measurement error

From 3x biological reps from Tenaillon 2011. Unfortunately these data were not published with the paper.

#make data frame
t_line=c(rep("b73",3),rep("lux",3)) #line names
t_g=c(5.87,5.99,6.02,9.21,9.15,8.86) #genome size
tenaillon_data=data.frame(t_line,t_g)

#model as only one effect -- line. here as random. fixed better?
t1 = lmer(tenaillon_data$t_g ~ 1|tenaillon_data$t_line)
t1_line = unlist(VarCorr(t1))[[1]]; 
t1_res=attr(VarCorr(t1), "sc")^2

Try with replicates from Bilinski

download.file("https://gist.githubusercontent.com/rossibarra/888b0eae4af606e82fffb6c902a8ca82/raw/a1813bf1f3a4ad7494c266367a7cdd5b7578f6e9/gistfile1.txt","~/Desktop/rep.txt")
bilinski_rep<-read.table("~/Desktop/rep.txt",header=T)
b1<-lmer(bilinski_rep$Gsize~1|bilinski_rep$Line)
b1_line = unlist(VarCorr(b1))[[1]]; 
b1_res=attr(VarCorr(b1), "sc")^2

Estimate \(h^2\)

Get data

download.file("https://gist.githubusercontent.com/rossibarra/1a70e69715188abb72532998af811986/raw/536b5b3a0e6d8cbdbf78a6af838b789a77635011/grote_data%2520bilinski","~/Desktop/bilinski_gsize.txt")
dat<-read.table("~/Desktop/bilinski_gsize.txt",header=T) %>% mutate(mom=Pop) %>% select(mom,genomesize) %>% filter(mom!="x")

Estimate \(h^2\) from half-sib families:

m1 = lmer(dat$genomesize ~ (1|dat$mom))
m1_line = unlist(VarCorr(m1))[[1]]; 
m1_res=attr(VarCorr(m1), "sc")^2
h2=4*(m1_line)/(4*m1_line+m1_res)
h2
## [1] 0.6811212

Incorporate error from Tenaillon to re-estimate \(h^2\)

m1_res_star<-m1_res-t1_res
h2_star=4*(m1_line)/(4*m1_line+m1_res_star)
h2_star
## [1] 0.8917228

Or from Bilinski

m1_res_star<-m1_res-b1_res
h2_star=4*(m1_line)/(4*m1_line+m1_res_star)
h2_star
## [1] 0.7116897