# [R] Help with a simple subroutine

Steven T. Yen @tyen @end|ng |rom ntu@edu@tw
Fri Sep 9 16:13:48 CEST 2022

```Thanks to all. It was just programming error. The following now works.
Essentially, to impose non-negative restrictions I estimated the natural
logs of two parameters and then do exponential transformation to uncover
the parameters, with mathematical transformation (delta method) for
their standard errors.

I believe deltaMethod {car} does that sort of things. I have yet to look
up if it does that for nonlinear regression object. But since the
exponential transformation is a simple transformation (with a derivative
equal to itself), I tried to program my own. Thanks to all.

> obj<-cbp11.pooled
> j<-grep("gamma",names(obj\$est),value=TRUE); j
[1] "log.gamma1" "log.gamma2"
> obj\$estimate[j]
log.gamma1 log.gamma2
-1.82378   -1.11313
> obj\$stat\$vb[j,j]
log.gamma1 log.gamma2
log.gamma1  0.0842252  0.0138778
log.gamma2  0.0138778  0.0793592
> mydelta <- function(obj,j){
+ # *******************************************
+ # Delta method for exponential transformation
+ # *******************************************
+   b<-obj\$estimate[j]
+   v<-obj\$stat\$vb[j,j]; v
+   gamma<-exp(b)
+   db<-gamma
+   vgamma<-db^2*v
+   sgamma<-sqrt(diag(vgamma))
+   t<-gamma/sgamma
+   df<-n<-obj\$stat\$n
+   p<-2*(1-pt(abs(t),df))
+   list(gamma=gamma,sgamma=sgamma,b=b,t=t,p=p)
+ }
> v<-mydelta(obj,j)
> v\$b
log.gamma1 log.gamma2
-1.82378   -1.11313
> v\$gamma
log.gamma1 log.gamma2
0.161414   0.328529
> v\$sgamma
log.gamma1 log.gamma2
0.0468449  0.0925490
> v\$t
log.gamma1 log.gamma2
3.44571    3.54978
> v\$p
log.gamma1  log.gamma2
0.000574108 0.000388996
>

On 9/9/2022 8:39 PM, Ebert,Timothy Aaron wrote:

> If t = 1/sqrt(v[2,2]) and there is no code to change the value of v[2,2] and no code to change to a different cell why would you get two different values?
>
> One approach to debugging is to make a small example and see if the code output matches (line-by-line) the output you get from doing a manual calculation. If they do not match you at least know where to start looking for a problem. Hand calculation can use pencil and paper or Excel or other tools. It is a tedious task but very effective.
>
> Tim
>
> -----Original Message-----
> From: R-help <r-help-bounces using r-project.org> On Behalf Of Ivan Krylov
> Sent: Friday, September 9, 2022 5:03 AM
> To: Steven T. Yen <styen using ntu.edu.tw>
> Cc: R-help Mailing List <r-help using r-project.org>
> Subject: Re: [R] Help with a simple subroutine
>
> [External Email]
>
> В Fri, 9 Sep 2022 16:46:00 +0800
> "Steven T. Yen" <styen using ntu.edu.tw> пишет:
>
>> I am expecting the line  t<-gamma/sgamma to produce two different
>> values. But I confirm that it is doing tt<-gamma[1]/sgamma[1]
> No, it just happens that gamma[1]/sgamma[1] is the same as gamma[2]/sgamma[2], subject to rounding errors:
>
>> +   gamma<-exp(b)
>> +   vgamma<-gamma^2*v[2,2]
>> +   sgamma<-sqrt(vgamma)
>> +   t<-gamma/sgamma
> t = gamma / sgamma = gamma / sqrt(gamma^2 * v[2,2]) =
>    = gamma / (abs(gamma) * sqrt(v[2,2])) = (given gamma = exp(b) > 0)
>    = 1 / sqrt(v[2,2]).
>
> --
> Best regards,
> Ivan
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see