My data: I have raw data points that form a logit style curve as if they were a time series. Which is to say they form 3 distinct lines with 3 distinct slopes in backwards z pattern. A certain class of my data looks essentially flat to the eye with marginal oscillation. What is important to me is the x value at which the state change is occurring, in other words, the break point Use of segmented(): Segmented does a very good job of capturing the breakpoints and fitting three distinct slopes, i.e. linear models. It does not, however give me Pr(>|t|) for the break point coefficients. I need to answer the question H:0 Beta0=Beta with a certainty metric, i.e. probability statistic. This is especially important for my, flat looking data class. davies.test() question: davies.test() only excepts lm() or glm() objects as input. If I run segmented to find 1 breakpoint instead of 2, I get a totally bogus answer. Without knowing the breakpoints, how is this test able to assess the proper breakpiont? It appears to only give 1 best breakpoint, which is not consistent with the breakpoints found by segmented(). Also, is K the number of breakpoints or the number of iterations that it evaluates the breakpoint? Thanks in advance. lmfit<-glm(TotRad_KW~HRRPUA_kWm2,data=d1) davies.test(lmfit,seg.Z=~HRRPUA_kWm2,k=1000,alternative="less", beta0=0,dispersion=NULL) Davies' test for a change in the slope data: Model = gaussian , link = identity formula = TotRad_KW ~ HRRPUA_kWm2 segmented variable = HRRPUA_kWm2 `Best' at = 561.205, n.points = 1000, p-value < 2.2e-16 alternative hypothesis: less segments <- segmented(lmfit, seg.Z=~HRRPUA_kWm2,psi=c(475,550)) summary(segments) ***Regression Model with Segmented Relationship(s)*** Call: segmented.glm(obj = lmfit, seg.Z = ~HRRPUA_kWm2, psi = c(475, 550)) Estimated Break-Point(s): Est. St.Err psi1.HRRP 430.2 4.087 psi2.HRRP 484.6 3.077 t value for the gap-variable(s) V: 0 0 Meaningful coefficients of the linear terms: Estimate Std. Error t value Pr(>|t|) (Intercept) -38.6993 274.7666 -0.141 0.8891 HRRPUA_kWm2 1.4297 0.7472 1.914 0.0668 . U1.HRRP 42.2884 4.7696 8.866 NA U2.HRRP -40.5897 4.7123 -8.614 NA --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 6934.706) Null deviance: 70776718 on 31 degrees of freedom Residual deviance: 180302 on 26 degrees of freedom AIC: 377.19 Convergence attained in 2 iterations with relative change -1.662839e-14 -- Greg Cohn Forestry Technician USDA Forest Service Fire Lab [[alternative HTML version deleted]]
Vito Muggeo
2012-Nov-19 10:49 UTC
[R] Interpretation of davies.test() in segmented package
dear Greg, > It does not, however give me > Pr(>|t|) for the break point coefficients. Of course a test H_0: breakpoint=0 is meaningless.. > I need to answer the question > H:0 Beta0=Beta with a certainty metric, sorry, who is your "Beta0"? davies.test() tests the hypothesis H0: leftSlope=rightSlope which implies "diffSlope=0" and then "the breakpoint does not exist". K in davies.test() means the number of evaluation points used to compute the approximate p-value.. Please contact me off list if you need more details (given detailed questions) vito Il 16/11/2012 20.57, Greg Cohn ha scritto:> My data: > I have raw data points that form a logit style curve as if they were a time > series. Which is to say they form 3 distinct lines with 3 distinct slopes > in backwards z pattern. A certain class of my data looks essentially flat > to the eye with marginal oscillation. What is important to me is the x > value at which the state change is occurring, in other words, the break > point > > Use of segmented(): > Segmented does a very good job of capturing the breakpoints and fitting > three distinct slopes, i.e. linear models. It does not, however give me > Pr(>|t|) for the break point coefficients. I need to answer the question > H:0 Beta0=Beta with a certainty metric, i.e. probability statistic. This is > especially important for my, flat looking data class. > > davies.test() question: > davies.test() only excepts lm() or glm() objects as input. If I run > segmented to find 1 breakpoint instead of 2, I get a totally bogus answer. > Without knowing the breakpoints, how is this test able to assess the proper > breakpiont? It appears to only give 1 best breakpoint, which is not > consistent with the breakpoints found by segmented(). Also, is K the number > of breakpoints or the number of iterations that it evaluates the breakpoint? > > > Thanks in advance. > > > > lmfit<-glm(TotRad_KW~HRRPUA_kWm2,data=d1) > davies.test(lmfit,seg.Z=~HRRPUA_kWm2,k=1000,alternative="less", > beta0=0,dispersion=NULL) > > Davies' test for a change in the slope > > data: Model = gaussian , link = identity > formula = TotRad_KW ~ HRRPUA_kWm2 > segmented variable = HRRPUA_kWm2 > `Best' at = 561.205, n.points = 1000, p-value < 2.2e-16 > alternative hypothesis: less > > > > > segments <- segmented(lmfit, seg.Z=~HRRPUA_kWm2,psi=c(475,550)) > summary(segments) > > ***Regression Model with Segmented Relationship(s)*** > > Call: > segmented.glm(obj = lmfit, seg.Z = ~HRRPUA_kWm2, psi = c(475, > 550)) > > Estimated Break-Point(s): > Est. St.Err > psi1.HRRP 430.2 4.087 > psi2.HRRP 484.6 3.077 > > t value for the gap-variable(s) V: 0 0 > > Meaningful coefficients of the linear terms: > Estimate Std. Error t value Pr(>|t|) > (Intercept) -38.6993 274.7666 -0.141 0.8891 > HRRPUA_kWm2 1.4297 0.7472 1.914 0.0668 . > U1.HRRP 42.2884 4.7696 8.866 NA > U2.HRRP -40.5897 4.7123 -8.614 NA > --- > Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 > (Dispersion parameter for gaussian family taken to be 6934.706) > > Null deviance: 70776718 on 31 degrees of freedom > Residual deviance: 180302 on 26 degrees of freedom > AIC: 377.19 > > Convergence attained in 2 iterations with relative change -1.662839e-14 > > > > > ______________________________________________ > R-help at r-project.org mailing list > stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >-- ===================================Vito M.R. Muggeo Dip.to Sc Statist e Matem `Vianelli' Universit? di Palermo viale delle Scienze, edificio 13 90128 Palermo - ITALY tel: 091 23895240 fax: 091 485726 dssm.unipa.it/vmuggeo