Canola Black stem (Leptosphaeria maculans / Phoma macdonaldii)
Experimento 1
El objetivo de este trabajo fue determinar aquellos estados fenológicos de colza en que aparición de manchas foliares de Phoma implicarÃan una máxima severidad de cancro en estados reproductivos avanzados.
Recicle este grafico porque no encuentro ni el archivo de datos original ni el script. Me parece que no habia mucho para agregar. En caso de necesidad de modificaciones podria hacer el script nuevamente.
Experimento 2
El objetivo de este experimento fue relacionar la epidemia inicial de manchas foliares con el cancro en la base del tallo en estados reproductivos avanzados.
Datos de epidemia de manchas foliares a través del tiempo.
head(datmf, 20)## Grouped Data: inc ~ tt | bq/par
## par tr bq tt inc
## 1 101 1 1 15 3.00
## 2 101 1 1 58 30.00
## 3 101 1 1 95 43.00
## 4 101 1 1 146 46.00
## 5 101 1 1 165 44.00
## 6 201 1 2 15 3.00
## 7 201 1 2 58 29.00
## 8 201 1 2 95 38.00
## 9 201 1 2 146 43.00
## 10 201 1 2 165 40.00
## 11 301 1 3 15 3.00
## 12 301 1 3 58 31.00
## 13 301 1 3 95 40.00
## 14 301 1 3 146 48.00
## 15 301 1 3 165 46.00
## 16 102 2 1 15 0.90
## 17 102 2 1 58 10.00
## 18 102 2 1 95 20.37
## 19 102 2 1 146 18.52
## 20 102 2 1 165 15.74Exploración gráfica
tiempot <- "Thermal time (°c.d)"
ls_inc <- "Phoma leaf spot incidence (%)"
xyplot(inc~tt|par, data=datmf,layout=c(11,3), type=c("g","l","p"),
xlab=tiempot, ylab=ls_inc,
panel=function(x,y,pn=panel.number(),...){
panel.xyplot(x,y,...)}
)
Ajustes de modelos:
Modelos con efectos fijos atribuyen un peso a cada unidad experimental siendo este el inverso de la varianza (1/v). Es decir que la heterogeneidad entre unidades experimentales es aleatoria.
Modelos con efectos mixtos (fijos + aleatorios) atribuyen un peso a cada unidade experimental: inverso de la varianza mas la heterogeneidad (1/v + h). Es decir que incorporan variabilidad entre unidades experimentales, en nuestro caso, las parcelas. De este modo, modelos mixtos generan resultados con mayor intervalo de confianza. Este tipo de modelos, en general, serian apropiados para experimentos de campo con presencia de heterogenidad espacial, debido a multiples factores ambientales, como propiedades de suelo, relieve, ambientales, etc.
La justificativa que encuentro de porque usamos este tipo de abordaje, a diferencia de simplemente registrar la incidencia de manchas foliares máxima es porque particularmente en este patosistema las temperaturas óptimas de crecimiento son diferentes para el hongo y para la planta. Esta ultima continua emitiendo nuevas hojas en bajas temperaturas en que el hongo deja de diseminarse y entra en estado asintomático. Por eso se puede ver una caida en la incidencia, cosa que no es real. Este tipo de modelos puede continuar prediciendo valores máximos un poco más a pesar de la baja en la actividad diseminadora.
En esta primera instancia intentaremos lograr convergencia global (todos los tratamientos juntos) de los siguientes modelos.
LogÃstico
The logistic model assumes that the rate of change of the disease is proportional both to the amount of disease and to the proportion of plant material not affected by the disease:
\[dY/dt = r_l Y (1-Y)\]
\[Y = A / (1 + exp ( B - C * x))\]
Monomolecular
The mono-molecular (negative exponential) model assumes that the rate of change of the disease is proportional to the proportion of plant material not affected by the disease:
\[dY/dt = r_m (1-Y)\]
\[Y = A ∗ (1 − B ∗ (exp −C x)\]
Gompertz
The Gompertz model assumes that the rate of change of the disease is proportional both to the amount of disease and to the log of the amount not affected by the disease.
\[dY/dt = r_g Y [-log Y]\]
\[Y = A ∗ exp(-exp (B − (C ∗ x))\]
Primer intento de ajuste
Modelo Gompertz con 2 parametros (inoculo inicial lo considero fijo, usando el valor medio de la primera observaión de sÃntomas). En el modelo a ajustar se observa:
- Efectos fijos: tratamientos (asÃntota y tasa)
- Efecto aleatorio: bloque y parcela (asÃntota)
\[y_{ijkl} = (A_i + a_{j/k}) ∗ exp(-exp (B_i− (C_i ∗ x_l) )\]
y: es la incidencia de manchas foliares (nivel de parcela) con el i-ésimo tratamiento de fungicida, j-esimo bloque, en la k-esima parcela, en el l-esimo tiempo.
A: es el efecto fijo del tratamiento en la asintota. a: efecto aleatorio, en la k-esima parcela anidada en el j-esimo bloque.
B: efecto fijo del tratamiento en el intercepto, o sea la incidencia de phoma cuando el tiempo es x=0. En nuestro caso no tenemos interes en estimar su valor ya que las parcelas tenian la misma presion de inoculo, el cual fue esparcido homogeneamente previo a la siembra. Por ello, será fijado en un valor aproximado a la media de incidencia al momento de primera observacion de sintomas.
C: efecto fijo del tratamiento de fungicida en la tasa de progreso de la enfermedad.
starters = c(
A=c(40,-27,-17,-14,-6,-4.5,-27,-12.5,-3.5,3,3.5),
C=c(0.026,-0.006,-0.016,0.006,-0.002,0.005,-0.019,0.009,0.028,0.005,0.003))
gomp1<- nlme(inc~ A*(exp(-exp(2-C*tt))),fixed=A+C~tr, random=A~1,
data=datmf, start=starters)
p1 = plot(gomp1, cex=0.7)
p2 = qqnorm(gomp1, ~ resid(., type = "p"), abline = c(0, 1), , cex=0.7)
grid.arrange(p1, p2, ncol = 2)
Se observa heterocedasiticidad » la varianza aumenta con el tiempo térmico
- Incorporar una varianza que depende el valor ajustado
gomp2 <- update(gomp1, weights=varPower(form=~fitted(.)))
p3 = plot(gomp2, cex=0.7)
p4 = qqnorm(gomp2, ~ resid(., type = "p"), abline = c(0, 1), cex=0.7)
grid.arrange(p3, p4, ncol = 2)
plot(augPred(gomp2), xlab=tiempot, ylab=ls_inc, layout=c(11,3))
- Comparacion de modelos
anova(gomp1,gomp2) ## Model df AIC BIC logLik Test L.Ratio
## gomp1 1 25 987.1184 1064.7671 -468.5592
## gomp2 2 26 916.8698 997.6243 -432.4349 1 vs 2 72.24866
## p-value
## gomp1
## gomp2 <.0001Incremento em logLik, menor BIC y AIC » la modificación fue significativo para estimar los parámetros.
- Ajustar modelo monomolecular
\[Y_{ijkl} = A_i + a_{j/k} ∗ (1 − B ∗ (exp −C_i ∗ x_l )+ ε_{ijkl}\]
mono1<-nlme(inc~ A*(1-2*exp(-C*tt)),fixed=A+C~tr, random=A~1, data=datmf,start=starters,
weights=varPower(form=~fitted(.)))
p5 = plot(mono1, cex=0.7)
p6 = qqnorm(mono1, ~ resid(., type = "p"), abline = c(0, 1), cex=0.7)
grid.arrange(p5, p6, ncol = 2)
- Comparar el mejor Gompertz con el mejor Monomolecular
AIC(gomp2, mono1)## df AIC
## gomp2 26 916.8698
## mono1 26 922.6920Gompertz=gomp2; Monomolecular=mono1Continua siendo mejor el gompertz (menor AIC)
plot(comparePred(Gompertz, Monomolecular), xlab=tiempot, ylab=ls_inc, layout=c(11,3))
- Probar Gompertz con A + C random
A priori no seria correcto ya que incurririamos en sobreparametrización del modelo considerando las pocas evaluaciones realizadas y la gran cantidad (4) de parámetros a estimar (2 fijos y 2 aleatorios).
\[y_{ijkl} = (A i + a_{j/k}) ∗ exp(-exp(B_i − (C_i + c_{j/k} ∗ x_l) )\]
gomp2AC<-update(gomp2,random=A + C ~1)
print(gomp2AC, correlation=F)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: inc ~ A * (exp(-exp(2 - C * tt)))
## Data: datmf
## Log-likelihood: -439.2405
## Fixed: A + C ~ tr
## A.(Intercept) A.tr2 A.tr3 A.tr4
## 4.329568e+01 -2.658452e+01 -2.026196e+01 -1.433822e+01
## A.tr5 A.tr6 A.tr7 A.tr8
## -5.254011e+00 -1.791213e+00 3.832322e+03 -7.803276e+00
## A.tr9 A.tr10 A.tr11 C.(Intercept)
## -1.681995e+00 4.494090e+00 2.806137e+00 5.562050e-02
## C.tr2 C.tr3 C.tr4 C.tr5
## -1.003911e-02 -2.101188e-02 1.076329e-02 -4.939191e-03
## C.tr6 C.tr7 C.tr8 C.tr9
## -1.616778e-03 -5.402127e-02 -3.298786e-03 1.431406e-02
## C.tr10 C.tr11
## 6.055060e-04 4.204952e-03
##
## Random effects:
## Formula: list(A ~ 1, C ~ 1)
## Level: bq
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## A.(Intercept) 1.105156810 A.(In)
## C.(Intercept) 0.002988929 0.999
##
## Formula: list(A ~ 1, C ~ 1)
## Level: par %in% bq
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## A.(Intercept) 1.491548e-06 A.(In)
## C.(Intercept) 7.392477e-03 -0.98
## Residual 8.817813e-01
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 0.4458256
## Number of Observations: 165
## Number of Groups:
## bq par %in% bq
## 3 33p7 = plot(gomp2AC, cex=0.7)
p8 = qqnorm(gomp2AC, ~ resid(., type = "p"), abline = c(0, 1), cex=0.7)
grid.arrange(p7, p8, ncol = 2)
plot(comparePred(gomp2, gomp2AC), xlab=tiempot, ylab=ls_inc, layout=c(11,3))
AIC(gomp2, gomp2AC)## df AIC
## gomp2 26 916.8698
## gomp2AC 30 938.4810anova(gomp2,gomp2AC) # Mejor el gomp2! sobreparametrizado el gomp2AC## Model df AIC BIC logLik Test
## gomp2 1 26 916.8698 997.6243 -432.4349
## gomp2AC 2 30 938.4810 1031.6594 -439.2405 1 vs 2
## L.Ratio p-value
## gomp2
## gomp2AC 13.61127 0.0086Se puede observar mejor acople de la curva ajustada a los valores observados pero el costo de incrementar un parámetro se vió reflejado en el AIC, criterio por el cual continuamos con el modelo mas parcimonioso.
Con esta evidencia continuamos con el modelo de Gompertz de 3 parámetros: Tasa de crecimiento (C) = fija AsÃntota (A) = fijo y aleatorio
Vale mencionar que el modelo logistico no logró convergencia global
Modelo para extraer los coeficientes de cada tratamiento
(variación: fixed=A+B+C~0+tr)
sta=c(A=c(40,40-27,40-17,40-14,40-6,40-4.5,40-27,40-12.5,40-3.5,40+3,40+3.5),
C=c(0.026,0.026-0.006,0.026-0.016,0.026+0.006,0.026-0.002,0.026+0.005,0.026-0.019,0.026+0.009,
0.026+0.028,0.026+0.005,0.026+0.003))
gompi <- nlme(inc~A*(exp(-exp(1.25-C*tt))),fixed=A+C~0+tr, random=A~1, data=datmf,
start = sta,
weights=varPower(form=~fitted(.))) Efectos fijos
summary(gompi)$tTable## Value Std.Error DF t-value p-value
## A.tr1 45.94900610 3.232383918 111 14.215207 9.515800e-27
## A.tr2 17.12143629 2.280425055 111 7.508002 1.616775e-11
## A.tr3 24.97611930 4.031195929 111 6.195710 1.009183e-08
## A.tr4 30.18312106 2.450097872 111 12.319149 1.676825e-22
## A.tr5 39.62887694 3.119145740 111 12.705042 2.234439e-23
## A.tr6 41.59437355 2.930873562 111 14.191801 1.071365e-26
## A.tr7 12.28029709 2.242334771 111 5.476567 2.729972e-07
## A.tr8 35.13487477 2.738519012 111 12.829882 1.167034e-23
## A.tr9 42.84773983 2.724899841 111 15.724519 5.182010e-30
## A.tr10 48.82822498 3.133132853 111 15.584473 1.029192e-29
## A.tr11 48.40773283 3.146579463 111 15.384240 2.756362e-29
## C.tr1 0.03316331 0.004265814 111 7.774203 4.172053e-12
## C.tr2 0.03017098 0.006340011 111 4.758822 5.907718e-06
## C.tr3 0.02122373 0.004518197 111 4.697389 7.598033e-06
## C.tr4 0.03920182 0.006304091 111 6.218473 9.063752e-09
## C.tr5 0.03183108 0.004411070 111 7.216179 7.023033e-11
## C.tr6 0.03604453 0.004890987 111 7.369583 3.252289e-11
## C.tr7 0.02611450 0.006701292 111 3.896935 1.670139e-04
## C.tr8 0.03578094 0.005272432 111 6.786419 5.893066e-10
## C.tr9 0.04246875 0.005828743 111 7.286092 4.948023e-11
## C.tr10 0.03624852 0.004547322 111 7.971401 1.517223e-12
## C.tr11 0.03576988 0.004500357 111 7.948232 1.709247e-12cigompi <- intervals(gompi)$fixed
cigompi <- cigompi[,c(2,1,3)]
colnames(cigompi) <- c("est","lwr","upr")
cigomp <- as.data.frame(rbind(cigompi))
cigomp$par <- rownames(cigomp)
rownames(cigomp) <- NULL
cigomp$para <- gl(2, 11, labels=c("A","C"))
cigomp$trat <- rep(seq(1:11), times=2)
cigomp## est lwr upr par para trat
## 1 45.94900610 39.98610383 51.91190838 A.tr1 A 1
## 2 17.12143629 12.91464895 21.32822363 A.tr2 A 2
## 3 24.97611930 17.53961795 32.41262065 A.tr3 A 3
## 4 30.18312106 25.66333178 34.70291034 A.tr4 A 4
## 5 39.62887694 33.87486947 45.38288442 A.tr5 A 5
## 6 41.59437355 36.18767896 47.00106814 A.tr6 A 6
## 7 12.28029709 8.14377635 16.41681783 A.tr7 A 7
## 8 35.13487477 30.08302398 40.18672556 A.tr8 A 8
## 9 42.84773983 37.82101285 47.87446682 A.tr9 A 9
## 10 48.82822498 43.04841494 54.60803502 A.tr10 A 10
## 11 48.40773283 42.60311732 54.21234835 A.tr11 A 11
## 12 0.03316331 0.02529400 0.04103262 C.tr1 C 1
## 13 0.03017098 0.01847532 0.04186664 C.tr2 C 2
## 14 0.02122373 0.01288884 0.02955862 C.tr3 C 3
## 15 0.03920182 0.02757242 0.05083122 C.tr4 C 4
## 16 0.03183108 0.02369380 0.03996835 C.tr5 C 5
## 17 0.03604453 0.02702194 0.04506713 C.tr6 C 6
## 18 0.02611450 0.01375237 0.03847663 C.tr7 C 7
## 19 0.03578094 0.02605468 0.04550719 C.tr8 C 8
## 20 0.04246875 0.03171625 0.05322126 C.tr9 C 9
## 21 0.03624852 0.02785991 0.04463714 C.tr10 C 10
## 22 0.03576988 0.02746790 0.04407186 C.tr11 C 11segplot(trat~lwr+upr|para, data=cigomp,
scales=list(x="free"), #layout=c(NA,1),
draw.bands=FALSE, centers=est,
ylab="Treatment",
xlab=list(c(
"A",
"C"),
segments.fun=panel.arrows, ends="both",
angle=90, length=1, unit="mm"))
coefi = fixef(gompi)
xyplot(inc~tt|tr,data=datmf,layout=c(11,1),
panel=function(x,y,pn=panel.number(),...){
panel.xyplot(x,y,...)
panel.curve(coefi[pn]*(1-1.25*exp(-coefi[pn+11]*x)),...)
})
Efectos aleatorios (asintota en parcelas anidadas en los bloques)
ranef(gomp1)## Level: bq
## A.(Intercept)
## 1 1.2548300
## 2 0.4663248
## 3 -1.7211548
##
## Level: par %in% bq
## A.(Intercept)
## 1/101 -0.12382703
## 1/102 0.43306667
## 1/103 -0.97799278
## 1/104 -0.38146672
## 1/105 -1.39022977
## 1/106 1.66855108
## 1/107 -0.05005714
## 1/108 1.43676320
## 1/109 -3.67557021
## 1/110 3.33915238
## 1/111 3.38584563
## 2/201 -1.91954441
## 2/202 0.38069715
## 2/203 1.89488388
## 2/204 1.25427231
## 2/205 0.68585748
## 2/206 1.78862032
## 2/207 -0.15610814
## 2/208 1.44423383
## 2/209 2.84793560
## 2/210 -4.74814821
## 2/211 -2.11098242
## 3/301 2.04337144
## 3/302 -0.81376382
## 3/303 -0.91689110
## 3/304 -0.87280560
## 3/305 0.70437230
## 3/306 -3.45717140
## 3/307 0.20616528
## 3/308 -2.88099703
## 3/309 0.82763461
## 3/310 1.40899583
## 3/311 -1.27486321Contraste de coeficientes A (asintotas) entre tratamiento
summary(glht(gompi,linfct=matriz),test = adjusted("none"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: nlme.formula(model = inc ~ A * (exp(-exp(1.25 - C * tt))), data = datmf,
## fixed = A + C ~ 0 + tr, random = A ~ 1, start = sta, weights = varPower(form = ~fitted(.)))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## A1vs2 == 0 28.8276 3.3675 8.561 < 2e-16 ***
## A1vs3 == 0 20.9729 4.5731 4.586 4.52e-06 ***
## A1vs4 == 0 15.7659 3.4698 4.544 5.53e-06 ***
## A1vs5 == 0 6.3201 3.9076 1.617 0.105792
## A1vs6 == 0 4.3546 3.7792 1.152 0.249217
## A1vs7 == 0 33.6687 3.3444 10.067 < 2e-16 ***
## A1vs8 == 0 10.8141 3.6520 2.961 0.003065 **
## A1vs9 == 0 3.1013 3.6433 0.851 0.394640
## A1vs10 == 0 -2.8792 3.9174 -0.735 0.462349
## A1vs11 == 0 -2.4587 3.9267 -0.626 0.531214
## A2vs3 == 0 -7.8547 4.0426 -1.943 0.052019 .
## A2vs4 == 0 -13.0617 2.7359 -4.774 1.80e-06 ***
## A2vs5 == 0 -22.5074 3.2734 -6.876 6.16e-12 ***
## A2vs6 == 0 -24.4729 3.1192 -7.846 4.22e-15 ***
## A2vs7 == 0 4.8411 2.5706 1.883 0.059667 .
## A2vs8 == 0 -18.0134 2.9636 -6.078 1.21e-09 ***
## A2vs9 == 0 -25.7263 2.9531 -8.712 < 2e-16 ***
## A2vs10 == 0 -31.7068 3.2853 -9.651 < 2e-16 ***
## A2vs11 == 0 -31.2863 3.2964 -9.491 < 2e-16 ***
## A3vs4 == 0 -5.2070 4.1301 -1.261 0.207402
## A3vs5 == 0 -14.6528 4.5042 -3.253 0.001141 **
## A3vs6 == 0 -16.6183 4.3935 -3.782 0.000155 ***
## A3vs7 == 0 12.6958 4.0208 3.158 0.001591 **
## A3vs8 == 0 -10.1588 4.2844 -2.371 0.017734 *
## A3vs9 == 0 -17.8716 4.2773 -4.178 2.94e-05 ***
## A3vs10 == 0 -23.8521 4.5130 -5.285 1.26e-07 ***
## A3vs11 == 0 -23.4316 4.5211 -5.183 2.19e-07 ***
## A4vs5 == 0 -9.4458 3.3787 -2.796 0.005179 **
## A4vs6 == 0 -11.4113 3.2294 -3.534 0.000410 ***
## A4vs7 == 0 17.9028 2.7071 6.613 3.76e-11 ***
## A4vs8 == 0 -4.9518 3.0795 -1.608 0.107845
## A4vs9 == 0 -12.6646 3.0692 -4.126 3.69e-05 ***
## A4vs10 == 0 -18.6451 3.3901 -5.500 3.80e-08 ***
## A4vs11 == 0 -18.2246 3.4009 -5.359 8.38e-08 ***
## A5vs6 == 0 -1.9655 3.6958 -0.532 0.594848
## A5vs7 == 0 27.3486 3.2494 8.417 < 2e-16 ***
## A5vs8 == 0 4.4940 3.5655 1.260 0.207524
## A5vs9 == 0 -3.2189 3.5566 -0.905 0.365448
## A5vs10 == 0 -9.1993 3.8369 -2.398 0.016504 *
## A5vs11 == 0 -8.7789 3.8465 -2.282 0.022470 *
## A6vs7 == 0 29.3141 3.0942 9.474 < 2e-16 ***
## A6vs8 == 0 6.4595 3.4244 1.886 0.059252 .
## A6vs9 == 0 -1.2534 3.4151 -0.367 0.713613
## A6vs10 == 0 -7.2339 3.7061 -1.952 0.050955 .
## A6vs11 == 0 -6.8134 3.7160 -1.834 0.066724 .
## A7vs8 == 0 -22.8546 2.9371 -7.781 7.11e-15 ***
## A7vs9 == 0 -30.5674 2.9269 -10.444 < 2e-16 ***
## A7vs10 == 0 -36.5479 3.2618 -11.205 < 2e-16 ***
## A7vs11 == 0 -36.1274 3.2729 -11.038 < 2e-16 ***
## A8vs9 == 0 -7.7129 3.2737 -2.356 0.018474 *
## A8vs10 == 0 -13.6934 3.5763 -3.829 0.000129 ***
## A8vs11 == 0 -13.2729 3.5865 -3.701 0.000215 ***
## A9vs10 == 0 -5.9805 3.5674 -1.676 0.093652 .
## A9vs11 == 0 -5.5600 3.5776 -1.554 0.120158
## A10vs11 == 0 0.4205 3.8564 0.109 0.913173
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)Etapa final
Datos de cancro de la base del tallo:
dcbt = read.csv("cbt.csv", sep=",", h=TRUE ); head(dcbt)## ia sym tr bq par s.f i.f s.m i.m sp.f sp.m
## 1 test 0 1 1 101 21 70 63 96 38 82
## 2 az 1 2 1 102 3 15 10 36 25 36
## 3 az 1 3 1 103 5 25 15 52 25 37
## 4 az 1 4 1 104 9 35 26 76 32 43
## 5 az 1 5 1 105 11 40 34 80 34 54
## 6 az 1 6 1 106 9 40 32 80 28 50datmf<- datmf[order(datmf$bq, datmf$par),]# efectos fijos del tratamiento (11) Tr:1-11
af <- fixef(gompi)[1:11]
afn <- strsplit(names(af), "tr")
Af<- NULL # ef fijos (165)
for (i in 1:dim(datmf)[1]) {
for (j in 1:length(af)) if(datmf$tr[i]==afn[[j]][2]) Af[i]<- af[j] }
A_trat = Af[c(1,6, 11, 16, 21, 26, 31, 36, 41, 46, 51, 56, 61, 66, 71, 76, 81, 86, 91, 96, 101,106, 111, 116, 121, 126, 131, 136, 141, 146, 151, 156, 161)] # a_bq: Efecto aleatorio de bq (33) Tr: 11111111111,222,333...
a_bq <- c(rep(ranef(gompi)$"bq"[1,1], times=11),
rep(ranef(gompi)$"bq"[2,1], times=11),
rep(ranef(gompi)$"bq"[3,1], times=11)) # a_par: efecto aleatorio de parcela (33) [ordenado por bq]
a_par = ranef(gompi)$"par"[,1]Suma de los efectos para estimar la incidencia de cada parcela:
efectos = data.frame(dcbt$tr, dcbt$bq, dcbt$par, A_trat, a_bq, a_par)
efectos$Ass = A_trat + a_bq + a_par
names(efectos)[1:3] = c("tr", "bq", "par")
efectos## tr bq par A_trat a_bq a_par Ass
## 1 1 1 101 45.94901 1.336013 -1.309597e-08 47.285019
## 2 2 1 102 17.12144 1.336013 6.192002e-08 18.457449
## 3 3 1 103 24.97612 1.336013 -7.559157e-08 26.312132
## 4 4 1 104 30.18312 1.336013 -3.201459e-08 31.519134
## 5 5 1 105 39.62888 1.336013 -9.697363e-08 40.964889
## 6 6 1 106 41.59437 1.336013 1.084788e-07 42.930386
## 7 7 1 107 12.28030 1.336013 -5.779253e-08 13.616310
## 8 8 1 108 35.13487 1.336013 1.137608e-07 36.470887
## 9 9 1 109 42.84774 1.336013 -2.010800e-07 44.183752
## 10 10 1 110 48.82822 1.336013 1.775321e-07 50.164238
## 11 11 1 111 48.40773 1.336013 1.630028e-07 49.743746
## 12 1 2 201 45.94901 1.014381 -1.228781e-07 46.963387
## 13 2 2 202 17.12144 1.014381 4.103106e-08 18.135817
## 14 3 2 203 24.97612 1.014381 1.603773e-07 25.990500
## 15 4 2 204 30.18312 1.014381 9.007341e-08 31.197502
## 16 5 2 205 39.62888 1.014381 5.909760e-08 40.643258
## 17 6 2 206 41.59437 1.014381 1.167716e-07 42.608754
## 18 7 2 207 12.28030 1.014381 -4.772909e-08 13.294678
## 19 8 2 208 35.13487 1.014381 1.141068e-07 36.149256
## 20 9 2 209 42.84774 1.014381 1.234681e-07 43.862121
## 21 10 2 210 48.82822 1.014381 -2.895816e-07 49.842605
## 22 11 2 211 48.40773 1.014381 -1.322557e-07 49.422113
## 23 1 3 301 45.94901 -2.350393 1.359741e-07 43.598613
## 24 2 3 302 17.12144 -2.350393 -1.029511e-07 14.771043
## 25 3 3 303 24.97612 -2.350393 -8.478574e-08 22.625726
## 26 4 3 304 30.18312 -2.350393 -5.805882e-08 27.832728
## 27 5 3 305 39.62888 -2.350393 3.787602e-08 37.278484
## 28 6 3 306 41.59437 -2.350393 -2.252503e-07 39.243980
## 29 7 3 307 12.28030 -2.350393 1.055216e-07 9.929904
## 30 8 3 308 35.13487 -2.350393 -2.278676e-07 32.784481
## 31 9 3 309 42.84774 -2.350393 7.761182e-08 40.497347
## 32 10 3 310 48.82822 -2.350393 1.120494e-07 46.477832
## 33 11 3 311 48.40773 -2.350393 -3.074714e-08 46.057339Combinación de datos: efectos (asintotas » incidencias maximas por parcela) + evaluaciones de severidad de cancro
finaldat = merge(dcbt, efectos, by = "par") # NA's match
finaldat = select(finaldat, par, ia, sym, s.m, Ass)
head(finaldat)## par ia sym s.m Ass
## 1 101 test 0 63 47.28502
## 2 102 az 1 10 18.45745
## 3 103 az 1 15 26.31213
## 4 104 az 1 26 31.51913
## 5 105 az 1 34 40.96489
## 6 106 az 1 32 42.93039Ajuste exponencial para severidad de cancro en la base del tallo en función de la incidencia máxima de manchas de phoma (A_Gompertz)
mod_fin = nls(s.m ~ a*exp(Ass*b), data=finaldat, start=list(a= 5, b=0.04))with(finaldat,
plot(Ass,s.m,
xlab='A-Gompertz',
ylab='Severidad Cancro de la base del tallo (%)',
xlim=c(9,55), ylim=c(0,65),
pch=sym
))
curve(coef(mod_fin)[[1]]*exp(x*coef(mod_fin)[[2]]), xlim = c(10, 50), col = "red", add = TRUE)
co <- format(abs(coef(mod_fin)), digits=2)
text(15, 50, label=substitute(y==a^(b*x), list(a=co[1], b=co[2])),cex=0.8)
p9 = plot(mod_fin, cex=0.7)
p10 = qqnorm(mod_fin, ~ resid(., type = "p"), abline = c(0, 1), cex=0.7)
grid.arrange(p9, p10, ncol = 2)