Re: [R] Automatic Knot selection in Piecewise linear splines

dear all, I apologize for my delay in replying you. Here my contribution, maybe just for completeness: Similar to "earth", "segmented" also fits piecewise linear relationships with the number of breakpoints being selected by the AIC or BIC (recommended). #code (example and code from Martin

[R] [R-pkgs] segmented 2.1-0 is released

dear R users, I am pleased to announce that segmented 2.1-0 is now available on CRAN. segmented focuses on estimation of breakpoints/changepoints of segmented, i.e. piecewise linear, relationships in (generalized) linear models. Starting with version 2.0-0, it is also possible to model stepmen

Re: [R] losing variable attributes when subsetting model.frame

dear Duncan, Thank you very much for your quick reply. Very useful and also very elegant, I would add. Thanks again, kind regards, Vito Il 08/05/2024 14:06, Duncan Murdoch ha scritto: I would guess that subsetting uses [] indexing on each column, with a logical argument.  If that is true, t

[R] losing variable attributes when subsetting model.frame

dear all, I have a simple function f() which, when included in model.frame() via the formula, returns the variable itself with some attributes. However when I specify the subset argument, the attributes get lost, apparently. I would like to extract the attributes also when specifying the subse

[R] OT: IWSM 2013

. best wishes, Vito Muggeo, on behalf of the IWSM2013 Scientific Committee -- == 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:

