I need to do some simulations based on a nls model. I have been able to
use nls to determine the model and I have almost gotten to the point at
which I can construct the function. When I get the function done I will
use uniroot to solve the function at a specific "y" value for x. I
found
an example in the archive using the search tools but I can seem to load
Args with the fitted parameters. The function fdose was constructed
manually based on the fitted parameters, ndose is the "dynamic"
version
of the function that will change when I fit the model using different
data sets. Thanks in advance for any help you can offer.
Here is the code to generate the model function:
library(boot)
library(lattice)
setwd("C:/Documents and Settings/mbock/My Documents/Mink")
#Mink <- read.csv("MinkStatsa.csv")
#Mink <- subset(Mink, QTEQW != "NR")
Mink <-
structure(list(Authors = structure(as.integer(c(1, 1, 1, 1, 1,
1, 1, 8, 3, 3, 3, 3, 6, 6, 6, 6, 6, 2, 12, 10, 10, 4, 4, 9, 7,
7, 7, 11, 11, 11, 11, 11, 11, 11, 11, 11, 4, 5, 5, 5, 5, 5)), .Label
c("Aulerich and Ringer 1977",
"Aulerich et al. 1985", "Bleavins et al. 1980",
"Brunstrom et al. 2001",
"Bursian et al. 2003", "Den Boer 1984", "Heaton
1992", "Jensen et al.
1977",
"Kakela et al. 2002", "Kihlstrom et al. 1992", "Restum
et al. 1999",
"Wren et al. 1987"), class = "factor"), Form =
structure(as.integer(c(3,
4, 4, 4, 1, 2, 4, 10, 2, 2, 2, 2, 11, 11, 11, 8, 8, 4, 4, 9,
4, 9, 9, 2, 6, 7, 5, 6, 7, 5, 6, 7, 5, 6, 7, 5, 17, 14, 13, 16,
15, 12)), .Label = c("Aroclor 1221", "Aroclor 1242",
"Aroclor
1242:1248:1254",
"Aroclor 1254", "Carp (high)", "Carp (low)",
"Carp (medium)",
"Clophen A30", "Clophen A50", "Clophen A50/A60",
"Clophen A60",
"Goldfish/carp (high)", "Goldfish/carp (low-med)",
"Goldfish/carp
(low)",
"Goldfish/carp (med-high)", "Goldfish/carp (medium)",
"Seal blubber
extract"
), class = "factor"), ADD = c(3913.0435, 723.9382, 1491.0537,
1956.5217, 260.8696, 260.8696, 260.8696, 3300, 937.5, 1875, 3750,
7500, 25.2, 6076, 2025, 6076, 2025, 307.4866, 180, 2094.2408,
1307.815, 81.2348, 267.0623, 826, 130, 260, 320, 58.8235, 117.6471,
235.2941, 58.8235, 117.6471, 235.2941, 58.7544, 120.048, 252.5253,
53.1856, 36, 63, 103, 169, 414), ADDTEQ = c(55.582, 18.2961,
37.6834, 49.4472, 0.0265, 1.1447, 6.593, 913.8127, 4.1139, 8.2277,
16.4555, 32.9109, 11.8593, 2859.4033, 952.9776, 19.8479, 6.6149,
7.7711, 4.5491, 62.0308, 33.0639, 2.3558, 7.4777, 45.4545, 4.6523,
7.8409, 12.118, 2.0393, 4.0785, 8.1571, 2.0393, 4.0785, 8.1571,
2.0369, 4.1618, 8.7544, 0.0914, 0.367, 0.588, 0.978, 1.686, 7.671
), WPCB = structure(as.integer(c(16, 22, 9, 10, 1, 2, 13, 20,
6, 11, 19, 23, 3, 24, 18, 21, 12, 17, 7, 15, 8, 4, 14, 5, 25,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
25)), .Label = c("0.2678", "0.2998", "0.6144",
"0.8402", "0.9493",
"1.0775", "1.854", "12.7837", "15.4173",
"19.8375", "2.1549",
"2.3273", "2.6976", "2.7621", "20.2492",
"25.8584", "3.0412",
"37.8885", "4.3098", "52.5134", "6.983",
"7.4854", "8.6196",
"93.9869", "NR"), class = "factor"), QTEQW =
c(0.5168, 0.1787,
0.3682, 0.479, 2e-04, 0.008, 0.0644, 6.8409, 0.0287, 0.0575,
0.1149, 0.2295, 0.0858, 19.8437, 6.76, 0.1473, 0.0491, 0.0744,
0.0444, 0.6962, 0.3149, 0.0251, 0.0805, 0.123, 0.053, 0.0806,
0.1311, 0.0183, 0.0366, 0.0732, 0.0183, 0.0367, 0.0734, 0.0183,
0.0374, 0.0787, 0.0016, 0.0047, 0.0072, 0.0119, 0.0205, 0.0907
), WTEQa = c(0.2251, 0.0873, 0.1797, 0.235, 1e-04, 0.0033, 0.0314,
0.9427, 0.0118, 0.0235, 0.0471, 0.0941, 0.0074, 1.4622, 0.524,
0.1515, 0.0505, 0.0368, 0.0217, 0.2226, 0.1561, 0.0109, 0.0331,
0.0334, 0.0262, 0.0352, 0.0552, 0.0079, 0.0158, 0.0316, 0.0079,
0.0159, 0.0317, 0.0079, 0.0162, 0.034, 4e-04, 0.0015, 0.0021,
0.003, 0.0048, 0.0204), Surv = c(0, 0, 0, 0, 235, 203, 0, 0,
0, 0, 0, 0, 99, 0, 0, 0, 0, 0, 111, 0, 0, 56.78, 0, 73, 25, 12,
0, 147, 99, 37, 92, 4, 7, 71, 4, 0, 104, 92, 70, 114, 147, 52
), Surva = c(0, 0, 0, 0, 100, 100, 0, 0, 0, 0, 0, 0, 99, 0, 0,
0, 0, 0, 100, 0, 0, 56.77647059, 0, 73, 25, 12, 0, 100, 99, 37,
92, 4, 7, 71, 4, 0, 100, 92, 70, 100, 100, 52), Effect = c(100,
100, 100, 100, 0, 0, 100, 100, 100, 100, 100, 100, 1, 100, 100,
100, 100, 100, 0, 100, 100, 43.22352941, 100, 27, 75, 88, 100,
0, 1, 63, 8, 96, 93, 29, 96, 100, 0, 8, 30, 0, 0, 48)), .Names
c("Authors",
"Form", "ADD", "ADDTEQ", "WPCB",
"QTEQW", "WTEQa", "Surv", "Surva",
"Effect"), row.names = c("1", "2", "3",
"4", "5", "6", "7", "8",
"9", "10", "11", "12", "13",
"14", "15", "16", "17", "18",
"19",
"20", "21", "22", "23", "24",
"25", "26", "27", "28", "29",
"30",
"31", "32", "33", "34", "35",
"36", "37", "38", "39", "40",
"41",
"42"), class = "data.frame")
TS <- (nls( Effect/100 ~ (Po +((1-Po)/(1+exp(-(a+b*log(QTEQW)))))),
data = Mink,,start = list(a = 5,b = 2,Po = 0),trace = FALSE ) )
fdose = function (x)
{(-0.01899614
+((1-(-0.01899614))/(1+exp(-(5.88112869+1.71111*log(x))))))}
plot(Effect/100 ~ log(QTEQW), Mink)
tt <- seq(-8,2,length = 100)
lines(tt,predict(TS,list (QTEQW = exp(tt))))
#code to creat dynamic function
alist(all.vars(summary(TS)$formula[[3]]))
xx <- all.vars(summary(TS)$formula[[3]])
xxx <- vector("list", length(xx))
names(xxx) <- xx
Args <- do.call("alist", xxx)
ndose <- as.function(c(Args, summary(TS)$formula[[3]]))