[R] Constructing formula using nls

From: Bock, Michael <MBock_at_arcadis-us.com>
Date: Tue 21 Jun 2005 - 23:39:43 EST


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]]))

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Jun 21 23:50:38 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:55 EST