Can someone help me with this simple replicable subroutine?
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]
Thanks.
> b<-v$est[j]; b
log.gamma1 log.gamma2
? -1.82378?? -1.11313
> v<-v$stat$vb[j,j]; v
?????????? log.gamma1 log.gamma2
log.gamma1? 0.0842252? 0.0138778
log.gamma2? 0.0138778? 0.0793592
> delta <- function(b,v){
+ # ***********************************************
+ # Delta method for exponential transformation
+ # ***********************************************
+?? df<-5140; #df<-nrow(mydata)
+?? gamma<-exp(b)
+?? vgamma<-gamma^2*v[2,2]
+?? sgamma<-sqrt(vgamma)
+?? t<-gamma/sgamma
+?? tt<-gamma[1]/sgamma[1]
+?? p<-2*(1-pt(abs(t),df))
+ list(gamma=gamma,sgamma=sgamma,b=b,t=t,p=p,tt=tt)
+ }
> options(digits=6)
> delta(b,v)$b
log.gamma1 log.gamma2
? -1.82378?? -1.11313
> delta(b,v)$gamma
log.gamma1 log.gamma2
? 0.161414?? 0.328529
> delta(b,v)$sgamma
log.gamma1 log.gamma2
?0.0454716? 0.0925490
> delta(b,v)$t
log.gamma1 log.gamma2
?? 3.54978??? 3.54978
> delta(b,v)$tt
log.gamma1
?? 3.54978
> delta(b,v)$p
?log.gamma1? log.gamma2
0.000389002 0.000389002
Help with a simple subroutine
5 messages · Steven Yen, Andrew Simmons, Ivan Krylov +1 more
I think you must've written your formula incorrectly. t can be simplified like this: t gamma/sgamma gamma/sqrt(vgamma) gamma/sqrt(gamma^2 * v[2, 2]) gamma/sqrt(gamma^2)/sqrt(v[2, 2]) gamma/gamma /sqrt(v[2, 2]) rep(1, length(gamma)) /sqrt(v[2, 2]) so I'm not surprised that both values of t are equal.
On Fri., Sep. 9, 2022, 04:46 Steven T. Yen, <styen at ntu.edu.tw> wrote:
Can someone help me with this simple replicable subroutine? 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] Thanks.
> b<-v$est[j]; b
log.gamma1 log.gamma2 -1.82378 -1.11313
> v<-v$stat$vb[j,j]; v
log.gamma1 log.gamma2 log.gamma1 0.0842252 0.0138778 log.gamma2 0.0138778 0.0793592
> delta <- function(b,v){
+ # *********************************************** + # Delta method for exponential transformation + # *********************************************** + df<-5140; #df<-nrow(mydata) + gamma<-exp(b) + vgamma<-gamma^2*v[2,2] + sgamma<-sqrt(vgamma) + t<-gamma/sgamma + tt<-gamma[1]/sgamma[1] + p<-2*(1-pt(abs(t),df)) + list(gamma=gamma,sgamma=sgamma,b=b,t=t,p=p,tt=tt) + }
> options(digits=6) > delta(b,v)$b
log.gamma1 log.gamma2 -1.82378 -1.11313
> delta(b,v)$gamma
log.gamma1 log.gamma2 0.161414 0.328529
> delta(b,v)$sgamma
log.gamma1 log.gamma2 0.0454716 0.0925490
> delta(b,v)$t
log.gamma1 log.gamma2
3.54978 3.54978
> delta(b,v)$tt
log.gamma1
3.54978
> delta(b,v)$p
log.gamma1 log.gamma2 0.000389002 0.000389002
______________________________________________ R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
? Fri, 9 Sep 2022 16:46:00 +0800 "Steven T. Yen" <styen at 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
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 at r-project.org> On Behalf Of Ivan Krylov Sent: Friday, September 9, 2022 5:03 AM To: Steven T. Yen <styen at ntu.edu.tw> Cc: R-help Mailing List <r-help at 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 at 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 at r-project.org mailing list -- To UNSUBSCRIBE and more, see https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-help&data=05%7C01%7Ctebert%40ufl.edu%7C96597c0724d947d3c1ea08da9242e073%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637983113179778561%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=EnaoCEmoUJNXsukvk8jqZVZ1tIveOeUCIX%2Bic5tMRLM%3D&reserved=0 PLEASE do read the posting guide https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r-project.org%2Fposting-guide.html&data=05%7C01%7Ctebert%40ufl.edu%7C96597c0724d947d3c1ea08da9242e073%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637983113179778561%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=18MG1xlJnQE%2Bqo54jNxYAIAGqSQC%2FFQOgTOkl7Ysvc8%3D&reserved=0 and provide commented, minimal, self-contained, reproducible code.
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 at r-project.org> On Behalf Of Ivan Krylov Sent: Friday, September 9, 2022 5:03 AM To: Steven T. Yen <styen at ntu.edu.tw> Cc: R-help Mailing List <r-help at 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 at 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 at r-project.org mailing list -- To UNSUBSCRIBE and more, see https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-help&data=05%7C01%7Ctebert%40ufl.edu%7C96597c0724d947d3c1ea08da9242e073%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637983113179778561%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=EnaoCEmoUJNXsukvk8jqZVZ1tIveOeUCIX%2Bic5tMRLM%3D&reserved=0 PLEASE do read the posting guide https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r-project.org%2Fposting-guide.html&data=05%7C01%7Ctebert%40ufl.edu%7C96597c0724d947d3c1ea08da9242e073%7C0d4da0f84a314d76ace60a62331e1b84%7C0%7C0%7C637983113179778561%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=18MG1xlJnQE%2Bqo54jNxYAIAGqSQC%2FFQOgTOkl7Ysvc8%3D&reserved=0 and provide commented, minimal, self-contained, reproducible code.