Para o seguinte dataset contendo avaliações de severidade de mancha alvo (Corynespora casiicola) e rendimento ao nivel de parcelas, de 26 ensaios de eficiencia de controle de fungicidas, serão ajustados regressões lineares para obter finalmente uma tabela com os parámetros b0 (intercepto) e b1 (coeficiente angular), e os respetivos desvios padrões.

Este resultante data frame será analizado por medio de meta-analises para estimar uma reta única para os 26 ensaios, representando a perda de rendimento da soja devido à mancha alvo.

##   year local trt rep sev   prod
## 1 2012     1   1   1  45 2642.3
## 2 2012     1   1   2  40 2864.8
## 3 2012     1   1   3  45 2306.6
## 4 2012     1   1   4  40 2746.9
## 5 2012     1   2   1  30 2593.7
## 6 2012     1   2   2  30 2718.7
##      year local trt rep sev   prod
## 1028 2014    16   8   3 7.4 2697.2
## 1029 2014    16   8   4 7.2 2756.0
## 1030 2014    16   9   1 5.2 2761.4
## 1031 2014    16   9   2 4.8 3221.0
## 1032 2014    16   9   3 6.0 3030.8
## 1033 2014    16   9   4 5.4 3282.2

1- Criar ID para c/ ensaio

da$id <- as.factor(paste(da$year,da$local,sep="_")); head(da)
##   year local trt rep sev   prod     id
## 1 2012     1   1   1  45 2642.3 2012_1
## 2 2012     1   1   2  40 2864.8 2012_1
## 3 2012     1   1   3  45 2306.6 2012_1
## 4 2012     1   1   4  40 2746.9 2012_1
## 5 2012     1   2   1  30 2593.7 2012_1
## 6 2012     1   2   2  30 2718.7 2012_1

2- Estimar a media de cada ensaio * trt

tr_mean=aggregate(sev ~ id + trt, data = da, mean); tr_mean[1:27,] 
##         id trt    sev
## 1   2012_1   1 42.500
## 2   2012_2   1 36.250
## 3   2012_3   1 41.250
## 4   2012_4   1 30.625
## 5   2012_5   1  8.250
## 6   2012_6   1 27.000
## 7   2012_7   1 44.750
## 8   2012_8   1 47.750
## 9   2013_1   1 40.000
## 10 2013_10   1 26.500
## 11 2013_12   1 22.250
## 12  2013_3   1 48.500
## 13  2013_5   1 41.125
## 14  2013_8   1 45.675
## 15  2013_9   1 29.500
## 16  2014_1   1 28.125
## 17 2014_10   1 41.050
## 18 2014_11   1 11.000
## 19 2014_12   1 18.000
## 20 2014_14   1 24.750
## 21 2014_15   1 12.375
## 22 2014_16   1 29.350
## 23  2014_2   1 51.250
## 24  2014_5   1 18.500
## 25  2014_6   1 33.250
## 26  2014_8   1 12.475
## 27  2014_9   1 20.100

3- Calcular quantas filas têm cada ensaio e extrair o vector com filas/ensaio

nrows_exp=lapply(unique(da$id), function(x){ sum(da$id==x) })

a <- c()
for(i in 1:length(nrows_exp)) {a <- c(a,nrows_exp[[i]])}  
a
##  [1] 36 36 36 36 36 45 36 36 44 44 44 44 44 44 44 32 36 36 36 36 36 36 36 36 36 36 36

4- Adicionar coluna com media do tratamento controle

da$ctrl_sev=rep(tr_mean[1:length(levels(da$id)),3], t=a) 

# Criar coluna disease press com os niveis baixo (low: 10-30%) ou alto (high > 30%)

da$dis_pres <- ifelse(da$ctrl_sev >= 30,"h", "l"); head(da)
##   year local trt rep sev   prod     id ctrl_sev dis_pres
## 1 2012     1   1   1  45 2642.3 2012_1     42.5        h
## 2 2012     1   1   2  40 2864.8 2012_1     42.5        h
## 3 2012     1   1   3  45 2306.6 2012_1     42.5        h
## 4 2012     1   1   4  40 2746.9 2012_1     42.5        h
## 5 2012     1   2   1  30 2593.7 2012_1     42.5        h
## 6 2012     1   2   2  30 2718.7 2012_1     42.5        h

Reordenar dataset

dat= da[,c(7,5,6,8,9)]; head(dat)
##       id sev   prod ctrl_sev dis_pres
## 1 2012_1  45 2642.3     42.5        h
## 2 2012_1  40 2864.8     42.5        h
## 3 2012_1  45 2306.6     42.5        h
## 4 2012_1  40 2746.9     42.5        h
## 5 2012_1  30 2593.7     42.5        h
## 6 2012_1  30 2718.7     42.5        h

5- Descartar ensaios com severidade no controle < 10%

