I'm having a bit of trouble understanding how to write a function and extract parameters from a matrix in R. This is similar to Non linear fit in r but I haven't quite figured it out.
As an example, I've already used a linear fit to extract the y intercept and slope from the same data and apply those as a new columns in the matrix:
library(tidyverse)
library(lazyeval)
design.mat <- cbind(1,c(340,340.33,340.66,341,342,343,346,349,352,355,358,368,378,388,398)) #design matrix with x values of time points
response.mat <- t(ts_rw[,36:50]) # transform rows to columns and select relevant points
reg <- lm.fit(design.mat, response.mat)$coefficients #linear fit to time points, extract intercept and slope (coefficients)
ts_rw <- cbind(ts_rw, t(reg)) #re-transform (intercept and slope for each row are added in new columns)
where ts_rw is a large rowwise matrix with each column a time point and each row an individual
head(ts_rw)
individual X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18
1 <NA> 0 0.77150 1.17260 1.62550 2.6667 3.0663 3.2323 3.3331 3.4293 3.5010 3.6000 3.6099 3.6941 4.1558 4.2500 4.3215 4.3215
2 PI665092 0 0.71872 1.00040 1.25040 1.9385 2.3136 2.5397 2.7147 2.8328 2.9164 3.0014 3.0828 3.1824 3.6377 3.7900 3.9024 4.0000
3 PI665092 0 0.58076 0.84209 1.03010 2.0000 2.6109 2.9331 3.1266 3.2663 3.3763 3.4662 3.5404 3.6437 4.0007 4.0794 4.1032 4.0828
with many more rows and columns.
As a simplified bit of data, from the columns of interest:
timepoints <- c(340,340.33,340.66,341,342,343,346,349,352,355,358,368,378,388,398)
row2 <- as.numeric(as.vector(ts_rw[2,36:50]))
row3 <- as.numeric(as.vector(ts_rw[3,36:50]))
print(timepoints)
[1] 340.00 340.33 340.66 341.00 342.00 343.00 346.00 349.00 352.00 355.00 358.00 368.00 378.00 388.00 398.00
print(row2)
[1] 3.56900 1.69690 1.35020 1.19020 0.96450 0.86805 0.77572 0.75023 0.82184 0.85174 0.84638 0.84638 0.78600 0.74063 0.71467
print(row3)
[1] 5.5473 3.6941 3.3142 3.1771 2.7650 2.6260 2.4220 2.4710 2.6328 2.6328 2.5301 2.5301 2.4710 2.4748 2.4443
The function for the non-linear fit is:
Value = a*(exp(-b*time))+c
I would like to apply this in place of the lm.fit, in order to extract a, b, and c for each individual. Any ideas on a good way to do this? Thank you.
---update
I think the trouble lies in implementing correct start values; if I try something like this:
x <- c(340,340.33,340.66,341,342,343,346,349,352,355,358,368,378,388,398)
y <- as.numeric(as.vector(ts_rw[2,36:50]))
df <- data.frame(x,y)
print(df)
x y
1 340.00 3.56900
2 340.33 1.69690
3 340.66 1.35020
4 341.00 1.19020
5 342.00 0.96450
6 343.00 0.86805
7 346.00 0.77572
8 349.00 0.75023
9 352.00 0.82184
10 355.00 0.85174
11 358.00 0.84638
12 368.00 0.84638
13 378.00 0.78600
14 388.00 0.74063
15 398.00 0.71467
m <- nls(y ~ I(a*(exp(-b*x))+c), data=df, start=list(a=max(y), b=6.4, c=0.8), alg="plinear",trace=T)
y_est<-predict(m,df$x)
plot(x,y)
lines(x,y_est)
summary(m)
attempting to calculate estimates for coefficients a,b, and c as in this question
summary(m)
Error in chol2inv(object$m$Rmat()) :
element (2, 2) is zero, so the inverse cannot be computed
The fit doesn't work out
17/12/20 addition:
Is there good way to implement nlsList with the design and response matrix approach? Something like
design.mat <- cbind(1,c(340,340.33,340.66,341,342,343,346,349,352,355,358,368,378,388,398))
response.mat <- t(ts_rw[,36:50])
reg <- nlsList(y ~ exp(-b *(x + d)),data=c(design.mat, response.mat),start=list(b=0.3, c = 0.5, d = -339))$coefficients
ts_rw <- cbind(ts_rw, t(reg))
only in this the data argument is invalid as it's looking for a standard xy dataframe.
Update 17/12/20
Reshaped dataframe to use in nlslist model format. In case it may help someone, here's a summary of where things are at, starting with the initial few individuals from the dataframe (which had to be renamed individually to avoid the model grouping them as replicates):
ts_rw %>% head(5)
accession X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16
1 id_disc 0 0.77150 1.17260 1.62550 2.6667 3.0663 3.2323 3.3331 3.4293 3.5010 3.6000 3.6099 3.6941 4.1558 4.2500
2 PI665092_43 0 0.71872 1.00040 1.25040 1.9385 2.3136 2.5397 2.7147 2.8328 2.9164 3.0014 3.0828 3.1824 3.6377 3.7900
3 PI665092_44 0 0.58076 0.84209 1.03010 2.0000 2.6109 2.9331 3.1266 3.2663 3.3763 3.4662 3.5404 3.6437 4.0007 4.0794
4 NSL54162_43 0 0.41834 0.57771 0.73795 1.2189 1.6075 1.9582 2.3481 2.6006 2.7776 2.9443 3.0010 3.1266 3.6893 4.0000
5 NSL54162_44 0 0.23645 0.42549 0.63412 1.3080 1.6542 1.9132 2.1887 2.4286 2.6260 2.7776 2.9250 3.0012 3.6656 3.9297
X17 X18 X19 X20 X21 X22 X23 X24 X25 X26 X27 X28 X29 X30 X31 X32
1 4.3215 4.3215 4.3102 4.3102 4.2883 4.2737 4.3030 4.2275 4.1664 4.1664 4.0921 4.0794 4.0794 4.0794 4.0794 4.0794
2 3.9024 4.0000 4.1032 4.1068 4.1104 4.1104 4.1248 4.1104 4.1032 4.0976 4.0009 4.0000 3.9065 3.9024 3.8193 3.7280
3 4.1032 4.0828 4.0828 4.1558 4.1833 4.2330 4.3649 4.4561 4.5465 4.5828 4.6664 4.7289 4.8193 4.8328 4.9034 5.0017
4 4.5411 5.0029 5.5473 5.9043 6.2219 6.5003 6.9043 7.1249 7.4437 7.5710 7.7771 7.8885 8.2500 8.4985 8.7504 9.2484
5 4.2737 4.6664 5.0038 5.3039 5.5038 5.5552 5.6249 5.6249 5.6041 5.6667 5.7036 5.8039 5.9979 6.1041 6.1104 6.2219
X33 X34 X35 X36 X37 X38 X39 X40 X41 X42 X43 X44 X45 X46
1 4.0022 4.0022 4.0000 4.0022 1.4130 1.2761 1.2147 1.0331 0.93759 0.88252 0.86885 0.84830 0.80672 0.72985
2 3.7026 3.7026 3.6023 3.5690 1.6969 1.3502 1.1902 0.9645 0.86805 0.77572 0.75023 0.82184 0.85174 0.84638
3 5.1032 5.2036 5.3649 5.5473 3.6941 3.3142 3.1771 2.7650 2.62600 2.42200 2.47100 2.63280 2.63280 2.53010
4 9.8558 10.5000 11.3310 11.8000 11.0000 10.8330 10.8000 10.3330 10.28400 9.71380 9.42720 9.33280 9.14200 8.83280
5 6.2495 6.3333 6.5003 6.6667 6.2854 6.0000 6.0000 5.7503 5.62490 5.50380 5.30390 5.30390 5.30390 5.10320
X47 X48 X49 npq_end npq_max npq_up_yint npq_slope_up npq_slope_down
1 0.62518 0.5408 0.48702 0.45223 4.3215 0.10449764 1.5837233 -2.5429803
2 0.84638 0.7860 0.74063 0.71467 4.1248 0.14077229 1.2092617 -2.2405433
3 2.53010 2.4710 2.47480 2.44430 5.5473 0.11330366 1.0048922 -2.2426498
4 8.75040 8.5701 8.33130 8.25000 11.8000 0.07942392 0.7117107 -0.9479913
5 5.18330 5.0921 5.00000 5.00170 6.6667 0.01161035 0.6279491 -0.6844091
###
ts_dark <- ts_rw[-c(2:35,51:54)] #subset out columns of interest for this model
#rename columns by time point (minutes)
colnames(ts_dark)[2:16] <- c("340","340.33","340.66","341","342","343","346","349","352","355","358","368","378","388","398")
#reshape to long format
ts_long <- ts_dark %>%
melt(id.var="accession") %>%
arrange(accession, variable)
ts_long$variable <- as.numeric(as.character(ts_long$variable)) #convert variable column from factor to numeric for model
#run model on all individuals
reg <- nlsList(value ~ exp(-b *(variable + d)) + c |accession, data=ts_long, start=list(b=0.3, c = 0.5, d = -339)) #fit model to dark time points, extract intercept and slope (coefficients)
The model only throws a singular gradient error for one individual and three "missing value or infinity" errors, out out 81.