Re: [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=right

Re: [R] fit a "threshold" function with nls

Véronique, in addition to Bert's comments, I would like to bring to your attention that there are several packages that perform threshold/breakpoint/changepoint estimation in R, including cumSeg, segmented, strucchange, and bcp for a Bayesian approach Moreover some packages, such as cghFLasso

Re: [R] Comparing crossing survival curves

hi isabel, You have to decide if focus is on the survival curves or hazards.. Crossing hazards do not imply crossing survival curves If you are dealing with crossing hazards, and you are interested in testing for an effect of a covariate (presumably with a crossing hazard effect), then a stan

Re: [R] Piecewise Lasso Regression

tructure and I believe that fitting a piecewise regression would be of great benefit. Thanks, Lucas On Jun 6, 2012, at 4:54 AM, Vito Muggeo (UniPa) wrote: dear lucas, yes you are right, segmented does not handle 'lars' objects. Out of curisity, are you interested in selecting the nu

Re: [R] Piecewise Lasso Regression

dear lucas, yes you are right, segmented does not handle 'lars' objects. Out of curisity, are you interested in selecting the number of breakpoints or in selecting additional covariates with linear parameters? vito Il 06/06/2012 0.01, Lucas Santana dos Santos ha scritto: Hi All, I am tryi

Re: [R] Finding multiple breakpoints - 'segmented' ?

dear Peter, Currently segmented handles multiple breakpoints for several variables (the limit discussed in the msg 2006 has been fixed..). However you are looking for a somewhat complicated model where the breakpoint of the relationship 's.size' and 'R.AUC' depends on another covariate 'bedek

[R] getting the name of the working .Rdata file

dear all, I do not if it is a nonsense question.. Is it possible in the R session to get the name of the current .Rdata file that I ran? I mean: suppose I double click the file myfile.Rdata. ls() returns the names of the objects in the current workspace (that is saved in myfile.Rdata). In th

Re: [R] Breakpoint in logistic GLM with 'segmented' package - error: replacement length zero

dear Peter, Your code appears correct, so it is difficult to reply without the data.. If you are interested in further details, please contact me off-list vito Il 25/05/2012 15.34, Peter Hoitinga ha scritto: Hello all, I've been having trouble with assessing a breakpoint in a logistic GLM wi

Re: [R] Error with psi value for 'segmented' package for R

dear Szymon, what do you mean "it does not work for others.. that fit within similar range"? Each dataset has its own features and breakpoint estimation is not as simple as estimation of linear models even if your data "fit within similar range". I will contact you out of the list for detai

Re: [R] Error with the 'segmented' package for R

dear Szymon, it is a bug (in the new version), thanks. It depends on the flat underlying relationship you are trying to estimate with a small sample.. I will correct it as soon as possible. Meanwhile you can use o1<-glm(gpp ~ temp) os1<-segmented(o1, seg.Z=~temp, psi=15, control=seg.control(n.

Re: [R] variable dispersion in glm models

there are at least two alternatives 1) package dglm for Double generalized linear models or probably better 2)package gamlss for Generalized Additive Models for Location Scale and Shape. Here you can use alternative, sometimes more appropriate, families and also you can include additive (nonp

Re: [R] linear-by-linear association model in R?

dear Christofer, Try the following d<-expand.grid(a=1:3,b=1:4) d$freq<-rpois(12,5) o<-glm(freq~factor(a)+factor(b)+I(a*b), family=poisson, data=d) vito Il 02/04/2012 9.34, Christofer Bogaso ha scritto: Dear all, can somebody give me some pointer how I can fit a "linear-by-linear association

[R] glmnet() vs. lars()

dear all, It appears that glmnet(), when "selecting" the covariates entering the model, skips from K covariates, say, to K+2 or K+3. Thus 2 or 3 variables are "added" at the same time and it is not possible to obtain a ranking of the covariates according to their importance in the model. On t

Re: [R] breakpoints and nonlinear regression

dear Julian, Il 18/01/2012 14.36, crimsonengineer87 ha scritto: Thanks for the comments. Yes, I also had segmented and then I went away from that. I can't remember. I've tried using it but I get some sort of strange error. Here's some code ... it is difficult for me to help you without knowin

[R] looking for an ecological dataset

dear all, apologizes for this off-topic question. I am looking for a "ecological" dataset (n>100, say) including measurements of one or more growth variable and age. Could anyone to suggest the R package/URL where I can find it? many thanks, vito -- Vit

Re: [R] Problem with segmented

dear Phil, I am not able to read the error message.. did you forget it? However: does x exist in the workspace? The following lines work: myreg2 = lm(y ~ x, data=xy) mysegmented = segmented(myreg2, seg.Z=~x, psi=c(245000)) myreg2 = lm(xy$y ~ xy$x) x<-xy$x mysegmented = segmented(myreg2, seg.Z=

Re: [R] mgcv gam predict problem

dear Philip, I am not able to solve your problem, however the error message you get does not depends on mgcv::gam, therefore gam(,..outer.ok=TRUE) or predict.gam(,outer.ok=TRUE) do not make sense. The error message comes from the function splines::splineDesign which is called when the option

Re: [R] plot not generic

dear Nick, getAnywhere("plot.glmnet") Note the message you get when you type methods(plot) ... >>> Non-visible functions are asterisked Il 28/01/2011 14.26, Nick Sabbe ha scritto: Hello list. I was trying to see some of the code for plot.glmnet in package glmnet (this function name i

Re: [R] Choosing statistical test - Fisher's Exact Test?

It appears that you have a 2x2 table coming from paired binary data.. If this is the case the McNemar test is appropriate. See ?mcnemar.test or even better the package exact2x2, function mcnemar.exact() for an "exact" approach, vito Il 18/01/2011 14.40, debz ha scritto: Hi I was wonderi

Re: [R] Convert a matrix's columns to list

hi feng, a possible solution is b1<-apply(a,2,list) and possibly lapply(b1,unlist) if you want exactly the output equal to "list(a[, 1], a[, 2])" best, vito Il 18/01/2011 13.53, Feng Li ha scritto: Dear R, Is there an efficient way to make a list that each element is from the corresponding

Re: [R] confidence interval for logistic joinpoint regression from package ljr

dear M., I do not know how to get the SE for the joinpoint (or breakpoint) from your ljr fit. However you can find useful the segmented package which works for any GLM (including the logistic one) and it returns (approximate) StErr (and Conf Int) also for the joinpoint (breakpoint in the segme

Re: [R] One silly question about "tapply output"

dear Vincy, Firstly, a suggestion: to increase the probability of getting help, you should provide reproducible code (people can do "copy-and-paste" of your code and to modify the code to obtain the response.. ) However a possible solution (not tested, of course..) could be simply a<-tapply(r

[R] building lme call via call()

dear all, I would like to get the lme call without fitting the relevant model. library(nlme) data(Orthodont) fm1 <- lme(distance ~ age, random=list(Subject=~age),data = Orthodont) To get fm1$call without fitting the model I use call(): my.cc<-call("lme.formula", fixed= distance ~ age, random =

Re: [R] Interactions in GAMMs

Hi Karen, I think you should decide what you mean for "interaction". s(x:y) is meaningless If you want to fit a surface you should use s(x,y). If you want to fit a varying coefficient model (interaction between a linear and a smooth term) you should use the argument by in s(). The help file

Re: [R] Peak Over Threshold values

dear Tonja, By plotting your data plot(df) it seems to me that you are looking for a piecewise linear relationships. If this is the case, have a look to the package segmented. You have to specify or not the number and the starting values for the breakpoints library(segmented) olm<-lm(waleve

[R] only actual variable names in all.names()

dear all, When I use all.vars(), I am interest in extracting only the variable names.. Here a simple example all.vars(as.formula(y~poly(x,k)+z)) returns [1] "y" "x" "k" "z" and I would like to obtain "y" "x" "z" Where is the trick? many thanks vito -- Vi

[R] partial matching with grep()

dear all, This is a probably a silly question. If I type > grep("x",c("a.x" ,"b.x","a.xx"),value=TRUE) [1] "a.x" "b.x" "a.xx" Instead, I would like to obtain only "a.x" "b.x" How is it possible to get this result with grep()? many thanks for your attention, best, vito --

Re: [R] Estimation in a changepoint regression with R

There are at least two R packages dealing with changepoint estimation, segmented and strucchange. Two possible relevant papers are available: 1)Journal of Statistical Software for strucchange (2002, Vol.7, Issue2) 2)Rnews for segmented (2008, 8/1: 20-25) Hope this helps you vito FMH ha scritt

Re: [R] Fitting a linear model with a break point

dear Dan, As far as I know, the strucchange package can be helpful for you.. On the other hand, if your regression function is continuous at the unknown break points to be estimated, you could try the segmented package. Hope this helps you, vito Daniel Brewer ha scritto: Hello, I would lik

Re: [R] Zinb for Non-interger data

I think that the (impressive) gamlss package (see http://www.gamlss.com) may be helpful. If I remember correctly, in gamlss you can fit model with zero-inflated continuous distributions hope this helps you, vito Alain Zuur ha scritto: JPS2009 wrote: Sorry bit of a Newbie question, and I

Re: [R] Linear Regression Problem

dear Alex, I think your problem with a large number of predictors and a relatively small number of subjects may be faced via some regularization approach (ridge or lasso regression..) hope this helps you, vito Alex Roy ha scritto: Dear All, I have a matrix say, X ( 100 X 40

Re: [R] zero cells in one variable in logistic regression

dear anna, if you are not interested in point estimate and SE of the parameter of the aforementioned categorical variable, I believe the conventional glm(..,family=binomial) is correct. In particular, the returned deviance is reliable and also it is the relevant likelihood ratio test.. hope t

Re: [R] gradually switching regression

Hi, Although Bacon&Watts (1971) assume a transition, the aim is to estimate the breakpoint where the linear relationship changes. The start- and end-point of the transition phase are not parameters to estimate ; it is a trick to estimate the model. As a possible alternative, you could have a

Re: [R] Linear Models: Explanatory variables with uncertainties

Probably you are looking for EIV (errors-in-variables) or ME (measurement errors) models. "simex" is a possible package which needs to know the error variance.. Also RSiteSearch() may be helpful.. hope this helps, vito kpal ha scritto: One of the assumptions, on which the (General) Linear Mo

Re: [R] Discover significant change in sorted vector

dear Hans, As pointed out by Petr Pikal, segmented allows to estimate breakpoints of (continuous) piecewise relationships. That is, the mean lines are assumed to be connected and segmented tries to join them at the estimated breakpoints. The estimated breakpoints may be any value in the range

Re: [R] Piecewise

Hi, In addition to useful Ben's suggestion, you have a, possibly simpler, alternative. If you are willing to assume to know the power of you piecewise polynomial (beta parameter according to code of Ben) you can use the package segmented. Using the data generated by Ben in his previous email

Re: [R] piecewise linear regression

Sorry for my delay.. If you do not know the breakpoint, I would suggest to estimate it.. Have a look to the segmented package. The relevant code here is attach(d) m0<-lm(percent~ year , weights=1/se) library(segmented) mseg<-segmented(m0,seg.Z=~year,psi=1995) points(year, fitted(mseg)) Hope th

Re: [R] harmonic function fiting? how to do

dear Yogesh It appears that your model based on parametric terms is too inflexible.. A better alternative to parametric harmonic terms is a spline-based approach, may be cyclic splines.. Have a look to the mgcv package.. vito Yogesh Tiwari ha scritto: Dear R Users, I have a CO2 time series

Re: [R] R equivalent of SAS Cochran-Mantel-Haenszel tests?

Dear Michael, It sounds as a linear-by-linear loglinear model (and its variants) which uses scores for one or more variables in the table.. (see Agresti, 1990, Categorical Data Analysis. I do remember the pages and I have not the book here..) If this is the case, you can use standard call to

Re: [R] holidays effect

Gabor Grothendieck ha scritto: One possibility if you don't have to have days is to reduce it to a weekly or monthly series. Alternatively you can put a dummy variable (1=holiday and zero otherwise) in the regression model for your response. For instance, you could use the xreg argument of t

[R] OFF topic testing for positive coeffs

Dear all, This is off-topic, however I hope someone can give me useful suggestion.. Given the regression model y = b0 + b1*x + e I am interested in testing for positive coeffs, namely H0: b0>0 AND b1>0 H1: b0,b1 unconstrained It is simple to estimate the model under H0 and H1 (there are several

Re: [R] OT: (quasi-?) separation in a logistic GLM

dear Gavin, I do not know whether such comment may be still useful.. Why are you unsure about quasi-separation? I think that it is quite evident in the plot plot(analogs ~ Dij, data = dat) Also it may be useful to see the plot of the monotone (profile) deviance (or the log-lik) for the coef of

Re: [R] Fitting weibull, exponential and lognormal distributions to left-truncated data.

Hi Gough, A possible solution is to use the survreg() in the survival package without specifying the covariates, i.e. library(survival) survreg(Surv(..)~1, dist="weibull") where Surv(..) accepts information about "times", censoring/truncation variables and dist allows to specify alternative d

Re: [R] NA's in segmented

Dear Tyler, Yes the problem is with NA.. There are two solutions: 1) You can use lm() + segmented (you fit a gaussian model, so why do you use glm()?) 2)If you want to use glm()+ segmented(), use na.omit() to pass your dataframe to the data argument of glm, glm(.., data=na.omit()) Also, if

[R] difference between MASS::polr() and Design::lrm()

Dear all, It appears that MASS::polr() and Design::lrm() return the same point estimates but different st.errs when fitting proportional odds models, grade<-c(4,4,2,4,3,2,3,1,3,3,2,2,3,3,2,4,2,4,5,2,1,4,1,2,5,3,4,2,2,1) score<-c(525,533,545,582,581,576,572,609,559,543,576,525,574,582,574,471,59

[R] error in using by + median

dear all, Could anyone explain me the behaviour of median() within by()? (I am running R.2.7.0) thanks, vito > H<-cbind(rep(0:1,l=20),matrix(rnorm(20*2),20,2)) > by(H[,-1],H[,1],mean) INDICES: 0 V1 V2 -0.2101069 0.2954377 -

Re: [R] Comparing switchpoints from segmented

Dear Rob, Rob Knell ha scritto: > Hello everyone > Not strictly an R question but close... hopefully someone will be able > to help. I wish to compare the switchpoints in two switchpoint > regressions. The switchpoints were estimated using the segmented > library It's a package. running in

Re: [R] R square for Monotone regression

Dear Thierry, I do not know fdrtool::monoreg(), however I think that the simple R^2 is meaningless in thise case. An unconstrained (non-monotonic) fit will have always a higher R^2 and more degrees-of-freedom. A better approach would use any goodness-of-fit criterion penalized for the degrees o

[R] strange behaviour of is.factor()

Dear all, It appears that the function is.factor() returns different results when used inside the apply() function: that is, is.factor() fails to recognize a factor.. Where is the trick? many thanks, vito > df1<-data.frame(y=1:10,x=rnorm(10),g=factor(c(rep("A",6),rep("B",4 > is.factor(df1

[R] differences in using source() or console

Dear all, Is there *any* reason explaining what I describe below? I have the following line myfun(x) If I type them directly in R (or copy/past), it works.. However if I type in R 2.6.1 > source("code.R") ##code.R includes the above line Error in inherits(x, "data.frame") : object "d" not foun

Re: [R] Segmented regression

Dear Brendan, I am not sure to understand your code.. It seems to me that your are interested in fitting a one-breakpoint segmented relationship in each level of your grouping variable If this is the case, the correct code is below. In order to fit a segmented relationship in each group you have

[R] using names with functions..

Dear all, I have the following (rather) strange problem.. For some reasons, I finally work with a variable whose name includes an R function, "a.log(z)", say. And that is a problem when I call it in a formula, for instance: > myname<-"a.log(z)" > dd<-data.frame("a.log(z)"=1:10,y=rnorm(10)) >

Re: [R] how to extract the original data of a glm object

See and *read* the help file ?glm the object returned by glm() includes the `data' component hence: aa<-glm(..) aa$data or also eval(aa$call$data) leffgh ha scritto: > my function is > glm(a~log(b)+c+d+e,family=binomial,data=f)->aa > > > I want to extract the original data set of a

Re: [R] R segmented package

Dear Maura, For package-specific questions, it is probably advisable to contact privately the author.. BTW: > However I cannot understand if a break-point per independent variable > (regressor) is handled when there are many independent variables in the > model, or if "segmented" can handle ma

Re: [R] Splines

As reported in ?spline, you should use splinefun() instead. ff<-splinefun(x,y) ff(x0) where x0=5.25 is your case. best, vito stat stat ha scritto: > I want to fit a cubic spline of x on y. where : > > x > [1] 467 468 460 460 450 432 419 420 423 423 > y > [1] 1 2 3 4 5 6 7 8 9