dat=subset(dat, ctrl_sev >10)

6- Separar o dataframe por ensaio! Pronto para estimar as regressões…

dat_id <- split(dat, f=dat$id)

Gráficos exploratórios

theme_set(theme_bw())
explot = ggplot(dat, aes(sev, prod)) + 
  geom_point(shape=20, size=2) + 
  facet_wrap(~id, ncol=4, scales="fixed") + coord_cartesian(xlim = c(0, 60)) +
  scale_x_continuous("Disease severity (%)") +
  scale_y_continuous("Yield (kg/ha)") +
  geom_smooth(method="lm", color="red")
explot

Regressão linear por ensaio

7- Criar função para extair os coeficientes e desvio padrões

myFun <-  function(lm)
  {
    out <- c(lm$coefficients[1],
             lm$coefficients[2],
             summary(lm)$coefficients[2,2],
             summary(lm)$coefficients[1,2]
    )
    names(out) <- c("b0","b1","b0.SE","b1.SE") 
    return(out)}

8- Calcular a regressão linear por ensaio

a <- list()
for (i in 1:length(dat_id)) {
  a[[names(dat_id[i])]] <- lm(dat_id[[i]]$prod ~ dat_id[[i]]$sev)
}

9- Criar um data frame com os coeficientes e desvios estimados

results <- list()
for (i in 1:length(a)) results[[names(a)[i]]] <- myFun(a[[i]])
b <- as.data.frame(results)
b=t(b)
id<-sapply(strsplit(rownames(b), "X"), "[[", 2)
reglin_coef = as.data.frame(b, id)
reglin_coef$id=row.names(reglin_coef)
rownames(reglin_coef) <- NULL
reglin_coef=reglin_coef[,c(5,1,2,3,4)]

10- Salvar o dataframe em un arquivo .txt para utilizar na metaregressão

da1=subset(da, trt=="1" & rep=="1")[,6:9]
dat_coef= as.data.frame(merge(da1,reglin_coef, "id")); dat_coef
##         id   prod ctrl_sev dis_pres       b0         b1     b0.SE     b1.SE
## 1   2012_1 2642.3   42.500        h 2898.667  -7.585391  7.085634 252.10767
## 2   2012_2 3135.0   36.250        h 3927.859 -24.477332  6.810604 201.26204
## 3   2012_3 2492.7   41.250        h 3251.704 -16.110989  7.418998 227.57640
## 4   2012_4 2957.8   30.625        h 3445.446  -4.364499  3.941473  86.87648
## 5   2012_6 3744.0   27.000        l 4218.295  -9.249662  7.425488 129.15512
## 6   2012_7 2992.5   44.750        h 3610.327  -7.228888  2.185961  68.90476
## 7   2012_8 3400.0   47.750        h 3848.940 -12.178492  6.589507 223.06366
## 8   2013_1 3493.4   40.000        h 3882.550  -4.588254  4.537889 123.84894
## 9  2013_10 2550.1   45.675        h 3330.020 -36.060893  5.574162  57.07209
## 10 2013_12 1208.4   29.500        l 1969.408 -13.027761  6.169600  72.30311
## 11  2013_3 2499.7   26.500        l 3751.011 -23.093483  2.966246 100.48090
## 12  2013_5 2833.9   22.250        l 4249.448 -35.079457  5.220999 149.42064
## 13  2013_8 2728.7   48.500        h 3617.205 -17.593282  2.909147  76.17176
## 14  2013_9 1284.3   41.125        h 1603.815 -15.716396  5.834118  81.44291
## 15  2014_1 3633.9   28.125        l 3883.763 -10.744523  3.000276  46.69593
## 16 2014_10 1770.0   29.350        l 4118.405 -60.802565  8.080833 269.20942
## 17 2014_11 2151.7   51.250        h 2740.186 -34.516343 11.252236  90.36883
## 18 2014_12 1486.0   18.500        l 2428.840 -49.601622  9.013081  92.18541
## 19 2014_14 2084.0   33.250        h 3670.118 -45.800428 11.329478 184.75890
## 20 2014_15 3034.9   12.475        l 4281.515 -80.603137 25.597768 143.58791
## 21 2014_16 2205.8   20.100        l 3005.912 -28.396429  5.929205  92.13784
## 22  2014_2 3345.1   41.050        h 3541.663  -5.668744  2.877861  93.80435
## 23  2014_5 3102.9   11.000        l 3920.497 -52.275505 14.465461 156.75904
## 24  2014_6 3819.4   18.000        l 3828.731  -6.238477  4.140976  75.97199
## 25  2014_8 3790.0   24.750        l 4559.966 -55.480689 11.510726  84.68645
## 26  2014_9 3750.0   12.375        l 4584.237 -49.166246  7.341558  74.38404
#write.table(dat_coef, file="dat_coef.txt",dec=",",sep="\t" )

Continua com a Metaregressão