This is an automated email from the git hooks/post-receive script. tille pushed a commit to branch master in repository r-cran-segmented.
commit 9b5e86ef167ff65d291642b51ee83719c6047c62 Author: Andreas Tille <[email protected]> Date: Fri Sep 29 10:24:23 2017 +0200 New upstream version 0.5-2.2 --- DESCRIPTION | 12 +- MD5 | 62 ++++++----- NAMESPACE | 4 +- NEWS | 37 ++++++- R/davies.test.r | 3 +- R/intercept.r | 8 +- R/plot.segmented.R | 36 +++--- R/pscore.test.r | 277 +++++++++++++++++++++++++++++++++++++++++++++++ R/seg.Ar.fit.r | 3 + R/seg.control.R | 4 +- R/seg.def.fit.r | 62 +++++++++-- R/seg.glm.fit.r | 2 + R/seg.lm.fit.r | 6 +- R/segmented.Arima.r | 3 +- R/segmented.default.r | 43 ++++++-- R/segmented.glm.R | 4 +- R/segmented.lm.R | 44 +++++--- R/slope.R | 55 +++++++--- R/vcov.segmented.R | 2 +- data/down.rda | Bin 401 -> 401 bytes data/plant.rda | Bin 664 -> 664 bytes data/stagnant.rda | Bin 352 -> 352 bytes inst/CITATION | 4 +- man/broken.line.Rd | 2 +- man/davies.test.Rd | 5 +- man/intercept.Rd | 8 +- man/plot.segmented.Rd | 27 +++-- man/points.segmented.Rd | 3 +- man/pscore.test.rd | 82 ++++++++++++++ man/seg.control.Rd | 3 +- man/segmented-package.Rd | 14 ++- man/segmented.Rd | 5 +- man/slope.Rd | 11 +- 33 files changed, 686 insertions(+), 145 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3c76936..b01395a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,15 +1,15 @@ Package: segmented Type: Package -Title: Regression Models with Breakpoints/Changepoints Estimation -Version: 0.5-1.4 -Date: 2015-11-04 +Title: Regression Models with Break-Points / Change-Points Estimation +Version: 0.5-2.2 +Date: 2017-09-18 Authors@R: c(person(given = c("Vito","M.","R."), family = "Muggeo", role = c("aut", "cre"), email = "[email protected]")) Author: Vito M. R. Muggeo [aut, cre] Maintainer: Vito M. R. Muggeo <[email protected]> -Description: Given a regression model, segmented `updates' the model by adding one or more segmented (i.e., piecewise-linear) relationships. Several variables with multiple breakpoints are allowed. +Description: Given a regression model, segmented `updates' the model by adding one or more segmented (i.e., piece-wise linear) relationships. Several variables with multiple breakpoints are allowed. License: GPL NeedsCompilation: no -Packaged: 2015-11-04 10:27:33 UTC; user +Packaged: 2017-09-18 11:10:10 UTC; enea Repository: CRAN -Date/Publication: 2015-11-04 17:33:57 +Date/Publication: 2017-09-19 00:28:24 UTC diff --git a/MD5 b/MD5 index a1cd6fa..7812a1a 100644 --- a/MD5 +++ b/MD5 @@ -1,55 +1,57 @@ -76de6f0eaa3fa57f8cd0a5baf29d3a6b *DESCRIPTION -bcd7c2c3f0c346d4b5acb8627e138d57 *NAMESPACE -fd592e80edbb5f22268e54aeaf1069ec *NEWS +69f0f65bfbcc946d9586bac2c49d2543 *DESCRIPTION +0ffe8f412060912587d634afef274654 *NAMESPACE +a3d677cb26ef4ce49683d183048d9ac6 *NEWS cf61d7b238288be0a75ea4c8eb35f7fb *R/broken.line.r 249a33f124608ea43e541b016a3d8110 *R/confint.segmented.R -ea2089afb8de71b205b6775af7e26347 *R/davies.test.r +82cb289649e5db67189c39cb67cccc1f *R/davies.test.r a42345d81421b74cdb8c8159fb2db83d *R/draw.history.R -70957e886612b4982d6aeab29457411c *R/intercept.r +811b0c28b99a4fe9b03a3aa9edf69ea7 *R/intercept.r 5e225c969dc3f328224c702ef14971e1 *R/lines.segmented.R -b2fd57a9cc6ec4f42b59b4714aa5eb08 *R/plot.segmented.R +3ddeff1d297cd6bf639b93a7d411ccf4 *R/plot.segmented.R 5f850c48c5c18919b2524b68be9c93c1 *R/points.segmented.r 764fb065475ee2f14297cd2a92c687ad *R/predict.segmented.r 65e331633c5afec47f0c98ccc871f058 *R/print.segmented.R 50650e747dedabc28e07a8e1ab09fe16 *R/print.summary.segmented.R +c34eb7ace2718dded6772bb386b37fa3 *R/pscore.test.r 3133f13e7a7ff6dc277ff3bbf8321c79 *R/seg.Ar.fit.boot.r -46a33c2c9cd789be7d796b44bd46a0fd *R/seg.Ar.fit.r -1897a3de17c45e04d74c97ba862551e8 *R/seg.control.R +c66b36ae348b6fac6c3117ca41318f95 *R/seg.Ar.fit.r +f127c92f7aaeaf603a70bccec59e5d02 *R/seg.control.R 388cbb0ef6faf923beb3b13df9c4098e *R/seg.def.fit.boot.r -5ca50277ebd588c37d55da46aa075e71 *R/seg.def.fit.r +33419bf83d6f7cf8540b7f93a7b59eb4 *R/seg.def.fit.r 92ed650a2dba26a38b3c51997f26a22c *R/seg.glm.fit.boot.r -10fabbc8bdf949e91f18afd80266b90e *R/seg.glm.fit.r +4e13424552caf91acad5028bc8968365 *R/seg.glm.fit.r 195d6dd356b9bc87584210e07d877ce4 *R/seg.lm.fit.boot.r -5bb2070b1c3c8368a17c2a19e6b9da4b *R/seg.lm.fit.r -e0d61b7c8115d3344d67a5c0eaf82a49 *R/segmented.Arima.r +93c700743d40f3ddbd0ceb10f83c4b2e *R/seg.lm.fit.r +dc05717bcef9538f666068c5af1b64e5 *R/segmented.Arima.r 6784ac55ef0763e1f1923bab5a59835d *R/segmented.R -f7dc1b7f68fa31ed619cd6b06f355f2a *R/segmented.default.r -d17b1174c4a7b00a3248353d7be520e2 *R/segmented.glm.R -c6206e2eed21e6e727515233471bb8cd *R/segmented.lm.R -0b5ecf873f23c77c3c8d07bc29b1ef47 *R/slope.R +223899d10d8f492a9f2e54f90a2a83c8 *R/segmented.default.r +1a0fbd6d9fd490fb790da3806c68ba67 *R/segmented.glm.R +13b7af2baabab3d658e8bb924f0af600 *R/segmented.lm.R +e4f667cabee25328e891045bb8786c82 *R/slope.R ce83cd92e7ae3b7783257bea4a6fa3cb *R/summary.segmented.R -42fcc7e2cabc1a73538a0a9f819cbb85 *R/vcov.segmented.R -e7b6f31aec1edfcc1bbf74672f23682c *data/down.rda -810a8a84a7aa8a3314ae5e2caaa06bc2 *data/plant.rda -618d00e8ec040138ec9cc793fa2cae5a *data/stagnant.rda -6227e4f7236f49f3debc821d4c98f7d1 *inst/CITATION -9d64af4e32959d9c47aadb932600d285 *man/broken.line.Rd +fa2b6ae145168ca63999617dd351bde1 *R/vcov.segmented.R +c4a5d60070a20fc1472a4e0efc2e8f55 *data/down.rda +2eedec2d09709515e3e5891351381172 *data/plant.rda +f7f474367bcb4ef9a00c80aba9f2596b *data/stagnant.rda +3158df5caf091ede33d58eed15bdb230 *inst/CITATION +a890933c26d412e555795aff354c326c *man/broken.line.Rd e38b11eb6601e9544619dfac019d96a6 *man/confint.segmented.Rd -2b3c2cde4bb5c1aa330c0582f5af4899 *man/davies.test.Rd +16464e04983759a33507e1bca4567b63 *man/davies.test.Rd a1fd6bbde564db5be2a9dde1ecedfb2d *man/down.Rd bb4e4f96d8eee8acb68b5773ab2a5c5a *man/draw.history.Rd -ba04128336ebed0bba0c5d4d3b2007eb *man/intercept.Rd +5b8fa6f1223e85dfe081169cf96b19a2 *man/intercept.Rd 7a543b5123d9c64b5a4da6aba4e1d3ea *man/lines.segmented.Rd 925cd1c6b4a05d3fca632a7f93999d00 *man/plant.Rd -9c9f96b038c8ae276bd4a7faa4ef0f2d *man/plot.segmented.Rd -607f76fbc9ece591960b9bb7086d79ea *man/points.segmented.Rd +138e297b73d53673fabd9025142aa524 *man/plot.segmented.Rd +b4da8d785ea0e901bf670dada5dd9a56 *man/points.segmented.Rd 1ba7c089bae8411532d0f40ff4091c7b *man/predict.segmented.Rd c5ff81f292c40cdc317fe2ddfc99c4cd *man/print.segmented.Rd -43f6ec26b9c92365a4da177eb2a9ae1c *man/seg.control.Rd +ea04af66300b8bf636d6a08f2ef573d9 *man/pscore.test.rd +1af889f836c98e6f631ecd058f675df6 *man/seg.control.Rd 55f09cc33b69788ffaaad8aad83faa3a *man/seg.lm.fit.Rd -33f289104c87bb682f8fe8365da73b8e *man/segmented-package.Rd -039ac0d1beb750d74cd3f1b8d90b94cf *man/segmented.Rd -866f3920741a92dd16afa46d31a38ca8 *man/slope.Rd +fb7088872dfa30aeff1a55f778e10cc4 *man/segmented-package.Rd +563ffd883ecca7937eec4f2c2277f0ab *man/segmented.Rd +33cff252508df1d97eec9ecf9ec1f12f *man/slope.Rd 1c8095f4472b4cabf3dc3cae5132260f *man/stagnant.Rd 706dec93d7feda2b0e6ef5f878529db2 *man/summary.segmented.Rd 99a8c632f6534f724d479cd4bca71824 *man/vcov.segmented.Rd diff --git a/NAMESPACE b/NAMESPACE index 240507f..a4064df 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,9 +10,11 @@ importFrom("stats", "approx", "as.formula", "coef", "contrasts", "family", "residuals", "runif", "summary.glm", "summary.lm", "update", "update.formula", "vcov", "weights") importFrom("utils", "flush.console") +importFrom("stats", "complete.cases") +importFrom("stats", "add1") export(segmented, segmented.default, segmented.lm, segmented.glm, segmented.Arima, - broken.line ,confint.segmented,davies.test,draw.history, + broken.line ,confint.segmented,davies.test,pscore.test,draw.history, intercept,lines.segmented,plot.segmented,print.segmented, seg.control,seg.lm.fit,seg.glm.fit,seg.lm.fit.boot,seg.glm.fit.boot, seg.def.fit,seg.def.fit.boot, seg.Ar.fit,seg.Ar.fit.boot, diff --git a/NEWS b/NEWS index a14d813..155413e 100644 --- a/NEWS +++ b/NEWS @@ -5,13 +5,44 @@ ************************************* +=============== +version 0.5-2.2 +=============== + +* When there is a single covariate in the starting (g)lm, seg.Z can be missing when calling the segmented methods. +* bug fixed: plot.segmented(.., link=FALSE) did not work correctly (sometimes it returned an error) for glm fits with multiple breakpoints. +Weights were not handled appropriately by segmented.lm. + + +=============== +version 0.5-2.1 +=============== +* pscore.test() now works also for "glm" fits +* plot.segmented() now plots the partial residuals as "component + working residuals" (rather than Pearson residuals, relevant only for glm fits). +* segmented.default() now works for fits obtained by MASS::rlm(). + + +=============== +version 0.5-2.0 +=============== +* pscore.test() introduced. The function tests for a breakpoint using a (pseudo) score statistic which is more powerful than davies.test(), especially when the breakpoint is expected to be in the middle of the covariate range and the signal-to-noise ratio is high. +* argument 'digits' added in seg.control() to fix the number of digits of the breakpoint estimate during the iterative estimation algorithm. +* bug fixed: conf.level>0 in plot.segmented() did not work for objects returned by segmented.default(). + + +=============== +version 0.5-1.5 (not on CRAN) +=============== +* arguments 'gap' and 'show.gap' removed in intercept() and in plot.segmented(). (they are meaningless, as segmented() always returns joined piecewise lines, i.e. with no gaps). +* slope() and broken.line() (and then plot.segmented() which uses them) did not work for objects returned by segmented.default() (Thanks to Marcos Krull for reporting). + =============== version 0.5-1.4 =============== -* segmented.Arima() should be slightly faster as starting values are passed in arima() (via 'init') throughout the iterative process. +* segmented.Arima() should be slightly faster, as starting values are passed in arima() (via 'init') throughout the iterative process. * plot.segmented() is expected to work for objects returned by segmented.Arima. -* print.summary.segmented() does not print anymore the t-values for the gap coefficients (this information is meaningless as the gap coeffs are always set to zero in the returned model) +* print.summary.segmented() does not print anymore the t-values for the gap coefficients (this information is meaningless as the gap coeffs are always set to zero in the returned model). * Bug fixed: intercept() ignored argument 'rev.sgn'; points.segmented() missed argument 'transf'. @@ -20,7 +51,7 @@ version 0.5-1.3 (not on CRAN) =============== * plot.segmented() gains argument 'transf' to plot 'transf(values)' rather 'values' on the current plot. * print.summary.segmented() now uses round() rather than signif() when displaying the breakpoint estimate. -* Bug fixed: psi=NA was not working in the segmented.* methods; this bug was incidentally introduced in the last version (thanks to Bertrand Sudre for reporting). +* Bug fixed: psi=NA was not working in the segmented.* methods; this bug was incidentally introduced in the last version (thanks to Bertrand Sudre for first reporting that). =============== diff --git a/R/davies.test.r b/R/davies.test.r index bf14e03..38590c6 100644 --- a/R/davies.test.r +++ b/R/davies.test.r @@ -183,7 +183,8 @@ daviesGLM<-function(y, z, xreg, weights, offs, values=NULL, k, list.glm, alterna if(k<10) warnings("k>=10 is recommended") alternative <- match.arg(alternative) type <- match.arg(type) - if(length(all.vars(seg.Z))>1) warning("multiple segmented variables ignored in 'seg.Z'",call.=FALSE) + #if(length(all.vars(seg.Z))>1) warning("multiple segmented variables ignored in 'seg.Z'",call.=FALSE) + if(length(all.vars(seg.Z))>1) stop("Only a single segmented variable can be specified in 'seg.Z' ") isGLM<-"glm"%in%class(obj) Call<-mf<-obj$call mf$formula<-formula(obj) diff --git a/R/intercept.r b/R/intercept.r index e3b331a..5218f4b 100644 --- a/R/intercept.r +++ b/R/intercept.r @@ -1,4 +1,4 @@ -intercept<-function (ogg, parm, gap=TRUE, rev.sgn = FALSE, var.diff = FALSE, +intercept<-function (ogg, parm, rev.sgn = FALSE, var.diff = FALSE, digits = max(3, getOption("digits") - 3)){ #corregge in caso di no model intercept -- CHE VOLEVO DIRE?? #forse che adesso funziona se nel modello non c'e' l'interc. #-- @@ -44,7 +44,7 @@ intercept<-function (ogg, parm, gap=TRUE, rev.sgn = FALSE, var.diff = FALSE, nomi <- names(coef(ogg)) nomi <- nomi[-match(nomepsi, nomi)] Allpsi <- index <- vector(mode = "list", length = length(nomeZ)) - gapCoef<-summary.segmented(ogg)$gap +# gapCoef<-summary.segmented(ogg)$gap ##eliminato 10/11/15 Ris <- list() rev.sgn <- rep(rev.sgn, length.out = length(nomeZ)) if("(Intercept)"%in%names(coef(ogg))){ @@ -63,10 +63,10 @@ intercept<-function (ogg, parm, gap=TRUE, rev.sgn = FALSE, var.diff = FALSE, cof <- coef(ogg)[ind] alpha <- vector(length = length(ind)) #gapCoef.i<-gapCoef[grep(paste("\\.",nomeZ[i],"$",sep=""), rownames(gapCoef), value = FALSE),"Est."] - gapCoef.i<-gapCoef[f.U(rownames(gapCoef), nomeZ[i]) ,"Est."] +# gapCoef.i<-gapCoef[f.U(rownames(gapCoef), nomeZ[i]) ,"Est."] ###eliminato 10/11/15 for (j in 1:length(cof)) { alpha[j] <- alpha0 - Allpsi[[i]][j] * cof[j] - if(gap) alpha[j] <- alpha[j] - gapCoef.i[j] +# if(gap) alpha[j] <- alpha[j] - gapCoef.i[j] ###eliminato 10/11/15 alpha0 <- alpha[j] } diff --git a/R/plot.segmented.R b/R/plot.segmented.R index ba36cab..e903d34 100644 --- a/R/plot.segmented.R +++ b/R/plot.segmented.R @@ -1,7 +1,7 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, interc=TRUE, link = TRUE, res.col = 1, rev.sgn = FALSE, const = 0, shade=FALSE, rug=TRUE, dens.rug=FALSE, dens.col = grey(0.8), - show.gap=FALSE, transf=I, ...){ + transf=I, ...){ #funzione plot.segmented che consente di disegnare anche i pointwise CI f.U<-function(nomiU, term=NULL){ #trasforma i nomi dei coeff U (o V) nei nomi delle variabili corrispondenti @@ -26,7 +26,9 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, } #-------------- #se l'oggetto e' segmented.Arima il nome dell'eventuale interc va sostituito.. - if((all(class(x)==c("segmented", "Arima")))) names(x$coef)<-gsub("intercept", "(Intercept)", names(coef(x))) +# if((all(class(x)==c("segmented", "Arima")))) names(x$coef)<-gsub("intercept", "(Intercept)", names(coef(x))) + if(all(c("segmented", "Arima") %in% class(x))) names(x$coef)<-gsub("intercept", "(Intercept)", names(coef(x))) + #-------------- linkinv <- !link if (inherits(x, what = "glm", which = FALSE) && linkinv && !is.null(x$offset) && res) stop("residuals with offset on the response scale?") @@ -42,10 +44,10 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, else { term <- x$nameUV$Z } - } - else { - if (!term %in% x$nameUV$Z) - stop("invalid `term'") + } else { + dterm<- deparse(substitute(term)) + if(dterm %in% x$nameUV$Z) term<-dterm + if (! isTRUE(term %in% x$nameUV$Z)) stop("invalid `term'") } opz <- list(...) cols <- opz$col @@ -69,7 +71,8 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, ylabs <- opz$ylab if (length(ylabs) <= 0) ylabs <- paste("Effect of ", term, sep = " ") - a <- intercept(x, term, gap = show.gap)[[1]][, "Est."] + #a <- intercept(x, term, gap = show.gap)[[1]][, "Est."] + a <- intercept(x, term)[[1]][, "Est."] #Poiche' intercept() restituisce quantita' che includono sempre l'intercetta del modello, questa va eliminata se interc=FALSE if(!interc && ("(Intercept)" %in% names(coef(x)))) a<- a-coef(x)["(Intercept)"] b <- slope(x, term)[[1]][, "Est."] @@ -85,14 +88,17 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, n<-length(x$residuals) #fitted.values - Arima non ha "fitted.values", ma ha "residuals".. tipo<- if(inherits(x, what = "glm", which = FALSE) && link) "link" else "response" - vall<-sort(c(seq(min(val), max(val), l=120), est.psi)) + vall<-sort(c(seq(min(val), max(val), l=150), est.psi)) #ciValues<-predict.segmented(x, newdata=vall, se.fit=TRUE, type=tipo, level=conf.level) vall.list<-list(vall) names(vall.list)<-term - ciValues<-broken.line(x, vall.list, link=link, interc=interc, se.fit=TRUE) + if(conf.level>0) { - k.alpha<-if(inherits(x, what = c("glm","Arima"), which = FALSE)) abs(qnorm((1-conf.level)/2)) else abs(qt((1-conf.level)/2, x$df.residual)) + #k.alpha<-if(inherits(x, what = c("glm","Arima"), which = FALSE)) abs(qnorm((1-conf.level)/2)) else abs(qt((1-conf.level)/2, x$df.residual)) + #cambiato nella 0.5-2.0: + k.alpha<-if(inherits(x, what = "lm", which = FALSE)) abs(qt((1-conf.level)/2, x$df.residual)) else abs(qnorm((1-conf.level)/2)) + ciValues<-broken.line(x, vall.list, link=link, interc=interc, se.fit=TRUE) ciValues<-cbind(ciValues$fit, ciValues$fit- k.alpha*ciValues$se.fit, ciValues$fit + k.alpha*ciValues$se.fit) #---> transf... ciValues<-apply(ciValues, 2, transf) @@ -109,7 +115,9 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, b.ok1 <- c(b, b[length(b)]) y.val <- y.val1 <- a.ok1 + b.ok1 * val + const s <- 1:(length(val) - 1) - xvalues <- if(all(class(x)==c("segmented", "Arima"))) x$Z[,1] else x$model[, term] +# xvalues <- if(all(class(x)==c("segmented", "Arima"))) x$Z[,1] else x$model[, term] + xvalues <- if(all(c("segmented", "Arima") %in% class(x))) x$Z[,1] else x$model[, term] + if (rev.sgn) { val <- -val xvalues <- -xvalues @@ -128,8 +136,10 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, #broken.line(x, term, gap = show.gap, link = link) + resid(x, "response") + const fit0 + resid(x, "response") + const else x$family$linkinv(c(y.val, y.val1)) - xout <- sort(c(seq(val[1], val[length(val)], l = 120), val[-c(1, length(val))])) +# xout <- sort(c(seq(val[1], val[length(val)], l = 150), val[-c(1, length(val))])) + xout <- sort(c(seq(val[1], val[length(val)], l = 150), val[-c(1, length(val))],val[-c(1, length(val))]*1.005)) l <- approx(as.vector(m[, c(1, 3)]), as.vector(m[, c(2, 4)]), xout = xout) + val[length(val)]<-max(l$x) #aggiunto 11/09/17 id.group <- cut(l$x, val, FALSE, TRUE) yhat <- l$y xhat <- l$x @@ -190,7 +200,7 @@ plot.segmented<-function (x, term, add = FALSE, res = FALSE, conf.level = 0, fit <- c(y.val, y.val1) if (res) { ress <- if (inherits(x, what = "glm", which = FALSE)) - residuals(x, "working") * sqrt(x$weights) + residuals(x, "working") #* sqrt(x$weights) mgcv::gam() usa " ..*sqrt(x$weights)/mean(sqrt(x$weights))" else resid(x) #if(!is.null(x$offset)) ress<- ress - x$offset #fit <- broken.line(x, term, gap = show.gap, link = link, interc = TRUE) + ress + const diff --git a/R/pscore.test.r b/R/pscore.test.r new file mode 100644 index 0000000..bbcc912 --- /dev/null +++ b/R/pscore.test.r @@ -0,0 +1,277 @@ +`pscore.test` <- function(obj, seg.Z, k = 10, alternative = c("two.sided", "less", "greater"), + values=NULL, dispersion=NULL, df.t=NULL, more.break=FALSE) { +#------------------------------------------------------------------------------- +test.Sc2<-function(y, z, xreg, sigma=NULL, values=NULL, fn="pmax(x-p,0)", df.t="Inf", alternative, w=NULL){ +#xreg: la matrice del disegno del modello nullo. Se mancante viene assunta solo l'intercetta. +#Attenzione che se invXtX e xx vengono entrambe fornite, non viene fatto alcun controllo +#invXtX: {X'X}^{-1}. if missing it is computed from xreg +#sigma: the sd. If missing it is computed from data (under the *null* model) +#values: the values with respect to ones to compute the average term. If NULL 10 values from min(z) to max(z) are taken. + n<-length(y) + if(missing(xreg)) xreg<-cbind(rep(1,n)) + id.ok<-complete.cases(cbind(y,z,xreg)) + y<-y[id.ok] + z<-z[id.ok] + xreg<-xreg[id.ok,,drop=FALSE] + n<-length(y) + k=ncol(xreg) #per un modello ~1+x + if(is.null(values)) values<-seq(min(z), max(z), length=10) + n1<-length(values) + PSI<-matrix(values, nrow=n, ncol=n1, byrow=TRUE) #(era X2) matrice di valori di psi + if(is.matrix(z)) { + X1<-matrix(z[,1], nrow=n, ncol=n1, byrow=FALSE) + X2<-matrix(z[,2], nrow=n, ncol=n1, byrow=FALSE) + X<-eval(parse(text=fn), list(x=X1, y=X2, p=PSI)) #X<-pmax(X1-X2,0) + pmaxMedio<-rowMeans(X) + } else { + X1<-matrix(z, nrow=n, ncol=n1, byrow=FALSE) #matrice della variabile Z + if(length(fn)<=1){ + X<-eval(parse(text=fn), list(x=X1, p=PSI)) #X<-pmax(X1-X2,0) + pmaxMedio<-rowMeans(X) + } else { + pmaxMedio<-matrix(,n,length(fn)) + #list.X<-vector("list", length=length(fn)) + for(j in 1:length(fn)){ + #list.X[[j]]<-eval(parse(text=fn[j]), list(x=X1, p=PSI)) + X<-eval(parse(text=fn[[j]]), list(x=X1, p=PSI)) + pmaxMedio[,j]<-rowMeans(X) + } + } + } + if(is.null(w)){ + invXtX<-solve(crossprod(xreg)) + IA<-diag(n) - xreg%*%tcrossprod(invXtX, xreg) #I- hat matrix + r<-IA%*%y + sc<- sum(pmaxMedio*r) + v.s<-crossprod(pmaxMedio,IA)%*% pmaxMedio + } else { + invXtX<-solve(crossprod(sqrt(w)*xreg)) + IA<-diag(n) - xreg%*%tcrossprod(invXtX, xreg*w) #I-hat matrix + sc<-t(pmaxMedio*w) %*% IA %*% y + v.s<- t(pmaxMedio*w) %*% crossprod(t(IA)/sqrt(w))%*%(w*pmaxMedio) + } + +# ris<-if(length(fn)<=1) sc/(sigma*sqrt(v.s)) else drop(crossprod(sc,solve(v.s,sc)))/(sigma^2) +# if(length(fn)<=1 && cadj) ris<- sign(ris)*sqrt((ris^2)*(1-(3-(ris^2))/(2*n))) + ris<- drop(sc/(sigma*sqrt(v.s))) + df.t<-eval(parse(text=df.t)) + + pvalue<- switch(alternative, + less = pt(ris, df=df.t, lower.tail =TRUE) , + greater = pt(abs(ris), df=df.t, lower.tail =FALSE) , + two.sided = 2*pt(abs(ris), df=df.t, lower.tail =FALSE) + ) + #pvalue<- 2*pt(abs(ris), df=df.t, lower.tail =FALSE) + r<-c(ris, pvalue, pmaxMedio) + r + } +#------------------------------------------------------------------------------- +scGLM<-function(y, z, xreg, family, values = NULL, size=1, weights.var, + fn="pmax(x-p,0)", alternative=alternative){ +#score test for GLM +#size (only if family=binomial()) +#weights.var: weights to be used for variance computations. If missing the weights come from the null fit + output<-match.arg(output) + n<-length(y) + if(missing(xreg)) xreg<-cbind(rep(1,n)) + id.ok<-complete.cases(cbind(y,z,xreg)) + y<-y[id.ok] + z<-z[id.ok] + xreg<-xreg[id.ok,,drop=FALSE] + n<-length(y) + if(family$family=="poisson") size=1 + if(length(size)==1) size<-rep(size,n) + yN<-y/size + k=ncol(xreg) #per un modello ~1+x + if(is.null(values)) values<-seq(min(z), max(z), length=10) + n1<-length(values) + PSI<-matrix(values, nrow=n, ncol=n1, byrow=TRUE) #(era X2) matrice di valori di psi + X1<-matrix(z, nrow=n, ncol=n1, byrow=FALSE) #matrice della variabile Z + X<-eval(parse(text=fn), list(x=X1, p=PSI)) #X<-pmax(X1-X2,0) + pmaxMedio<-rowMeans(X) + + o<-glm.fit(yN, x=xreg, weights=size, family=family) + r<-y-(o$fitted*size) + sc<-drop(crossprod(r, pmaxMedio)) +# if(output=="unst.score") return(drop(sc)) + + p <- o$rank + Qr <- o$qr + COV <- chol2inv(Qr$qr[1:p, 1:p, drop = FALSE]) #vcov(glm(y~x, family=poisson)) + A<-xreg%*%COV%*%crossprod(xreg, diag(o$weights)) + h<- drop(tcrossprod(pmaxMedio, diag(n)- A)) + if(missing(weights.var)) weights.var<-o$weights + v.s<- drop(crossprod(h*sqrt(weights.var))) #t(h)%*%diag(exp(lp))%*%h + + ris<-if(length(fn)<=1) sc/sqrt(v.s) else drop(crossprod(sc,solve(v.s,sc))) +# if(output=="score") return(drop(ris)) + + pvalue<- switch(alternative, + less = pnorm(ris, lower.tail =TRUE) , + greater = pnorm(abs(ris), lower.tail =FALSE) , + two.sided = 2*pnorm(abs(ris), lower.tail =FALSE) + ) + +# pvalue<- if(length(fn)<=1) 2*pnorm(abs(ris), lower.tail =FALSE) else pchisq(ris,df=length(fn), lower.tail =FALSE) + # NB: se calcoli ris<-drop(t(sc)%*%solve(v.s,sc))/(length(fn)*sigma^2) devi usare pf(ris,df1=length(fn),df2=df.t, lower.tail =FALSE) + return(c(ris, pvalue)) + } +#---------------------------------------------------- + + fn="pmax(x-p,0)" +# if(inherits(obj, "glm")) stop("Currently only 'lm', or 'segmented-lm' models are allowed") + if(!inherits(obj, "lm")) stop("A 'lm', or 'segmented-lm' model is requested") + if(inherits(obj, "segmented") && length(obj$nameUV$Z)==1) seg.Z<- as.formula(paste("~", obj$nameUV$Z )) + if(!inherits(obj, "segmented") && length(all.vars(formula(obj)))==2) seg.Z<- as.formula(paste("~", all.vars(formula(obj))[2])) + if(class(seg.Z)!="formula") stop("'seg.Z' should be an one-sided formula") + name.Z <- all.vars(seg.Z) + if(length(name.Z)>1) stop("Only a single segmented variable can be specified in 'seg.Z' ") + + if(k<=1) stop("k>1 requested! k>=10 is recommended") + if(k<10) warnings("k>=10 is recommended") + alternative <- match.arg(alternative) + #if(length(all.vars(seg.Z))>1) warning("multiple segmented variables ignored in 'seg.Z'",call.=FALSE) + isGLM<-"glm"%in%class(obj) + + if(isGLM){ + if(is.null(dispersion)) dispersion<- summary.glm(obj)$dispersion + if(inherits(obj, "segmented")){ + mf<-model.frame(obj) + X0<-model.matrix(obj) + Z<-X0[,name.Z] + n<-length(Z) + if(is.null(values)) values<-seq(min(Z), max(Z), length=k) #values<-seq(sort(Z)[2], sort(Z)[(n - 1)], length = k) + n1<-length(values) + #escludi *tutte* le variabili psi (sia da X0 che dalla formula) + #X0<-X0[, -match(obj$nameUV$V, colnames(X0))] + #formulaNull <-update.formula(formula(obj),paste("~.-",paste(obj$nameUV$V, collapse="-"))) #togli tutti i termini "V" + X1<-matrix(Z, nrow=n, ncol=n1, byrow=FALSE) #matrice della variabile Z + fn="pmax(x-p,0)" + PSI<-matrix(values, nrow=n, ncol=n1, byrow=TRUE) + X<-eval(parse(text=fn), list(x=X1, p=PSI)) # fn t.c. length(fn)<=1 + mf$pmaxMedio<-rowMeans(X) + nc<-obj$orig.call #se c'e' un break e vuoi saggiare uno in piu' devi aggiustare la call + if(more.break){ + nc$formula<-update.formula(nc$formula, paste("~.+",paste(obj$nameUV$U, collapse="+"))) + } + formulaNull<-nc$formula + nc$data=quote(mf) + a<-eval(nc) + #assign("mf", mf, envir=sys.frame()) #funziona ma R CMD check da problemi.. + #ne <- new.env(parent = baseenv()) + pos<-1 + assign("mf", mf, envir=as.environment(pos)) + #r<-as.numeric(add1(a, ~.+pmaxMedio, scale=dispersion, test="Rao")[c("scaled Rao sc.", "Pr(>Chi)")][2,]) + r<- as.numeric(add1(a, ~.+pmaxMedio, scale=dispersion, test="Rao")[2,4:5]) + } else { + Call<-mf<-obj$call + mf$formula<-formula(obj) + m <- match(c("formula", "data", "subset", "weights", "na.action","offset"), names(mf), 0L) + mf <- mf[c(1, m)] + mf$drop.unused.levels <- TRUE + mf[[1L]] <- as.name("model.frame") + mf$formula<-update.formula(mf$formula,paste(seg.Z,collapse=".+")) + formulaNull <- formula(obj) + mf <- eval(mf, parent.frame()) + mt <- attr(mf, "terms") + XREG <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) + n <- nrow(XREG) + Z<- XREG[,match(name.Z, colnames(XREG))] + if(!name.Z %in% names(coef(obj))) XREG<-XREG[,-match(name.Z, colnames(XREG)),drop=FALSE] + if(is.null(values)) values<-seq(min(Z), max(Z), length=k) #values<-seq(sort(Z)[2], sort(Z)[(n - 1)], length = k) + n1<-length(values) + PSI<-matrix(values, nrow=n, ncol=n1, byrow=TRUE) #(era X2) matrice di valori di psi + X1<-matrix(Z, nrow=n, ncol=n1, byrow=FALSE) #matrice della variabile Z + fn="pmax(x-p,0)" + X<-eval(parse(text=fn), list(x=X1, p=PSI)) # fn t.c. length(fn)<=1 + pmaxMedio<-rowMeans(X) + #r<-as.numeric(add1(update(obj, data=mf), ~.+pmaxMedio, scale=dispersion, test="Rao")[c("scaled Rao sc.", "Pr(>Chi)")][2,]) + r<-as.numeric(add1(update(obj, data=mf), ~.+pmaxMedio, scale=dispersion, test="Rao")[2,4:5]) + } + } else { + #SE E' un LM + if(is.null(dispersion)) dispersion<- summary(obj)$sigma^2 + if(is.null(df.t)) df.t <- obj$df.residual + #df.ok<- if(!is.null(df.t)) df.t else obj$df.residual + # se e' LM-segmented + if(inherits(obj, "segmented")){ + if(!is.null(eval(obj$call$obj)$call$data)) mf$data <- eval(obj$call$obj)$call$data + y<- obj$res+obj$fitted + if(!is.null(obj$offset)) y<- y-obj$offset + weights<- obj$weights + X0<-model.matrix(obj) + Z<-X0[,name.Z] + #escludi *tutte* le variabili psi (sia da X0 che dalla formula) + X0<-X0[, -match(obj$nameUV$V, colnames(X0))] + formulaNull <-update.formula(formula(obj),paste("~.-",paste(obj$nameUV$V, collapse="-"))) #togli tutti i termini "V" + #se vuoi saggiare per un additional breakpoint + if(!more.break){ + idU<-grep(name.Z, obj$nameUV$U) + X0<-X0[, -match(obj$nameUV$U[idU], colnames(X0))] + formulaNull <-update.formula(formulaNull, paste("~.-",paste(obj$nameUV$U[idU], collapse="-"))) + } + #browser() + if(is.null(values)) values<-seq(min(Z), max(Z), length=k) #values<-seq(sort(Z)[2], sort(Z)[(n - 1)], length = k) +# formulaOrig<-mf$formula<-update.formula(mf$formula,paste("~.-",paste(obj$nameUV$V, collapse="-"))) #togli tutti i termini "V" +# update.formula(formulaOrig, paste("~.-",paste(obj$nameUV$U[idU], collapse="-"))) +# #for(i in 1:length(obj$nameUV$U)) assign(obj$nameUV$U[i], obj$model[,obj$nameUV$U[i]], envir=parent.frame()) +# formulaOrig<-update.formula(formulaOrig, paste("~.-",paste(obj$nameUV$U, collapse="-"))) + r<-test.Sc2(y=y, z=Z, xreg=X0, sigma=sqrt(dispersion), values=values, fn=fn, df.t=df.t, alternative=alternative, w=weights) + } else { + #SE l'oggetto obj NON e' "segmented".. + # ci possono essere mancanti nella variabile seg.Z, quindi devi costruire un nuovo dataframe.... + Call<-mf<-obj$call + mf$formula<-formula(obj) + m <- match(c("formula", "data", "subset", "weights", "na.action","offset"), names(mf), 0L) + mf <- mf[c(1, m)] + mf$drop.unused.levels <- TRUE + mf[[1L]] <- as.name("model.frame") + mf$formula<-update.formula(mf$formula,paste(seg.Z,collapse=".+")) + formulaNull <- formula(obj) + mf <- eval(mf, parent.frame()) + + weights <- as.vector(model.weights(mf)) + offs <- as.vector(model.offset(mf)) + + mt <- attr(mf, "terms") + interc<-attr(mt,"intercept") + y <- model.response(mf, "any") + XREG <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) + n <- nrow(XREG) + if (is.null(weights)) weights <- rep(1, n) + if (!is.null(offs)) y<-y-offs + Z<- XREG[,match(name.Z, colnames(XREG))] + if(!name.Z %in% names(coef(obj))) XREG<-XREG[,-match(name.Z, colnames(XREG)),drop=FALSE] + if(is.null(values)) values<-seq(min(Z), max(Z), length=k) #values<-seq(sort(Z)[2], sort(Z)[(n - 1)], length = k) + r<-test.Sc2(y=y, z=Z, xreg=XREG, sigma=sqrt(dispersion), values=values, fn=fn, df.t=df.t, alternative=alternative, w=weights) + } + } + if(is.null(obj$family$family)) { + famiglia<-"gaussian" + legame<-"identity" + } else { + famiglia<-obj$family$family + legame<-obj$family$link + } + + out <- list(method = "Score test for one change in the slope", +# data.name=paste("Model = ",famiglia,", link =", legame, +# "\nformula =", as.expression(formulaOrig), +# "\nsegmented variable =", name.Z), + data.name=paste("formula =", as.expression(formulaNull), ", method =", obj$call[[1]] , + "\nmodel =",famiglia,", link =", legame, NULL , + "\nsegmented variable =", name.Z), + statistic = c(`observed value` = r[1]), + parameter = c(n.points = length(values)), p.value = r[2], + alternative = alternative) + class(out) <- "htest" + return(out) +} + + + + + + + + diff --git a/R/seg.Ar.fit.r b/R/seg.Ar.fit.r index 7db6288..05319fc 100644 --- a/R/seg.Ar.fit.r +++ b/R/seg.Ar.fit.r @@ -10,6 +10,7 @@ dpmax<-function(x,y,pow=1){ c2 <- apply((Z >= PSI), 2, all) if(sum(c1 + c2) != 0 || is.na(sum(c1 + c2))) stop("psi out of the range") # + digits<-opz$digits pow<-opz$pow nomiOK<-opz$nomiOK toll<-opz$toll @@ -82,6 +83,8 @@ dpmax<-function(x,y,pow=1){ psi.values[[length(psi.values) + 1]] <- psi.old <- psi # if(it>=old.it.max && h<1) H<-h psi <- psi.old + h*gamma.c/beta.c + if(!is.null(digits)) + psi<-round(psi, digits) PSI <- matrix(rep(psi, rep(nrow(Z), length(psi))), ncol = length(psi)) #check if psi is admissible.. a <- apply((Z <= PSI), 2, all) #prima era solo < diff --git a/R/seg.control.R b/R/seg.control.R index 969d33f..03e4868 100644 --- a/R/seg.control.R +++ b/R/seg.control.R @@ -1,7 +1,7 @@ `seg.control` <- function(toll=.0001, it.max=10, display=FALSE, stop.if.error=TRUE, K=10, quant=FALSE, last=TRUE, maxit.glm=25, h=1, - n.boot=20, size.boot=NULL, gap=FALSE, jt=FALSE, nonParam=TRUE, random=TRUE, powers=c(1,1), seed=NULL, fn.obj=NULL){ + n.boot=20, size.boot=NULL, gap=FALSE, jt=FALSE, nonParam=TRUE, random=TRUE, powers=c(1,1), seed=NULL, fn.obj=NULL, digits=NULL){ list(toll=toll,it.max=it.max,visual=display,stop.if.error=stop.if.error, K=K,last=last,maxit.glm=maxit.glm,h=h,n.boot=n.boot, size.boot=size.boot, gap=gap, jt=jt, - nonParam=nonParam, random=random, pow=powers, seed=seed, quant=quant, fn.obj=fn.obj)} + nonParam=nonParam, random=random, pow=powers, seed=seed, quant=quant, fn.obj=fn.obj, digits=digits)} diff --git a/R/seg.def.fit.r b/R/seg.def.fit.r index 56a8ce6..fbc9a59 100644 --- a/R/seg.def.fit.r +++ b/R/seg.def.fit.r @@ -1,15 +1,38 @@ seg.def.fit<-function(obj, Z, PSI, mfExt, opz, return.all.sol=FALSE){ #----------------- +fn.costr<-function(n.psi,isLeft=1,isInterc=1){ +#build the constraint matrix +#isLeft: TRUE (or 1) if there is the left slope +#isInterc: TRUE (or 1) if there is the intercept.. +IU<- -diag(n.psi) +sumU<- diag(n.psi) #n. of diff slopes +sumU[row(sumU)>col(sumU)]<-1 +if(isLeft) { + sumU<-cbind(1, sumU) + IU<-diag(c(1, -rep(1, n.psi))) + } + A<-rbind(IU,sumU) + if(isInterc) { + A<-rbind(0,A) + A<-cbind(c(1, rep(0,nrow(A)-1)), A) + } + #add zeros for the V coeffs + A<-cbind(A, matrix(0,nrow(A), n.psi)) + A + } +#----------------- dpmax<-function(x,y,pow=1){ #deriv pmax if(pow==1) ifelse(x>y, -1, 0) else -pow*pmax(x-y,0)^(pow-1) } #----------- + vincoli<- FALSE c1 <- apply((Z <= PSI), 2, all) c2 <- apply((Z >= PSI), 2, all) if(sum(c1 + c2) != 0 || is.na(sum(c1 + c2))) stop("psi out of the range") # + digits<-opz$digits pow<-opz$pow nomiOK<-opz$nomiOK toll<-opz$toll @@ -28,11 +51,25 @@ dpmax<-function(x,y,pow=1){ epsilon <- 10 dev.values<-psi.values <- NULL id.psi.ok<-rep(TRUE, length(psi)) - + nomiU<- opz$nomiU nomiV<- opz$nomiV call.ok <- opz$call.ok call.noV <- opz$call.noV + +#browser() + if(is.null(opz$constr)) opz$constr<-0 + if((opz$constr %in% 1:2) && class(obj)=="rq"){ + vincoli<-TRUE + call.ok$method<-"fnc" + call.ok$R<-quote(R) + call.ok$r<-quote(r) + + call.noV$method<-"fnc" + call.noV$R<-quote(R.noV) + call.noV$r<-quote(r) + } + fn.obj<-opz$fn.obj toll<-opz$toll k<-ncol(Z) @@ -44,6 +81,9 @@ dpmax<-function(x,y,pow=1){ mfExt[nomiU[i]] <- U[,i] mfExt[nomiV[i]] <- V[,i] } + R<- fn.costr(ncol(U),1,1) + R.noV<-R[,-((ncol(R)-1)+seq_len(ncol(U))),drop=FALSE] + r<- rep(0, nrow(R)) obj <- suppressWarnings(eval(call.ok, envir=mfExt)) dev.old<-dev.new dev.new <- dev.new1 <- eval(parse(text=fn.obj), list(x=obj)) #control$f.obj should be something like "sum(x$residuals^2)" or "x$dev" @@ -74,6 +114,8 @@ dpmax<-function(x,y,pow=1){ psi.values[[length(psi.values) + 1]] <- psi.old <- psi # if(it>=old.it.max && h<1) H<-h psi <- psi.old + h*gamma.c/beta.c + if(!is.null(digits)) + psi<-round(psi, digits) PSI <- matrix(rep(psi, rep(nrow(Z), length(psi))), ncol = length(psi)) #check if psi is admissible.. a <- apply((Z <= PSI), 2, all) #prima era solo < @@ -96,8 +138,7 @@ dpmax<-function(x,y,pow=1){ nomiOK<-nomiOK[id.psi.ok] #salva i nomi delle U per i psi ammissibili nomiU<-nomiU[id.psi.ok] #salva i nomi delle U per i psi ammissibili nomiV<-nomiV[id.psi.ok] #salva i nomi delle V per i psi ammissibili - - + id.psi.group<-id.psi.group[id.psi.ok] names(psi)<-id.psi.group if(ncol(PSI)<=0) return(0) @@ -105,7 +146,8 @@ dpmax<-function(x,y,pow=1){ #aggiorna la call, altrimenti il modello avra' sempre lo stesso numero di termini anche se alcuni psi vengono rimossi!!! Fo <- update.formula(opz$formula.orig, as.formula(paste(".~.+", paste(c(nomiU, nomiV), collapse = "+")))) Fo.noV <- update.formula(opz$formula.orig, as.formula(paste(".~.+", paste(nomiU, collapse = "+")))) - + + call.ok <- update(obj, formula = Fo, evaluate=FALSE, data = mfExt) call.noV <- update(obj, formula = Fo.noV, evaluate=FALSE, data = mfExt) @@ -114,7 +156,6 @@ dpmax<-function(x,y,pow=1){ #obj$psi <- psi } #end while -#browser() psi<-unlist(tapply(psi, id.psi.group, sort)) names(psi)<-id.psi.group PSI <- matrix(rep(psi, rep(nrow(Z), length(psi))), ncol = length(psi)) @@ -151,7 +192,14 @@ dpmax<-function(x,y,pow=1){ obj$it <- it #fino a qua.. obj<-list(obj=obj,it=it,psi=psi,psi.values=psi.values,U=U,V=V,rangeZ=rangeZ, - epsilon=epsilon,nomiOK=nomiOK, SumSquares.no.gap=SS.new, id.psi.group=id.psi.group, nomiV=nomiV, nomiU=nomiU,mfExt=mfExt) #inserire id.psi.ok? - return(obj) + epsilon=epsilon,nomiOK=nomiOK, SumSquares.no.gap=SS.new, id.psi.group=id.psi.group, + nomiV=nomiV, nomiU=nomiU, mfExt=mfExt) #inserire id.psi.ok? + #browser() + if(vincoli) { + obj$R<-R + obj$R.noV<-R.noV + obj$r<-r + } + return(obj) } diff --git a/R/seg.glm.fit.r b/R/seg.glm.fit.r index 07f016c..abc5e14 100644 --- a/R/seg.glm.fit.r +++ b/R/seg.glm.fit.r @@ -9,6 +9,7 @@ dpmax<-function(x,y,pow=1){ c1 <- apply((Z <= PSI), 2, all) #prima era solo < c2 <- apply((Z >= PSI), 2, all) #prima era solo > if(sum(c1 + c2) != 0 || is.na(sum(c1 + c2))) stop("psi out of the range") + digits<-opz$digits pow<-opz$pow eta0<-opz$eta0 fam<-opz$fam @@ -83,6 +84,7 @@ dpmax<-function(x,y,pow=1){ psi.values[[length(psi.values) + 1]] <- psi.old <- psi #if(it>=old.it.max && h<1) H<-h psi <- psi.old + h*gamma.c/beta.c + if(!is.null(digits)) psi<-round(psi, digits) PSI <- matrix(rep(psi, rep(nrow(Z), ncol(Z))), ncol = ncol(Z)) #check if psi is admissible.. a <- apply((Z <= PSI), 2, all) #prima era solo < diff --git a/R/seg.lm.fit.r b/R/seg.lm.fit.r index 99d6a54..64550b6 100644 --- a/R/seg.lm.fit.r +++ b/R/seg.lm.fit.r @@ -26,6 +26,7 @@ mylm<-function(x,y,w,offs=rep(0,length(y))){ c2 <- apply((Z >= PSI), 2, all) if(sum(c1 + c2) != 0 || is.na(sum(c1 + c2))) stop("psi out of the range") # + digits<-opz$digits pow<-opz$pow nomiOK<-opz$nomiOK toll<-opz$toll @@ -59,7 +60,7 @@ mylm<-function(x,y,w,offs=rep(0,length(y))){ obj <- lm.wfit(x = X, y = y, w = w, offset = offs) dev.old<-dev.new dev.new <- dev.new1 <-sum(obj$residuals^2) - if(return.all.sol) dev.new1 <- sum(mylm(x = cbind(XREG, U), y = y, w = w, offs = offs)$residuals^2) + if(return.all.sol) dev.new1 <- sum(mylm(x = cbind(XREG, U), y = y, w = w, offs = offs)$residuals^2*w/sum(w)) dev.values[[length(dev.values) + 1]] <- dev.new1 if (visual) { flush.console() @@ -88,6 +89,7 @@ mylm<-function(x,y,w,offs=rep(0,length(y))){ psi.values[[length(psi.values) + 1]] <- psi.old <- psi # if(it>=old.it.max && h<1) H<-h psi <- psi.old + h*gamma.c/beta.c + if(!is.null(digits)) psi<-round(psi, digits) PSI <- matrix(rep(psi, rep(nrow(Z), length(psi))), ncol = length(psi)) #check if psi is admissible.. a <- apply((Z <= PSI), 2, all) #prima era solo < @@ -124,7 +126,7 @@ mylm<-function(x,y,w,offs=rep(0,length(y))){ obj$epsilon <- epsilon obj$it <- it obj.new <- lm.wfit(x = cbind(XREG, U), y = y, w = w, offset = offs) - SS.new<-sum(obj.new$residuals^2) + SS.new<-sum(obj.new$residuals^2*w/sum(w)) if(!gap){ names.coef<-names(obj$coefficients) obj$coefficients<-c(obj.new$coefficients, rep(0,ncol(V))) diff --git a/R/segmented.Arima.r b/R/segmented.Arima.r index 47670a8..9566f69 100644 --- a/R/segmented.Arima.r +++ b/R/segmented.Arima.r @@ -25,6 +25,7 @@ dpmax<-function(x,y,pow=1){ } if(length(all.vars(seg.Z))!=n.Seg) stop("A wrong number of terms in `seg.Z' or `psi'") it.max <- old.it.max<- control$it.max + digits<-control$digits toll <- control$toll visual <- control$visual stop.if.error<-control$stop.if.error @@ -141,7 +142,7 @@ dpmax<-function(x,y,pow=1){ nomiOK<-nomiU opz<-list(toll=toll,h=h,stop.if.error=stop.if.error,dev0=dev0,visual=visual,it.max=it.max, - nomiOK=nomiOK, id.psi.group=id.psi.group, gap=gap, visualBoot=visualBoot, pow=pow) + nomiOK=nomiOK, id.psi.group=id.psi.group, gap=gap, visualBoot=visualBoot, pow=pow, digits=digits) opz$call.ok<-call.ok opz$call.noV<-call.noV diff --git a/R/segmented.default.r b/R/segmented.default.r index 8a3ad4a..5ef7d95 100644 --- a/R/segmented.default.r +++ b/R/segmented.default.r @@ -33,6 +33,7 @@ dpmax<-function(x,y,pow=1){ } if(length(all.vars(seg.Z))!=n.Seg) stop("A wrong number of terms in `seg.Z' or `psi'") it.max <- old.it.max<- control$it.max + digits<-control$digits toll <- control$toll visual <- control$visual stop.if.error<-control$stop.if.error @@ -171,8 +172,9 @@ dpmax<-function(x,y,pow=1){ Fo <- update.formula(formula(obj), as.formula(paste(".~.+", paste(nnomi, collapse = "+")))) Fo.noV <- update.formula(formula(obj), as.formula(paste(".~.+", paste(nomiU, collapse = "+")))) - call.ok <- update(obj, formula = Fo, evaluate=FALSE, data = mfExt) #objF <- update(obj0, formula = Fo, data = KK) - call.noV <- update(obj, formula = Fo.noV, evaluate=FALSE, data = mfExt) #objF <- update(obj0, formula = Fo, data = KK) + #ho tolto "formula = Fo" + call.ok <- update(obj, Fo, evaluate=FALSE, data = mfExt) #objF <- update(obj0, formula = Fo, data = KK) + call.noV <- update(obj, Fo.noV, evaluate=FALSE, data = mfExt) #objF <- update(obj0, formula = Fo, data = KK) if (it.max == 0) { if(!is.null(call.noV[["subset"]])) call.noV[["subset"]]<-NULL @@ -195,7 +197,7 @@ dpmax<-function(x,y,pow=1){ nomiOK<-nomiU opz<-list(toll=toll,h=h,stop.if.error=stop.if.error,dev0=dev0,visual=visual,it.max=it.max, - nomiOK=nomiOK, id.psi.group=id.psi.group, gap=gap, visualBoot=visualBoot, pow=pow) + nomiOK=nomiOK, id.psi.group=id.psi.group, gap=gap, visualBoot=visualBoot, pow=pow, digits=digits) opz$call.ok<-call.ok opz$call.noV<-call.noV @@ -204,6 +206,10 @@ dpmax<-function(x,y,pow=1){ opz$nomiV<-nomiV opz$fn.obj <- fn.obj + opz<-c(opz,...) #8/10/16... + +#browser() + if(n.boot<=0){ obj<-seg.def.fit(obj, Z, PSI, mfExt, opz) } else { @@ -215,7 +221,7 @@ dpmax<-function(x,y,pow=1){ warning("No breakpoint estimated", call. = FALSE) return(obj0) } - if(!is.null(obj$obj$df.residual)){ + if(!is.null(obj$obj$df.residual) && !is.na(obj$obj$df.residual)){ if(obj$obj$df.residual==0) warning("no residual degrees of freedom (other warnings expected)", call.=FALSE) } id.psi.group<-obj$id.psi.group @@ -234,6 +240,8 @@ dpmax<-function(x,y,pow=1){ V<-obj$V # return(obj) +#browser() + #if(any(table(rowSums(V))<=1)) stop("only 1 datum in an interval: breakpoint(s) at the boundary or too close") for(jj in colnames(V)) { VV<-V[, which(colnames(V)==jj), drop=FALSE] @@ -245,10 +253,13 @@ dpmax<-function(x,y,pow=1){ rangeZ<-obj$rangeZ mfExt<-obj$mfExt names(mfExt)[match(obj$nomiV, names(mfExt))]<-nomiVxb + R<-obj$R + R.noV<-obj$R.noV + r<-obj$r obj<-obj$obj k<-length(psi) # beta.c<-if(k == 1) coef(obj)["U"] else coef(obj)[paste("U", 1:ncol(U), sep = "")] -#browser() + beta.c<- coef(obj)[nomiOK] #nomiOK e' stato estratto da obj e contiene tutti i nomi delle variabili U inserite nel modello psi.values[[length(psi.values) + 1]] <- psi id.warn <- FALSE @@ -258,7 +269,8 @@ dpmax<-function(x,y,pow=1){ } Vxb <- V %*% diag(beta.c, ncol = length(beta.c)) - #se usi una procedura automatica devi cambiare ripetizioni, nomiU e nomiV, e quindi: + +# #se usi una procedura automatica devi cambiare ripetizioni, nomiU e nomiV, e quindi: # length.psi<-tapply(as.numeric(as.character(names(psi))), as.numeric(as.character(names(psi))), length) # forma.nomiU<-function(xx,yy)paste("U",1:xx, ".", yy, sep="") # forma.nomiVxb<-function(xx,yy)paste("psi",1:xx, ".", yy, sep="") @@ -278,11 +290,23 @@ dpmax<-function(x,y,pow=1){ # mfExt[nomiOK[i]]<-mf[nomiOK[i]]<-U[,i] # mfExt[nomiVxb[i]]<-mf[nomiVxb[i]]<-Vxb[,i] # } +#ottobre 2016: ci vuole altrimenti Vxb non viene inserita nel dataframe + for(i in 1:ncol(U)) { + mfExt[nomiU[i]]<-mf[nomiU[i]]<-U[,i] + mfExt[nomiVxb[i]]<-mf[nomiVxb[i]]<-Vxb[,i] + } nnomi <- c(nomiOK, nomiVxb) Fo <- update.formula(formula(obj0), as.formula(paste(".~.+", paste(nnomi, collapse = "+")))) - objF <- update(obj0, formula = Fo, evaluate=FALSE, data = mfExt) + objF <- update(obj0, Fo, evaluate=FALSE, data = mfExt) #tolto "formula =Fo" if(!is.null(objF[["subset"]])) objF[["subset"]]<-NULL + + if(is.null(opz$constr)) opz$constr<-0 + if((opz$constr %in% 1:2) && class(obj0)=="rq"){ + objF$method<-"fnc" + objF$R<-quote(R) + objF$r<-quote(r) + } objF<- eval(objF, envir=mfExt) #Puo' capitare che psi sia ai margini e ci sono 1 o 2 osservazioni in qualche intervallo. Oppure ce ne #sono di piu' ma hanno gli stessi valori di x @@ -294,7 +318,10 @@ dpmax<-function(x,y,pow=1){ call. = FALSE)} else { warning("some estimate is NA: too many breakpoints? 'var(hat.psi)' cannot be computed \n ..returning a 'lm' model", call. = FALSE) Fo <- update.formula(formula(obj0), as.formula(paste(".~.+", paste(nomiU, collapse = "+")))) - objF <- update(obj0, formula = Fo, evaluate=TRUE, data = mfExt) + objF <- if((opz$constr %in% 1:2) && class(obj0)=="rq") {update(obj0, formula = Fo, R=R.noV, r=r, method="fnc", evaluate=TRUE, data = mfExt) + } else { + update(obj0, Fo, evaluate=TRUE, data = mfExt) #ho tolto "formula = Fo" per farlo funzionare con gls() + } names(psi)<-nomiVxb objF$psi<-psi return(objF) diff --git a/R/segmented.glm.R b/R/segmented.glm.R index b428fdb..0274cc3 100644 --- a/R/segmented.glm.R +++ b/R/segmented.glm.R @@ -1,6 +1,7 @@ `segmented.glm` <- function(obj, seg.Z, psi, control = seg.control(), model = TRUE, ...) { n.Seg<-1 + if(missing(seg.Z) && length(all.vars(formula(obj)))==2) seg.Z<- as.formula(paste("~", all.vars(formula(obj))[2])) if(missing(psi)){if(length(all.vars(seg.Z))>1) stop("provide psi") else psi<-Inf} if(length(all.vars(seg.Z))>1 & !is.list(psi)) stop("`psi' should be a list with more than one covariate in `seg.Z'") if(is.list(psi)){ @@ -11,6 +12,7 @@ function(obj, seg.Z, psi, control = seg.control(), model = TRUE, ...) { if(length(all.vars(seg.Z))!=n.Seg) stop("A wrong number of terms in `seg.Z' or `psi'") maxit.glm <- control$maxit.glm it.max <- old.it.max<- control$it.max + digits<-control$digits toll <- control$toll visual <- control$visual stop.if.error<-control$stop.if.error @@ -240,7 +242,7 @@ function(obj, seg.Z, psi, control = seg.control(), model = TRUE, ...) { nomiOK<-nomiU opz<-list(toll=toll,h=h,stop.if.error=stop.if.error,dev0=dev0,visual=visual,it.max=it.max,nomiOK=nomiOK, fam=fam, eta0=obj$linear.predictors, maxit.glm=maxit.glm, id.psi.group=id.psi.group, gap=gap, - pow=pow, visualBoot=visualBoot) + pow=pow, visualBoot=visualBoot, digits=digits) if(n.boot<=0){ obj<-seg.glm.fit(y,XREG,Z,PSI,weights,offs,opz) diff --git a/R/segmented.lm.R b/R/segmented.lm.R index edd849d..dcf90be 100644 --- a/R/segmented.lm.R +++ b/R/segmented.lm.R @@ -1,6 +1,7 @@ `segmented.lm` <- function(obj, seg.Z, psi, control = seg.control(), model = TRUE, ...) { n.Seg<-1 + if(missing(seg.Z) && length(all.vars(formula(obj)))==2) seg.Z<- as.formula(paste("~", all.vars(formula(obj))[2])) if(missing(psi)){if(length(all.vars(seg.Z))>1) stop("provide psi") else psi<-Inf} if(length(all.vars(seg.Z))>1 & !is.list(psi)) stop("`psi' should be a list with more than one covariate in `seg.Z'") if(is.list(psi)){ @@ -10,6 +11,7 @@ function(obj, seg.Z, psi, control = seg.control(), model = TRUE, ...) { } if(length(all.vars(seg.Z))!=n.Seg) stop("A wrong number of terms in `seg.Z' or `psi'") it.max <- old.it.max<- control$it.max + digits<-control$digits toll <- control$toll visual <- control$visual stop.if.error<-control$stop.if.error @@ -28,7 +30,7 @@ function(obj, seg.Z, psi, control = seg.control(), model = TRUE, ...) { set.seed(employed.Random.seed) } if(visual) {visual<-FALSE; visualBoot<-TRUE}# warning("`display' set to FALSE with bootstrap restart", call.=FALSE)} - if(!stop.if.error) stop("Bootstrap restart only with a fixed number of breakpoints") +# if(!stop.if.error) stop("Bootstrap restart only with a fixed number of breakpoints") } last <- control$last K<-control$K @@ -251,7 +253,7 @@ if(!is.null(nomiNO)) mfExt$formula<-update.formula(mfExt$formula,paste(".~.-", p # psi.values <- NULL nomiOK<-nomiU opz<-list(toll=toll,h=h,stop.if.error=stop.if.error,dev0=dev0,visual=visual,it.max=it.max, - nomiOK=nomiOK,id.psi.group=id.psi.group,gap=gap,visualBoot=visualBoot,pow=pow) + nomiOK=nomiOK,id.psi.group=id.psi.group,gap=gap,visualBoot=visualBoot,pow=pow,digits=digits) if(n.boot<=0){ obj<-seg.lm.fit(y,XREG,Z,PSI,weights,offs,opz) } else { @@ -275,19 +277,7 @@ if(!is.null(nomiNO)) mfExt$formula<-update.formula(mfExt$formula,paste(".~.-", p psi.values<-if(n.boot<=0) obj$psi.values else obj$boot.restart U<-obj$U V<-obj$V -#browser() - - #if(any(table(rowSums(V))<=1)) stop("only 1 datum in an interval: breakpoint(s) at the boundary or too close") - for(jj in colnames(V)) { - VV<-V[, which(colnames(V)==jj), drop=FALSE] - sumV<-abs(rowSums(VV)) -# if( (any(diff(sumV)>=2) #se ci sono due breakpoints uguali -# || any(table(sumV)<=1)) && stop.if.error) stop("only 1 datum in an interval: breakpoint(s) at the boundary or too close each other") -# Tolto perche' se la variabile segmented non e' ordinata non ha senso.. -#magari potresti fare un abs(diff(psi))<=.0001? ma clusterizzato.. - if(any(table(sumV)<=1) && stop.if.error) stop("only 1 datum in an interval: breakpoint(s) at the boundary or too close each other") - } rangeZ<-obj$rangeZ obj<-obj$obj k<-length(psi) @@ -336,16 +326,31 @@ if(!is.null(nomiNO)) mfExt$formula<-update.formula(mfExt$formula,paste(".~.-", p #eliminiamo subset, perche' se e' del tipo subset=x>min(x) allora continuerebbe a togliere 1 osservazione if(!is.null(objF[["subset"]])) objF[["subset"]]<-NULL objF<-eval(objF, envir=mfExt) + + +# #11/10/16 il controllo e' stato commentato in modo tale da restituire anche un oggetto lm in cui psi viene considerato fisso.. +# for(jj in colnames(V)) { +# VV<-V[, which(colnames(V)==jj), drop=FALSE] +# sumV<-abs(rowSums(VV)) +## if( (any(diff(sumV)>=2) #se ci sono due breakpoints uguali +## || any(table(sumV)<=1)) && stop.if.error) stop("only 1 datum in an interval: breakpoint(s) at the boundary or too close each other") +## Tolto perche' se la variabile segmented non e' ordinata non ha senso.. +##magari potresti fare un abs(diff(psi))<=.0001? ma clusterizzato.. +# if(any(table(sumV)<=1) && stop.if.error) stop("only 1 datum in an interval: breakpoint(s) at the boundary or too close each other") +# } + + #Puo' capitare che psi sia ai margini o molto vicini (e ci sono solo 1 o 2 osservazioni in qualche intervallo. Oppure ce ne #sono di piu' ma hanno gli stessi valori di x. In questo caso objF$coef puo' avere mancanti.. names(which(is.na(coef(objF)))) - + objF$offset<- obj0$offset isNAcoef<-any(is.na(objF$coefficients)) if(isNAcoef){ if(stop.if.error) {stop("at least one coef is NA: breakpoint(s) at the boundary? (possibly with many x-values replicated)", - call. = FALSE)} else { + call. = FALSE) + } else { warning("some estimate is NA: too many breakpoints? 'var(hat.psi)' cannot be computed \n ..returning a 'lm' model", call. = FALSE) Fo <- update.formula(formula(obj0), as.formula(paste(".~.+", paste(nomiU, collapse = "+")))) objF <- update(obj0, formula = Fo, evaluate=TRUE, data = mfExt) @@ -411,6 +416,13 @@ if(!is.null(nomiNO)) mfExt$formula<-update.formula(mfExt$formula,paste(".~.-", p objF$id.warn <- id.warn objF$orig.call<-orig.call if (model) objF$model <- mf #objF$mframe <- data.frame(as.list(KK)) + +# PSI <- matrix(rep(psi, rep(nrow(Z), length(psi))), ncol = length(psi)) +# SE.PSI <- matrix(rep(sqrt(vv), rep(nrow(Z), length(psi))), ncol = length(psi)) +# X.is<-model.matrix(Fo, data=objF$model) +# X.is[,nomiVxb]<-pnorm((Z-PSI)/SE.PSI)%*% diag(-beta.c, ncol = length(beta.c)) +# objF$cov.unscaled.is<-crossprod(X.is) +#browser() if(n.boot>0) objF$seed<-employed.Random.seed class(objF) <- c("segmented", class(obj0)) list.obj[[length(list.obj) + 1]] <- objF diff --git a/R/slope.R b/R/slope.R index c64f309..1628ca1 100644 --- a/R/slope.R +++ b/R/slope.R @@ -2,6 +2,7 @@ function(ogg, parm, conf.level=0.95, rev.sgn=FALSE, var.diff=FALSE, APC=FALSE, digits = max(3, getOption("digits") - 3)){ #-- + f.U<-function(nomiU, term=NULL){ #trasforma i nomi dei coeff U (o V) nei nomi delle variabili corrispondenti #and if 'term' is provided (i.e. it differs from NULL) the index of nomiU matching term are returned @@ -17,11 +18,23 @@ function(ogg, parm, conf.level=0.95, rev.sgn=FALSE, var.diff=FALSE, APC=FALSE, return(nomiU.ok) } #-- - if(!"segmented"%in%class(ogg)) stop("A 'segmented' model is requested") +# if(!"segmented"%in%class(ogg)) stop("A 'segmented' model is requested") if(var.diff && length(ogg$nameUV$Z)>1) { var.diff<-FALSE warning("var.diff set to FALSE with multiple segmented variables", call.=FALSE) } + #se e' un "newsegmented" + +# if(!is.null(ogg$R.slope)) { +# covv<-old.coef.var(ogg) +# ogg$coefficients<- covv$b +# covv<- covv$cov +# ogg$psi<-old.psi(ogg) +# ogg$nameUV<-old.nomi(ogg) +# } else { + covv<-try(vcov(ogg,var.diff=var.diff), silent=TRUE) +# } + nomepsi<-rownames(ogg$psi) #OK nomeU<-ogg$nameUV$U nomeZ<-ogg$nameUV$Z @@ -38,18 +51,29 @@ function(ogg, parm, conf.level=0.95, rev.sgn=FALSE, var.diff=FALSE, APC=FALSE, nomi<-nomi[-match(nomepsi,nomi)] #escludi i coef delle V index<-vector(mode = "list", length = length(nomeZ)) for(i in 1:length(nomeZ)) { - #id.cof.U<-grep(paste("\\.",nomeZ[i],"$",sep=""), nomi, value=FALSE) - #psii<-ogg$psi[grep(paste("\\.",nomeZ[i],"$",sep=""), rownames(ogg$psi), value=FALSE),2] - #id.cof.U<- match(grep(nomeZ[i], ogg$nameUV$U, value=TRUE), nomi) - #psii<-ogg$psi[grep(nomeZ[i], ogg$nameUV$V, value=TRUE),2] - #il paste con "$" (paste("\\.",nomeZ[i],"$",sep="")) e' utile perche' serve a distinguere variabili con nomi simili (ad es., "x" e "xx") - #Comunque nella versione dopo la 0.3-1.0 ho (FINALMENTE) risolto mettendo f.U - id.cof.U<- f.U(ogg$nameUV$U, nomeZ[i]) - #id.cof.U e' la posizione nel vettore ogg$nameUV$U; la seguente corregge per eventuali variabili che ci sono prima (ad es., interc) - id.cof.U<- id.cof.U + (match(ogg$nameUV$U[1], nomi)-1) - psii<- ogg$psi[f.U(ogg$nameUV$V, nomeZ[i]) , "Est."] - id.cof.U <- id.cof.U[order(psii)] +#---> DA RIMUOVERE E SOSTITUIRE CON QUELLI DI SUBITO DOPO? +# #id.cof.U<-grep(paste("\\.",nomeZ[i],"$",sep=""), nomi, value=FALSE) +# #psii<-ogg$psi[grep(paste("\\.",nomeZ[i],"$",sep=""), rownames(ogg$psi), value=FALSE),2] +# #id.cof.U<- match(grep(nomeZ[i], ogg$nameUV$U, value=TRUE), nomi) +# #psii<-ogg$psi[grep(nomeZ[i], ogg$nameUV$V, value=TRUE),2] +# #il paste con "$" (paste("\\.",nomeZ[i],"$",sep="")) e' utile perche' serve a distinguere variabili con nomi simili (ad es., "x" e "xx") +# #Comunque nella versione dopo la 0.3-1.0 ho (FINALMENTE) risolto mettendo f.U +# id.cof.U<- f.U(ogg$nameUV$U, nomeZ[i]) +# #id.cof.U e' la posizione nel vettore ogg$nameUV$U; la seguente corregge per eventuali variabili che ci sono prima (ad es., interc) +# id.cof.U<- id.cof.U + (match(ogg$nameUV$U[1], nomi)-1) +# psii<- ogg$psi[f.U(ogg$nameUV$V, nomeZ[i]) , "Est."] +# id.cof.U <- id.cof.U[order(psii)] +#---> + + #questi funzionano anche con oggetti da segreg + nomiPsi<-grep(paste(".", nomeZ[i], sep="") , ogg$nameUV$V, value=TRUE) + psii<- ogg$psi[nomiPsi , "Est."] + nomiU<-grep(paste(".", nomeZ[i], sep="") , ogg$nameUV$U, value=TRUE) + #cof<-coef(ogg)[nomiU] + id.cof.U<- match(nomiU, names(ogg$coefficients)) index[[i]]<-c(match(nomeZ[i],nomi), id.cof.U) + + } Ris<-list() digits <- max(3, getOption("digits") - 3) @@ -60,6 +84,7 @@ function(ogg, parm, conf.level=0.95, rev.sgn=FALSE, var.diff=FALSE, APC=FALSE, # if(transf=="APC") transf<-c("100*(exp(x)-1)", "100*exp(x)") # my.f<-function(x)eval(parse(text=transf[1])) # my.f.deriv<-function(x)eval(parse(text=transf[2])) + for(i in 1:length(index)){ ind<-as.numeric(na.omit(unlist(index[[i]]))) @@ -68,10 +93,10 @@ function(ogg, parm, conf.level=0.95, rev.sgn=FALSE, var.diff=FALSE, APC=FALSE, cof<-coef(ogg)[ind] cof.out<-M%*%cof - covv<-try(vcov(ogg,var.diff=var.diff), silent=TRUE) + if(class(covv)!="try-error"){ - covv<-covv[ind,ind] - cov.out<-M%*%covv%*%t(M) + cov.ok<-covv[ind,ind] + cov.out<-M%*%cov.ok%*%t(M) se.out<-sqrt(diag(cov.out)) k<-if("lm"%in%class(ogg)) abs(qt((1-conf.level)/2,df=ogg$df.residual)) else abs(qnorm((1-conf.level)/2)) k<-k*se.out diff --git a/R/vcov.segmented.R b/R/vcov.segmented.R index 75d35e6..353a3f2 100644 --- a/R/vcov.segmented.R +++ b/R/vcov.segmented.R @@ -16,7 +16,7 @@ vcov.segmented<-function (object, var.diff=FALSE, ...){ v<-summary.segmented(object, var.diff=TRUE, correlation = FALSE, ...)$cov.var.diff } else { so<-summary.segmented(object, var.diff=FALSE, correlation = FALSE, ...) - v<-so$sigma^2 * so$cov.unscaled + v<-so$sigma^2 * so$cov.unscaled #object$cov.unscaled.is } } return(v) diff --git a/data/down.rda b/data/down.rda index a735e71..b44d875 100644 Binary files a/data/down.rda and b/data/down.rda differ diff --git a/data/plant.rda b/data/plant.rda index e42e721..f492797 100644 Binary files a/data/plant.rda and b/data/plant.rda differ diff --git a/data/stagnant.rda b/data/stagnant.rda index 0abed91..e698515 100644 Binary files a/data/stagnant.rda and b/data/stagnant.rda differ diff --git a/inst/CITATION b/inst/CITATION index 7f21b38..e7f89e1 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -21,10 +21,10 @@ citEntry(entry="Article", volume = "8", number = "1", pages = "20--25", - url = "http://cran.r-project.org/doc/Rnews/";, + url = "https://cran.r-project.org/doc/Rnews/";, textVersion = paste("Vito M. R. Muggeo (2008).", "segmented: an R Package to Fit Regression Models with Broken-Line Relationships.", "R News, 8/1, 20-25.", - "URL http://cran.r-project.org/doc/Rnews/.";) + "URL https://cran.r-project.org/doc/Rnews/.";) ) diff --git a/man/broken.line.Rd b/man/broken.line.Rd index abb7e01..f92094d 100644 --- a/man/broken.line.Rd +++ b/man/broken.line.Rd @@ -44,7 +44,7 @@ y<-rpois(100,exp(2+1.8*pmax(z-.6,0))) o<-glm(y~z,family=poisson) o.seg<-segmented(o,seg.Z=~z,psi=.5) \dontrun{plot(z,y)} -\dontrun{points(z,broken.line(o.seg,link=FALSE)$fit,col=2,pch=20)} +\dontrun{points(z,broken.line(o.seg,link=FALSE)$fit,col=2)} #just to illustrate, use plot.segmented } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. diff --git a/man/davies.test.Rd b/man/davies.test.Rd index 1d3b4a4..f0f012c 100644 --- a/man/davies.test.Rd +++ b/man/davies.test.Rd @@ -14,8 +14,8 @@ davies.test(obj, seg.Z, k = 10, alternative = c("two.sided", "less", "greater"), \item{obj}{ a fitted model typically returned by \code{glm} or \code{lm}. Even an object returned by \code{segmented} can be set (e.g. if interest lies in testing for an additional breakpoint).} \item{seg.Z}{ a formula with no response variable, such as \code{seg.Z=~x1}, indicating the - (continuous) segmented variable being tested. Only a single variable may be tested and a - warning is printed when \code{seg.Z} includes two or more terms. } + (continuous) segmented variable being tested. Only a single variable may be tested and an + error is printed when \code{seg.Z} includes two or more terms. } \item{k}{ number of points where the test should be evaluated. See Details. } \item{alternative}{ a character string specifying the alternative hypothesis. } \item{type}{ the test statistic to be used (only for GLM, default to lrt. @@ -87,6 +87,7 @@ Davies, R.B. (2002) Hypothesis testing when a nuisance parameter is present only %%\section{Warning }{Currently \code{davies.test} does not work if the fitted model \code{ogg} %% does not include the segmented variable \code{term} being tested.} +\seealso{See also \code{\link{pscore.test}} which is more power, especially when the signal-to-noise ratio is low. } \examples{ \dontrun{ diff --git a/man/intercept.Rd b/man/intercept.Rd index 80d3e73..391cac4 100644 --- a/man/intercept.Rd +++ b/man/intercept.Rd @@ -8,7 +8,7 @@ Computes the intercepts of each `segmented' relationship in the fitted model. } \usage{ -intercept(ogg, parm, gap = TRUE, rev.sgn = FALSE, var.diff=FALSE, +intercept(ogg, parm, rev.sgn = FALSE, var.diff=FALSE, digits = max(3, getOption("digits") - 3)) } %- maybe also 'usage' for other objects documented here. @@ -19,9 +19,9 @@ intercept(ogg, parm, gap = TRUE, rev.sgn = FALSE, var.diff=FALSE, \item{parm}{ the segmented variable whose intercepts have to be computed. If missing all the segmented variables in the model are considered. } - \item{gap}{ - logical. should the intercepts account for the (possible) gaps? -} +% \item{gap}{ +% logical. should the intercepts account for the (possible) gaps? +%} \item{rev.sgn}{vector of logicals. The length should be equal to the length of \code{parm}, but it is recycled otherwise. When \code{TRUE} it is assumed that the current \code{parm} is `minus' the actual segmented variable, therefore the order is reversed before printing. This is useful when a null-constraint has been set on the last slope. diff --git a/man/plot.segmented.Rd b/man/plot.segmented.Rd index f0df8cb..202fddd 100644 --- a/man/plot.segmented.Rd +++ b/man/plot.segmented.Rd @@ -9,7 +9,7 @@ \usage{ \method{plot}{segmented}(x, term, add=FALSE, res=FALSE, conf.level=0, interc=TRUE, link=TRUE, res.col=1, rev.sgn=FALSE, const=0, shade=FALSE, rug=TRUE, - dens.rug=FALSE, dens.col = grey(0.8), show.gap=FALSE, transf=I, ...) + dens.rug=FALSE, dens.col = grey(0.8), transf=I, ...) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -34,9 +34,9 @@ at the foot of the plot.} \item{dens.rug}{when \code{TRUE} then smooth covariate distribution is plotted on the x-axis.} \item{dens.col}{if \code{dens.rug=TRUE}, it means the colour to be used to plot the density.} - \item{show.gap}{ when \code{FALSE} the (possible) gaps between the fitted lines at the estimated breakpoints - are hidden. When bootstrap restarting has been employed (default in \code{segmented}), \code{show.gap} is meaningless - as the gap coefficients are always set to zero in the fitted model.} +% \item{show.gap}{ when \code{FALSE} the (possible) gaps between the fitted lines at the estimated breakpoints +% are hidden. When bootstrap restarting has been employed (default in \code{segmented}), \code{show.gap} is meaningless +% as the gap coefficients are always set to zero in the fitted model.} \item{transf}{ A possible function to convert the fitted values before plotting. It is only effective if the fitted values refer to a linear or a generalized linear model (on the link scale) \emph{and} \code{res=FALSE}.} \item{\dots}{ other graphics parameters to pass to plotting commands: `col', `lwd' and `lty' (that @@ -47,13 +47,15 @@ \details{ Produces (or adds to the current device) the fitted segmented relationship between the response and the selected \code{term}. If the fitted model includes just a single `segmented' variable, - \code{term} may be omitted. Due to the parameterization of the segmented terms, sometimes - the fitted lines may not appear to join at the estimated breakpoints. If this is the case, the apparent - `gap' would indicate some lack-of-fit. However, since version 0.2-9.0, the gap coefficients are set to zero by default - (see argument \code{gap} in in \code{\link{seg.control}}). + \code{term} may be omitted. + %Due to the parameterization of the segmented terms, sometimes + %the fitted lines may not appear to join at the estimated breakpoints. If this is the case, the apparent + %`gap' would indicate some lack-of-fit. However, since version 0.2-9.0, the gap coefficients are set to zero by default + %(see argument \code{gap} in in \code{\link{seg.control}}). The partial residuals are computed as `fitted + residuals', where `fitted' are the fitted values of the - segmented relationship. Notice that for GLMs the residuals are the response residuals if \code{link=FALSE} and - the working residuals weighted by the IWLS weights if \code{link=TRUE}. + segmented relationship relevant to the covariate specified in \code{term}. + Notice that for GLMs the residuals are the response residuals if \code{link=FALSE} and + the working residuals if \code{link=TRUE}. %weighted by the IWLS weights [fino alla versione 0.5-2.0 i workRes were weighted by the IWLS weights] } \value{ None. @@ -74,13 +76,16 @@ set.seed(1234) z<-runif(100) y<-rpois(100,exp(2+1.8*pmax(z-.6,0))) o<-glm(y~z,family=poisson) -o.seg<-segmented(o,seg.Z=~z,psi=list(z=.5)) +o.seg<-segmented(o, ~z) #single segmented covariate and one breakpoint:'psi' can be omitted par(mfrow=c(2,1)) plot(o.seg, conf.level=0.95, shade=TRUE) +points(o.seg, link=FALSE, col=2) +## new plot plot(z,y) ## add the fitted lines using different colors and styles.. plot(o.seg,add=TRUE,link=FALSE,lwd=2,col=2:3, lty=c(1,3)) lines(o.seg,col=2,pch=19,bottom=FALSE,lwd=2) +points(o.seg,col=4, link=FALSE) } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. diff --git a/man/points.segmented.Rd b/man/points.segmented.Rd index 51b7fe8..6dcf1a2 100644 --- a/man/points.segmented.Rd +++ b/man/points.segmented.Rd @@ -66,8 +66,7 @@ joinpoints on the current plot. This could be useful to emphasize the changes of \examples{ \dontrun{ -#continues from ?plot.segmented -points(o.seg,col=2) +#see examples in ?plot.segmented } } \keyword{ nonlinear } diff --git a/man/pscore.test.rd b/man/pscore.test.rd new file mode 100644 index 0000000..12860e0 --- /dev/null +++ b/man/pscore.test.rd @@ -0,0 +1,82 @@ +\name{pscore.test} +\alias{pscore.test} +\title{ Testing for existence of one breakpoint} +\description{ + Given a (generalized) linear model, the (pseudo) Score statistic tests for the existence of one breakpoint. +} +\usage{ +pscore.test(obj, seg.Z, k = 10, alternative = c("two.sided", "less", "greater"), + values=NULL, dispersion=NULL, df.t=NULL, more.break=FALSE) +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{obj}{ a fitted model typically returned by \code{glm} or \code{lm}. Even an object returned + by \code{segmented} can be set. Offset and weights are allowed.} + \item{seg.Z}{ a formula with no response variable, such as \code{seg.Z=~x1}, indicating the + (continuous) segmented variable being tested. Only a single variable may be tested and an + error is printed when \code{seg.Z} includes two or more terms. \code{seg.Z} can be omitted if i)\code{obj} is a segmented fit with a single segmented covariate (and that variable is taken), or ii)if it is a "lm" or "glm" fit with a single covariate (and that variable is taken).} + \item{k}{ optional. Number of points used to compute the pseudo Score statistic. See Details. } + \item{alternative}{ a character string specifying the alternative hypothesis. } + \item{values}{ optional. The evaluation points where the Score test is computed. See Details for default values.} + \item{dispersion}{ optional. the dispersion parameter for the family to be used to compute the test statistic. + When \code{NULL} (the default), it is inferred from \code{obj}. Namely it is taken as \code{1} for the + Binomial and Poisson families, and otherwise estimated by the residual Chi-squared statistic in the model \code{obj} (calculated from cases with + non-zero weights divided by the residual degrees of freedom).} + \item{df.t}{ optional. The degress-of-freedom used to compute the p-value. When \code{NULL}, the df extracted from \code{obj} are used.} + \item{more.break}{ optional logical. If \code{obj} is a segmented fit, \code{more.break=FALSE} tests for the actual breakpoint in \code{obj}, + while \code{more.break=TRUE} tests for an \emph{additional} breakpoint. Ignored when \code{obj} is not a segmented fit.} + +} +\details{ + \code{pscore.test} tests for a non-zero difference-in-slope parameter of a segmented + relationship. Namely, the null hypothesis is \eqn{H_0:\beta=0}{H_0:beta=0}, where \eqn{\beta}{beta} is the difference-in-slopes, + i.e. the coefficient of the segmented function \eqn{\beta(x-\psi)_+}{beta*(x-psi)_+}. The hypothesis of interest + \eqn{\beta=0}{beta=0} means no breakpoint. Simulation studies have shown that such Score test is more powerful than the Davies test (see reference) when the alternative hypothesis is `one changepoint'. + + The \code{dispersion} value, if unspecified, is taken from \code{obj}. If \code{obj} represents the fit under the null hypothesis (no changepoint), the dispersion parameter estimate will be usually larger, leading to a (potentially severe) loss of power. + + The \code{k} evaluation points are \code{k} equally spaced values in the range of the segmented covariate. \code{k} should not be small. + Specific values can be set via \code{values}. However I have found no important difference due to number and location of the evaluation points, thus default is \code{k=10} equally-spaced points. +} +\value{ + A list with class '\code{htest}' containing the following components: + \item{method}{title (character)} + \item{data.name}{the regression model and the segmented variable being tested} + \item{statistic }{the point within the range of the covariate in \code{seg.Z} at which the maximum + (or the minimum if \code{alternative="less"}) occurs} + \item{parameter }{number of evaluation points} + \item{p.value }{the p-value} + \item{process}{the alternative hypothesis set} +} +\references{ +Muggeo, V.M.R. (2016) Testing with a nuisance parameter present only under the alternative: + a score-based approach with application to segmented modelling. + \emph{J of Statistical Computation and Simulation}, \bold{86}, 3059--3067. + } +\author{ Vito M.R. Muggeo } + +\seealso{See also \code{\link{davies.test}}. } + +\examples{ +\dontrun{ +set.seed(20) +z<-runif(100) +x<-rnorm(100,2) +y<-2+10*pmax(z-.5,0)+rnorm(100,0,3) + +o<-lm(y~z+x) + +#testing for one changepoint +#use the simple null fit +pscore.test(o,~z) #compare with davies.test(o,~z).. + +#use the segmented fit +os<-segmented(o, ~z) +pscore.test(os,~z) #smaller p-value, as it uses the dispersion under the alternative (from 'os') + +#test for the 2nd breakpoint in the variable z +pscore.test(os,~z, more.break=TRUE) + + } +} +\keyword{ htest } diff --git a/man/seg.control.Rd b/man/seg.control.Rd index 8f6690c..60014ce 100644 --- a/man/seg.control.Rd +++ b/man/seg.control.Rd @@ -10,7 +10,7 @@ seg.control(toll = 1e-04, it.max = 10, display = FALSE, stop.if.error = TRUE, K = 10, quant = FALSE, last = TRUE, maxit.glm = 25, h = 1, n.boot=20, size.boot=NULL, gap=FALSE, jt=FALSE, nonParam=TRUE, - random=TRUE, powers=c(1,1), seed=NULL, fn.obj=NULL) + random=TRUE, powers=c(1,1), seed=NULL, fn.obj=NULL, digits=NULL) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -70,6 +70,7 @@ seg.control(toll = 1e-04, it.max = 10, display = FALSE, stop.if.error = TRUE, \code{"coxph"} fits it should be \code{fn.obj="-x$loglik[2]"}. If \code{NULL} the `minus log likelihood' extracted from the object, namely \code{"-logLik(x)"}, is used. See \code{\link{segmented.default}}. } + \item{digits}{optional. If specified it means the desidered number of decimal points of the breakpoint to be used during the iterative algorithm.} } \details{ diff --git a/man/segmented-package.Rd b/man/segmented-package.Rd index 08cf621..2cdaf21 100644 --- a/man/segmented-package.Rd +++ b/man/segmented-package.Rd @@ -3,23 +3,23 @@ %\alias{segmented} \docType{package} \title{ -Segmented relationships in regression models with breakpoints/changepoints estimation +Segmented relationships in regression models with breakpoints / changepoints estimation } \description{ -Estimation of Regression Models with piecewise linear relationships having a +Estimation and Inference of Regression Models with piecewise linear relationships having a fixed number of break-points. } \details{ \tabular{ll}{ Package: \tab segmented\cr Type: \tab Package\cr -Version: \tab 0.5-1.4\cr -Date: \tab 2015-11-04\cr +Version: \tab 0.5-2.2\cr +Date: \tab 2017-09-18\cr License: \tab GPL\cr } Package \code{segmented} is aimed to estimate linear and generalized linear models (and virtually any regression model) having one or more segmented relationships in the linear predictor. Estimates of the slopes and - of the possibly multiple breakpoints are provided. The package includes testing/estimating + breakpoints are provided along with standard errors. The package includes testing/estimating functions and methods to print, summarize and plot the results. \cr The algorithm used by \code{segmented} is \emph{not} grid-search. It is an iterative procedure (Muggeo, 2003) @@ -39,6 +39,8 @@ Vito M.R. Muggeo <[email protected]> } \references{ +Muggeo, V.M.R. (2016) Testing with a nuisance parameter present only under the alternative: a score-based approach with application to segmented modelling. \emph{J of Statistical Computation and Simulation} \bold{86}, 3059--3067. + Davies, R.B. (1987) Hypothesis testing when a nuisance parameter is present only under the alternative. \emph{Biometrika} \bold{74}, 33--43. @@ -59,6 +61,8 @@ Muggeo, V.M.R., Adelfio, G. (2011) Efficient change point detection in genomic s Wood, S. N. (2001) Minimizing model fitting objectives that contain spurious local minima by bootstrap restarting. \emph{Biometrics} \bold{57}, 240--244. +Muggeo, V.M.R. (2010) Comment on `Estimating average annual per cent change in trend analysis' by Clegg et al., Statistics in Medicine; 28, 3670-3682. +\emph{Statistics in Medicine}, \bold{29}, 1958--1960. } \keyword{ regression } diff --git a/man/segmented.Rd b/man/segmented.Rd index 5f89a43..0d399c4 100644 --- a/man/segmented.Rd +++ b/man/segmented.Rd @@ -37,8 +37,9 @@ segmented(obj, seg.Z, psi, control = seg.control(), fit may be supplied (see 'Details').} \item{seg.Z}{ a formula with no response variable, such as \code{seg.Z=~x1+x2}, indicating the (continuous) explanatory variables having segmented relationships with the response. - Currently, formulas involving functions, such as \code{seg.Z=~log(x1)} or \code{seg.Z=~sqrt(x1)}, - or selection operators, such as \code{seg.Z=~d[,"x1"]} or \code{seg.Z=~d$x1}, are \emph{not} allowed. + Currently, formulas involving functions, such as \code{seg.Z=~log(x1)} or even \code{seg.Z=~sqrt(x1)}, + or selection operators, such as \code{seg.Z=~d[,"x1"]} or \code{seg.Z=~d$x1}, are \emph{not} allowed. + It can be missing when \code{obj} ("lm" or "glm" fit) includes only one covariate which is taken as segmented variable. } \item{psi}{ named list of vectors. The names have to match the variables in the \code{seg.Z} argument. Each vector includes starting values for the break-point(s) for the corresponding diff --git a/man/slope.Rd b/man/slope.Rd index 48b0910..baabdd9 100644 --- a/man/slope.Rd +++ b/man/slope.Rd @@ -33,11 +33,14 @@ slope(ogg, parm, conf.level = 0.95, rev.sgn=FALSE, \code{slope} returns a list of matrices. Each matrix represents a segmented relationship and its number of rows equal to the number of segments, while five columns summarize the results. } -\references{Muggeo, V.M.R. (2003) Estimating regression models with unknown break-points. +\references{ + Muggeo, V.M.R. (2003) Estimating regression models with unknown break-points. \emph{Statistics in Medicine} \bold{22}, 3055--3071. + + Muggeo, V.M.R. (2010) Comment on `Estimating average annual per cent change in trend analysis' by Clegg et al., + Statistics in Medicine; 28, 3670-3682. \emph{Statistics in Medicine}, \bold{29}, 1958--1960. } -\author{Vito M. R. Muggeo, \email{[email protected]} -} +\author{Vito M. R. Muggeo, \email{[email protected]} } \note{The returned summary is based on limiting Gaussian distribution for the model parameters involved in the computations. Sometimes, even with large sample sizes such approximations are questionable (e.g., with small difference-in-slope parameters) and the results returned by \code{slope} @@ -45,7 +48,7 @@ slope(ogg, parm, conf.level = 0.95, rev.sgn=FALSE, approximations. Anyway, the t values may be not assumed for testing purposes and they should be used just as guidelines to assess the estimate uncertainty. } -\seealso{See also \code{\link{davies.test}} to test for a nonzero difference-in-slope parameter. +\seealso{See also \code{\link{davies.test}} and \code{\link{pscore.test}} to test for a nonzero difference-in-slope parameter. } \examples{ set.seed(16) -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-segmented.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
