Search code examples
rvectordecomposition

Determination of k value using for (variable in vector) { in R


I am trying to find the value of k which is a constant in Decompostion models

Wt = W0 e^-kt

I am using R script to calculate it. However, I am having trouble in selection of the data set.

Example data

Site   Treatment    Month    Repli    WD
W2         0          0      2        1
W3         0          0      1        .9
W7         0          0      3        .7
W12        0          0      4        .7
W2         0          1      2        .4
W3         0          1      1         .3
W7         0          1      3        .3 
W12        0          1      4        .2 
W2         0          3      2        .2
W3         0          3      1        .1
W7         0          3      3        .1  
W12        0          3      4        .1 
W2         0          4      2        .9
W3         0          4      1         .9
W7         0          4      3         .8
W12        0          4      4         .8


k=numeric(max(Cs$Repli))

pred=array(0,dim=c(4,max(Cs$Repli)))  
resid=array(0,dim=c(4,max(Cs$Repli)))
obs=array(0,dim=c(4,max(Cs$Repli)))

Now the part I am having a problem Nonlinear reggression

for(i in 1:max(Cs$Repli)){
  s1=subset(Cs,Cs$Repli==i)

In my original data it selects only four data set like this (which is fine) May be because of having max Repli value?

Site      Treatment   Month    Repli
W12           0         0       4
W12           0         1       4
W12           0         3       4
W12           0         5       4

Then I run the below command then I get the value of K for these 4 data. Now I want to select W2 or any other. How can I do that? I tried different things but I couldnt get the results. I tried using %in% but (i) from above script is also being used below so I cant use %in%. Any help would be really appreciated.

   Month=s1$Month               
   WD=s1$WD         
   nonlin = nls(WD ~ 1*exp(-k*Month), trace=TRUE, start = list(k = .01))
   summary(nonlin)          
   pr=predict(nonlin)       
   res=residuals(nonlin)        
   k[i]=coef(nonlin)[1]
   pred[,i]=pr
   resid[,i]=res
   obs[,i]=WD   
}   #end of nls for loop

I can provide original data.


Solution

  • If you want to check which subsets you have processed, you may use an additionally variable W_subsets to store processed subsets:

    # before the nls loop
    W_subsets <- rep(list(NA), max(Cs$Repli))
    
    # inside the nls loop
    W_subsets[[i]] <- s1
    

    As I may see, the output of W_subsets is all right:

    [[1]]
        Site    Treatment   Month   Repli   WD
    2   W3          0           0       1   0.9
    6   W3          0           1       1   0.3
    10  W3          0           3       1   0.1
    14  W3          0           4       1   0.9
    
    [[2]]
        Site    Treatment   Month   Repli   WD
    1   W2          0           0       2   1.0
    5   W2          0           1       2   0.4
    9   W2          0           3       2   0.2
    13  W2          0           4       2   0.9
    ... etc. to W12