RE: Digest for forestr@googlegroups.com - 2 Messages in 1 Topic
Gregoire, Timothy
Hi Phil and Lauri,
I must be getting old (no, I am getting old!) , hence my acerbic response below.
Hell, if one cannot account for within subject correlations, then estimating a “best” fitting mean function is a ruse. What’s the point? What metric can one use that would be impervious to the within-subject correlation?
Tim
Timothy G. Gregoire
J. P. Weyerhaeuser Professor of Forest Management
School of Forestry & Environmental Studies, Yale University
360 Prospect Street, New Haven, CT 06511-2104 U.S.A.
office: 1.203.432.9398 mobile: 1.203.508.4014, fax: 1.203.432.3809
G&V sampling text: http://crcpress.com/product/isbn/9781584883708
From: for...@googlegroups.com [mailto:for...@googlegroups.com]
Sent: Thursday, August 09, 2012 3:45 AM
To: Digest Recipients
Subject: Digest for for...@googlegroups.com - 2 Messages in 1 Topic
Group: http://groups.google.com/group/forestr/topics
§ estimating change points [2 Updates]
Phil Radtke <pradt...@gmail.com> Aug 07 02:53PM -0700
Thanks Tim.
I'm not using mixed effects because it is forbidden in the challenge posted
by Lewis Jordan, part of his 2012 Mike Strub Challenge (attached). I
honestly don't think I have a chance at this sort of thing since A) I am
not an expert in taper modeling, B) my best answer comes from bagging,
which is also forbidden, and C) Quang Cao goes to Southern Mensurationist
meetings and he'll almost certainly win! Anyway, Sharma sent me his SAS
code and I was able to estimate join points with it. The *a priori* join
points suggested by Harold (0.1 and 0.75 values of h/HT) didn't work too
well with this species (Pinus wallichiana). The trees have a very steep
butt flare below 4.5 feet, with basal diameters more often exceeding 2*dbh
than not.
Not sure how far you'd have to go from Delhi to find one of these trees,
but at least you are on the right continent. Best.
Phil
On Tuesday, August 7, 2012 4:22:20 AM UTC-4, Tim wrote:
> forestr+u...@googlegroups.com <javascript:>.
> For more options, visit this group at
> http://groups.google.com/group/forestr?hl=en.
On Tuesday, August 7, 2012 4:22:20 AM UTC-4, Tim wrote:
"Lauri Mehtätalo" <lauri.m...@uef.fi> Aug 09 08:49AM +0300
Hi Phil and Tim,
I completerly agree on Tims point on using pre-defined change points.
However, the code below has an example how to technically fit the segmented Max & Burkharts function (using trunk to implement the segmentation). In addition, I included an example how to fit the switching model of Valentine and Gregoire by using log and logit-
transformations to parameters that have restrictions,
The codes are related to our recent paper de Miguel, S., Mehtätalo, L., Shater, Z., Kraid, B., and Pukkala, T. 2012. Evaluating marginal and conditional predictions of taper models in the absence of calibration data. Candadian Journal of Forest research 42(7):
1383-1394
maxburkhart<-function(D,H,h,a1,a2,b1,b2,b3,b4) {
D*sqrt(b1*(h/H-1)+b2*((h/H)^2-1)+b3*(a1-h/H)^2*trunc(a1-h/H+1)+
b4*(a2-h/H)^2*trunc(a2-h/H+1))
}
# Switching model of Valentine and greroire 2001
# Restrictions for the parameters
# theta4>= pi*D/4/(H*(H-BH))
Valentinegregoire<-function(D,H,h,theta1=0.326,
logittheta3=logit(0.0508*H/(H-(pi*(D/2)^2)/(theta4*(H-BH)))),
ltheta4=log(4.584-pi*D^2/4/(H*(H-BH))),
lambda1=1.881,lambda2=5.299,alpha1=1.203,alpha2=0.353,
BH=1.37) {
theta4<-exp(ltheta4)+pi*D^2/4/(H*(H-BH))
C<-H-(pi*(D/2)^2)/(theta4*(H-BH))
theta3<-inv.logit(logittheta3)*C/H
# cat(theta3," ",C/H,theta4,"\n")
S5<-theta1*H/(1+(h/(theta3*H))^lambda1)
S6<-(h/C)^lambda2/(1+(h/C)^lambda2)
a<-pi*(D/2)^2*((H-h)/(H-1.3))^(alpha1+S5)*((H-h)/(H-C))^(alpha2*S6)
sqrt(4*a/pi)
}
# MAx and Burkhart
# fit for diameter
fm1.mb<-nls(d~maxburkhart(Dbh,Hht,h,a1,a2,b1,b2,b3,b4),
data=dat,
start=list(b1=-4.07,b2=2.05,b3=20.15,b4=-1.90,a1=0.159,a2=0.780))
# fit for diameter squared
fm2.mb<-nls(I(d^2)~(maxburkhart(Dbh,Hht,h,a1,a2,b1,b2,b3,b4))^2,
data=dat,
start=list(b1=-4.07,b2=2.05,b3=20.15,b4=-1.90,a1=0.159,a2=0.780))
# fit for diameter
fm1.VG01<-nls(d~Valentinegregoire(Dbh,Hht,h,theta1,logittheta3,ltheta4,lambda1,lambda2,alpha1,alpha2,BH=1.3),
data=dat,
start=list(theta1=0.326,logittheta3=-2.713226,ltheta4=1.320155,
lambda1=1.881,lambda2=5.299,
alpha1=1.203,alpha2=0.353))
# fit for diameter squared
fm2.VG01<-nls(I(d^2)~(Valentinegregoire(Dbh,Hht,h,theta1,logittheta3,ltheta4,lambda1,lambda2,alpha1,alpha2,BH=1.3)^2),
data=dat,
start=list(theta1=0.326,logittheta3=-2.713226,ltheta4=1.320155,
lambda1=1.881,lambda2=5.299,
alpha1=1.203,alpha2=0.353))
Lauri
========================================================================
Lauri Mehtätalo Lauri Mehtätalo
Yliopistotutkija, Metsäsuunnittelu Senior researcher, Forest planning
Itä-Suomen yliopisto University of Eastern Finland
Joensuun kampus Joensuu campus
Metsätieteiden osasto School of Forest Sciences
PL 111 P.O.Box 111
80101 Joensuu 80101 Joensuu, Finland
lauri.m...@uef.fi<mailto:lauri.m...@joensuu.fi>
lauri.m...@uef.fi<mailto:lauri.m...@joensuu.fi>
http://joyx.joensuu.fi/~lamehtat/
http://joyx.joensuu.fi/~lamehtat/
050-442-2962 +358-50-442-2962
======================================================================
From: for...@googlegroups.com [mailto:for...@googlegroups.com] On Behalf Of Phil Radtke
Sent: Wednesday, August 08, 2012 12:53 AM
To: for...@googlegroups.com
Subject: Re: estimating change points
Thanks Tim.
I'm not using mixed effects because it is forbidden in the challenge posted by Lewis Jordan, part of his 2012 Mike Strub Challenge (attached). I honestly don't think I have a chance at this sort of thing since A) I am not an expert in taper modeling, B) my
best answer comes from bagging, which is also forbidden, and C) Quang Cao goes to Southern Mensurationist meetings and he'll almost certainly win! Anyway, Sharma sent me his SAS code and I was able to estimate join points with it. The a priori join points
suggested by Harold (0.1 and 0.75 values of h/HT) didn't work too well with this species (Pinus wallichiana). The trees have a very steep butt flare below 4.5 feet, with basal diameters more often exceeding 2*dbh than not.
Not sure how far you'd have to go from Delhi to find one of these trees, but at least you are on the right continent. Best.
Phil
On Tuesday, August 7, 2012 4:22:20 AM UTC-4, Tim wrote:
Phil,
There may be something in this biblio of help, esp. the Hinckley 1969 paper in Technometrics. If I recall correctly (and I may not be, as I have not looked at his paper in decades) he took a brute force ML approach. Calc the likelihood for each possible pair
of join points and chose the pair at which it is maximized.
Personally, I would prefer not to estimate them, but rather set them a priori, and condition other param estimates and subsequent inference on the selected join points. In trees, we know something happens at the base of the live crown, therefore I would set
your alpha2 there. Then if you really want to estimate alpha1 < alpha2, perhaps with one less param you will avoid singular grad problem.
Yet, how can you use nls anyway? Dbh-ht pairs within tree are not independent, so don't you need to use gnls (nlme pkg) with a corStruct... perhaps a CAR correlation function?
Let me know how you make out.
Tim (from Delhi)
Timothy G. Gregoire
J. P. Weyerhaeuser Professor of Forest Management
School of Forestry & Environmental Studies, Yale University
360 Prospect Street, New Haven, CT 06511-2104 U.S.A.
office: 1.203.432.9398 mobile: 1.203.508.4014, fax: 1.203.432.3809
G&V sampling text: http://crcpress.com/product/isbn/9781584883708
From: for...@googlegroups.com<javascript:> [mailto:for...@googlegroups.com<javascript:>]
Sent: Tuesday, August 07, 2012 3:00 AM
To: Digest Recipients
Subject: Digest for for...@googlegroups.com<javascript:> - 2 Messages in 2 Topics
Today's Topic Summary
Group: http://groups.google.com/group/forestr/topics
* estimating join points in segmented taper model [1 Update]
* Question on manuscript methods [1 Update]
estimating join points in segmented taper model<http://groups.google.com/group/forestr/t/15227f4fd2d1309a>
Phil Radtke <pradt...@gmail.com<javascript:>> Aug 06 10:41AM -0700
Hello Forest-Rs,
Can anyone tell me how to estimate break/join points of a segmented taper
function (maybe using nls?). I'm trying to estimate join points alpha1 and
alpha2 that connect model segments following Mahadev Sharma's model (8) in
Forest Science 49(2):324. <https://lh5.googleusercontent.com/-CMb71JOgInQ/UB__5_-K2iI/AAAAAAAAAAY/C-lNT8PnarI/s1600/Sharma8.gif>
With my code below, all I get is a "singular gradient" error message. :(
-Phil
(data set "strub.dat" attached)
# R-script ###################################################
# read data for the Mike Strub 2012 Challenge
pt <- read.table(file="strub.dat",header=TRUE)
pt$x <- pt$h/pt$HT
pt$y <- (pt$d/pt$DBH)^2
# Join point definitions coded as binary indicators
i1 <- function(x,alpha) { ifelse(x >= alpha,1.0,0) }
i2 <- function(x,alpha) { ifelse(x < alpha,1.0,0) }
i3 <- function(x,alpha) { ifelse(x < alpha,1.0,0) }
# model 8 from Sharma & Burkhart (For Sci 49(2): pg324)
fit <- nls(y ~ (a1*(x-1) + a2*(x^2-1))*i1(x,alpha1) +
a4*(alpha1-x)^2*i2(x,alpha1) + a6*(alpha2-x)^2*i3(x,alpha2),
data=pt,start=list(a1=-3.3, a2=1.6, a4=-1.7, a6=41,alpha1=0.75,alpha2=0.1))
# end R-script ###################################################
Thanks!
Question on manuscript methods<http://groups.google.com/group/forestr/t/4e3cb6c2ad800203>
"Saunders, Michael R" <msau...@purdue.edu<javascript:>> Aug 06 02:12PM -0400
All:
I have not had any help from other R-listservs on this question, so I was hopeful the Forestry group could provide me some guidance. I hate to cause undue problems for authors as a reviewer, but this analysis just does not pass the "straight-face" test for
me.
The experimental design that was presented included sampling sections of a log (3 sections per tree) from 2-3 trees per treatment (2 or 3 treatments per stand) in 3 stands. The independent variables were counts of internal structures in the log stems (e.g.,
structure A, structure B, etc).
The authors used SAS to "link" the structures to one another. So for example,
Structure A = a + b*Structure B
They used random effects for stand, treatment (nested in stand), and tree (nested in treatment and experiment). They present that only the intercept had random effects associated with it (although I question that approach from a biological standpoint).
So my question:
With only 18 trees total (54 log sections), is it even possible to fit a 3-level mixed-effects model? My experience with multi-level models is that you have to have a large enough sample size within the innermost experiment units (maybe not all of them, but
some of them) to estimate the variance of the parameter(s). Here, there is only 3 in that innermost unit and even with a linear model, I find it hard to believe that the variance would be estimated well.
Thanks in advance,
Mike R. Saunders
Assistant Professor of Hardwood Silviculture
Department of Forestry and Natural Resources
Purdue University
715 State Street
West Lafayette, IN 47907
765-430-1440
You received this message because you are subscribed to the Google Group forestr.
You can post via email<javascript:>.
To unsubscribe from this group, send<javascript:> an empty message.
For more options, visit<http://groups.google.com/group/forestr/topics> this group.
--
You received this message because you are subscribed to the Google Groups "Forest-R" group.
To post to this group, send email to
for...@googlegroups.com<javascript:>.
To unsubscribe from this group, send email to
forestr+u...@googlegroups.com<javascript:>.
For more options, visit this group at
http://groups.google.com/group/forestr?hl=en.
On Tuesday, August 7, 2012 4:22:20 AM UTC-4, Tim wrote:
Phil,
There may be something in this biblio of help, esp. the Hinckley 1969 paper in Technometrics. If I recall correctly (and I may not be, as I have not looked at his paper in decades) he took a brute force ML approach. Calc the likelihood for each possible pair
of join points and chose the pair at which it is maximized.
Personally, I would prefer not to estimate them, but rather set them a priori, and condition other param estimates and subsequent inference on the selected join points. In trees, we know something happens at the base of the live crown, therefore I would set
your alpha2 there. Then if you really want to estimate alpha1 < alpha2, perhaps with one less param you will avoid singular grad problem.
Yet, how can you use nls anyway? Dbh-ht pairs within tree are not independent, so don't you need to use gnls (nlme pkg) with a corStruct... perhaps a CAR correlation function?
Let me know how you make out.
Tim (from Delhi)
Timothy G. Gregoire
J. P. Weyerhaeuser Professor of Forest Management
School of Forestry & Environmental Studies, Yale University
360 Prospect Street, New Haven, CT 06511-2104 U.S.A.
office: 1.203.432.9398 mobile: 1.203.508.4014, fax: 1.203.432.3809
timothy....@yale.edu<javascript:>
G&V sampling text: http://crcpress.com/product/isbn/9781584883708
From: for...@googlegroups.com<javascript:> [mailto:for...@googlegroups.com<javascript:>]
Sent: Tuesday, August 07, 2012 3:00 AM
To: Digest Recipients
Subject: Digest for for...@googlegroups.com<javascript:> - 2 Messages in 2 Topics
Today's Topic Summary
Group: http://groups.google.com/group/forestr/topics
* estimating join points in segmented taper model [1 Update]
* Question on manuscript methods [1 Update]
estimating join points in segmented taper model<http://groups.google.com/group/forestr/t/15227f4fd2d1309a>
Phil Radtke <pradt...@gmail.com<javascript:>> Aug 06 10:41AM -0700
Hello Forest-Rs,
Can anyone tell me how to estimate break/join points of a segmented taper
function (maybe using nls?). I'm trying to estimate join points alpha1 and
alpha2 that connect model segments following Mahadev Sharma's model (8) in
Forest Science 49(2):324. <https://lh5.googleusercontent.com/-CMb71JOgInQ/UB__5_-K2iI/AAAAAAAAAAY/C-lNT8PnarI/s1600/Sharma8.gif>
With my code below, all I get is a "singular gradient" error message. :(
-Phil
(data set "strub.dat" attached)
# R-script ###################################################
# read data for the Mike Strub 2012 Challenge
pt <- read.table(file="strub.dat",header=TRUE)
pt$x <- pt$h/pt$HT
pt$y <- (pt$d/pt$DBH)^2
# Join point definitions coded as binary indicators
i1 <- function(x,alpha) { ifelse(x >= alpha,1.0,0) }
i2 <- function(x,alpha) { ifelse(x < alpha,1.0,0) }
i3 <- function(x,alpha) { ifelse(x < alpha,1.0,0) }
# model 8 from Sharma & Burkhart (For Sci 49(2): pg324)
fit <- nls(y ~ (a1*(x-1) + a2*(x^2-1))*i1(x,alpha1) +
a4*(alpha1-x)^2*i2(x,alpha1) + a6*(alpha2-x)^2*i3(x,alpha2),
data=pt,start=list(a1=-3.3, a2=1.6, a4=-1.7, a6=41,alpha1=0.75,alpha2=0.1))
# end R-script ###################################################
Thanks!
Question on manuscript methods<http://groups.google.com/group/forestr/t/4e3cb6c2ad800203>
"Saunders, Michael R" <msau...@purdue.edu<javascript:>> Aug 06 02:12PM -0400
All:
I have not had any help from other R-listservs on this question, so I was hopeful the Forestry group could provide me some guidance. I hate to cause undue problems for authors as a reviewer, but this analysis just does not pass the "straight-face" test for
me.
The experimental design that was presented included sampling sections of a log (3 sections per tree) from 2-3 trees per treatment (2 or 3 treatments per stand) in 3 stands. The independent variables were counts of internal structures in the log stems (e.g.,
structure A, structure B, etc).
The authors used SAS to "link" the structures to one another. So for example,
Structure A = a + b*Structure B
They used random effects for stand, treatment (nested in stand), and tree (nested in treatment and experiment). They present that only the intercept had random effects associated with it (although I question that approach from a biological standpoint).
So my question:
With only 18 trees total (54 log sections), is it even possible to fit a 3-level mixed-effects model? My experience with multi-level models is that you have to have a large enough sample size within the innermost experiment units (maybe not all of them, but
some of them) to estimate the variance of the parameter(s). Here, there is only 3 in that innermost unit and even with a linear model, I find it hard to believe that the variance would be estimated well.
Thanks in advance,
Mike R. Saunders
Assistant Professor of Hardwood Silviculture
Department of Forestry and Natural Resources
Purdue University
715 State Street
West Lafayette, IN 47907
765-430-1440
You received this message because you are subscribed to the Google Group forestr.
You can post via email<javascript:>.
To unsubscribe from this group, send<javascript:> an empty message.
For more options, visit<http://groups.google.com/group/forestr/topics> this group.
--
You received this message because you are subscribed to the Google Groups "Forest-R" group.
To post to this group, send email to
for...@googlegroups.com<javascript:>.
To unsubscribe from this group, send email to
forestr+u...@googlegroups.com<javascript:>.
You received this message because you are subscribed to the Google Group forestr.
You can post via email.
To unsubscribe from this group,
send an empty message.
For more options, visit this group.
--
You received this message because you are subscribed to the Google Groups "Forest-R" group.
To post to this group, send email to for...@googlegroups.com.
To unsubscribe from this group, send email to
forestr+u...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/forestr?hl=en.
