estimating change points

Gregoire, Timothy

unread,

Aug 7, 2012, 4:22:20 AM8/7/12

to for...@googlegroups.com

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

G&V sampling text: http://crcpress.com/product/isbn/9781584883708

From: for...@googlegroups.com [mailto:for...@googlegroups.com]
Sent: Tuesday, August 07, 2012 3:00 AM
To: Digest Recipients
Subject: Digest for for...@googlegroups.com - 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

Phil Radtke <pradt...@gmail.com> 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

"Saunders, Michael R" <msau...@purdue.edu> 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.
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.

Two-phase_Change-point_Models.pdf

Phil Radtke

unread,

Aug 7, 2012, 5:53:22 PM8/7/12

to for...@googlegroups.com

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

2012 Southern Mensurationists Meeting V3.pdf

Lauri Mehtätalo

unread,

Aug 9, 2012, 1:49:00 AM8/9/12

to for...@googlegroups.com

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 lauri.m...@uef.fi

http://joyx.joensuu.fi/~lamehtat/ http://joyx.joensuu.fi/~lamehtat/

050-442-2962 +358-50-442-2962

======================================================================

To view this discussion on the web visit https://groups.google.com/d/msg/forestr/-/0xsuwtrzCb0J.

Reply all

Reply to author

Forward