CODEWITHSUNDEEP
r
Mateusz1981
Mateusz1981

Reputation: 1867

PROC MODEL in SAS write in R

I have a SAS code as below:

proc model;
  parms 
     X1 X2 X3 a2 a3 a4
;
  S = prfl1*X1 + prfl2*X2 + prfl3*X3 
;
      Z0= S-a2;
      R = Z0 + (Z0**2 + (2*a4*S)/(100**a3))**0.5;
      estH100 = S * (T**a3*(100**a3*R + a4))/(100**a3*(T**a3*R + a4));

      FIT est100
    START = (
   X1=20 X2=20 X3=20  a2=34.700396595 a3=1.7086613571 a4=8901.6700482
    )
        / OLS DATA=a startiter maxiter=10000 converge = 0.000001
          out=pns outresid outpredict outactual outest=pnsest4;

RUN;

I programed the same in R

nls(estH100 = (prfl1*X1 + prfl2*X2 + prfl3*X3) * (T**a3*(100**a3*
(((prfl1*X1 + prfl2*X2 + prfl3*X3)-a2) + (((prfl1*X1 + prfl2*X2 + prfl3*X3)-
a2)**2 + (2*a4*(prfl1*X1 + prfl2*X2 + prfl3*X3))/(100**a3))**0.5) + 
a4))/(100**a3*(T**a3*(((prfl1*X1 + prfl2*X2 + prfl3*X3)-a2) + (((prfl1*X1 + 
prfl2*X2 + prfl3*X3)-a2)**2 + (2*a4*(prfl1*X1 + prfl2*X2 + prfl3*X3))/(100**a3))**0.5) + a4))
, data = a, start = (a2=0.1, a3=2, a4=4, X1=20, X2 = 20, X3 = 20 ))

and it works. However, I do not like the structure as I have in R code. I have to replace S, Z0 and R whenever they appear in the function. Is it possible to write these in a similar way as it is in SAS????

Upvotes: 1

Views: 503

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 270195

1) function Define a function and put the calculations in it:

f <- function(a2, a3, a4, X1, X2, X3, prf11, prf12, prf13, T) {
    S <- prfl1*X1 + prfl2*X2 + prfl3*X3
    Z0 <- S-a2
    R <- Z0 + (Z0^2 + (2*a4*S)/(T0^a3))^0.5
    S * (T^a3*(100^a3*R + a4))/(100^a3*(T^a3*R + a4))
}
nls(estH100 ~ f(a2, a3, a4, X1, X2, X3, prf11, prf12, prf13, T),
  data = a, start = list(a2 = 0.1, a3 = 2, a4 = 4, X1 = 20, X2 = 20, X3 = 20 ))

2) substitute An alternative is to use substitute like this:

S <- quote(prfl1*X1 + prfl2*X2 + prfl3*X3)
Z0 <- substitute(S-a2, list(S = S))
R <- substitute(Z0 + (Z0^2 + (2*a4*S)/(T0^a3))^0.5, list(Z0 = Z0))
fo <- substitute(estH100 ~ S * (T^a3*(100^a3*R + a4))/(100^a3*(T^a3*R + a4)), 
 list(S = S, R = R))

At this point fo is the following and can be used in nls:

> fo
estH100 ~ (prfl1 * X1 + prfl2 * X2 + prfl3 * X3) * (T^a3 * (100^a3 * 
    (prfl1 * X1 + prfl2 * X2 + prfl3 * X3 - a2 + ((prfl1 * X1 + 
        prfl2 * X2 + prfl3 * X3 - a2)^2 + (2 * a4 * S)/(T0^a3))^0.5) + 
    a4))/(100^a3 * (T^a3 * (prfl1 * X1 + prfl2 * X2 + prfl3 * 
    X3 - a2 + ((prfl1 * X1 + prfl2 * X2 + prfl3 * X3 - a2)^2 + 
    (2 * a4 * S)/(T0^a3))^0.5) + a4))


> nls(fo, data = a, start = list(a2 = 0.1, a3 = 2, a4 = 4, X1 = 20, X2 = 20, X3 = 20 ))

Upvotes: 1

Related Questions