I am relatively new to R and am looking for advice, ideas or both...
I have a data set that consists of pathogen population sizes on
individual plant units in an experimental field plot. However, in
order to estimate the pathogen population sizes I had to destroy the
plant unit and could not determine if that plant unit became diseased
or to what extent it would have become diseased. I collected disease
data from each plot and have estimates of the disease in an
experimental plot several days after the pathogen population sizes were
estimated. From this data I would like to determine a dose-response
relationship in the field between the pathogen and host.
i. e.
P(disease) = sum P(disease | pathogen population size) x P(pathogen
population size) ... I think
I determined that the pathogen population sizes among plant units are
well described by the lognormal distribution. I am also assuming that
P(disease | pathogen population size) is could be described by the
normal distribution described by two parameters lam and tau...my dose
response.
I would like to use nonlinear least squares to estimate lam and tau. I
included the R code and some data below. I'm not sure if this code is
correct, but it give reasonable estimates of lam and tau for the first
data set. However, I included another data set where the estimates of
lam and tau do not make sense. I have never done nonlinear regression.
Is the code correct? Is this a sound approach to estimating a dose
response in the field? Does anyone have other ideas about how to
approach this data? I am currently doing this analysis ad hoc, but if
the results are promise I might like to do some planned experiments.
R 2.0.1
mac os 10.2
Thanks in advance for any responses
Ken
Example data...
mn is the mean of the log transformed pathogen population size
ss is the standard deviation of the log transformed pathogen population
size
dap27 is the disease 7 days after the pathogen population sizes were
estimated
dap31 is the disease 11 days after the pathogen population sizes were
estimated
mn ss dap27 dap31
6.762 0.492 0.582 0.567
4.382 1.63 0.393 0.394
2.472 2.449 0.181 0.19
6.66 0.495 0.698 0.541
4.282 1.401 0.345 0.435
1.08 2.423 0.16 0.196
6.636 0.362 0.704 0.526
3.638 1.445 0.12 0.509
1.854 2.07 0.148 0.075
fm1 <- nls(dap27 ~ pnorm((mn - lam) / ((ss^2) + (tau^2))) , data = dg,
start=c(lam = 6, tau = 2), trace = TRUE)
summary(fm1)
plot(dap27 ~ mn, dg)
kfc<-as.numeric(kf<-lm(dg$ss ~ dg$mn)$coefficients)
mm<-mean(dg$mn, na.rm = T)
vv<-var(dg$mn, na.rm = T)
for(i in 1:5){
n<-100
mn<-rnorm(n, mm, sqrt(vv))
xx<- mn * kfc[2] + kfc[1]; lter<-rnorm(n, 0, 1)
ss<-xx + lter
newdat<-data.frame(mn,ss)
lines(lowess(newdat$mn , predict(fm1, newdat), f = 2/3), col = i)
}
The nonlinear regression did not fit for this data set (or at least
dap27).
mn ss dap27 dap31
0.31 0.94 0.01 0.31
0.09 3.33 0.01 .
3.3 2.43 0.07 0.22
5.9 1.78 0.46 0.91
6.03 1.26 0.33 0.96
4.7 1.96 0.09 0.67
2.67 2.04 0.06 0.51
4.7 1.78 0.31 0.93
3.07 2.74 0.02 0.39
5.13 2.6 0.27 0.68
3.38 2.24 0.05 0.44
5.9 1.7 0.49 .
5.06 2.18 0.08 0.54
2.99 2.64 0.09 0.38
5.41 2.73 0.08 0.46
1.65 2.96 0.04 0.39
