[R] cycling k times a realization of a random walk.....problems..
hello R-masters. i have an R-issue here that i don't know if you'd wish to help me about it: briefly i'd like to generate many (say hundred) realizations of a random walk, execute a few operations on each of them (mean time of return), and graph each realization on the same plot. IN OTHER WORDS I'D LIKE TO IMPOSE A LOOPING CYCLE TO THE COMMAND NOT THE ARGUMENT OF THE COMMAND. for some of these questions i have already a partial answer: my main problem here is automatizing the process for 100 times. my random walk generating function is: rwalk <- function(n,p, x0=0) { x <- rbinom(n,1,p) x[which(x==0)] <- -1 y<-c(x0,x) y <- cumsum(y) } THIS IS THE CODE THAT I'D LIKE TO REITERATE a<- rwalk(n,p,x0=0) #whish to generate 100 of those #and for each calculate the following o<-which(a==0) N<-length(o) # number of times the walk returns to 0 Nn<-(o[-1]-N) # number of steps to get to 0 V<-mean(Nn) # mean time of return to 0 par(mfrow=c(1,1)) plot(0:n, a, type= "l", main="...", xlab="tempo", ylab="stato", ylim=c(...), lwd=3, col="red", lty=1) par(new=T) plot(0:n, b, type= "l", main="", xlab="", ylab="", ylim=c(), lwd=3, col="green", lty=2) #technically i've tryed to fuse all those parts in a function reiterating k=100 times with a for cycle, but with no succes, and i dare there must be something really disturbing in the code below function(k){ for (i in 1:k){ a[i]<-rwalk(1e+04,0.5,0) o[i]<-which(a[i]==0) N[i]<-length(o[i]) Nn[i]<-(o[i][-1]-N[i]) V[i]<-mean(Nn[i]) w<-c(V[1]:V[k]) M<-mean(w) #mean over all the realizations par(mfrow=c(1,1)) plot(0:n, a[1], type= "l", main="...", xlab="tempo", ylab="stato", ylim=c(...), lwd=3, col="red", lty=1) par(new=T) plot(0:n, a[i], type= "l", main="", xlab="", ylab="", ylim=c(...), lwd=3, col=c(1:i), lty=2) } } #NOTE: I'VE ALSO TRIED TO REITERATE THE WALK WITH THE rep() # t<-rep(rwalk(1e+04,0.5,0),100) # then of course i have the problem of splicing this ridicolously #one milion long vector in 100 bits of tenthousand, and still goig through all the calculations and plotting i meant before i wonder what.. THANK YOU FOR YOUR PRECIUOS ATTENTION!!! __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] the sample() function
hello R-Wizards! again i'm invoking your presence!! since all the fuss about the paradoxes on a computer algorithm generating "caos" or un-determinancy, i recently grew quite curious about the mechanism underlying the procedure of NUMBER RADOMIZATION. could anyone of you, masters, attach me the algorithm (or source code? is it right?) behind the *sample()* function, so i can inspect in detail the mechanism of this so gossiped "radomization"? thank you, sincerely federico bonofiglio, student of statistics at milano bicocca university, italia [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] dpois().......bizarre warning messages
Dear Masters, I have a question to submit consider the following script m<-4.95 obs<-rpois(36,m) # i generate 36 realization from a poisson(m) hist(obs,freq=F) curve(dpois(x,m),add=T,col="red") #i wish to overlay on the histogram the theorical poisson density function errors are returned saing the x vector doesn't contain integers really bizarre i can't give explanation (R version 2.11.1, no source codes) would u be so kind to suggest me a solution??? thank u FB student of statistics at milano bicocca [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] lowess() won't handle NAs
Dear Masters, I'm driving crazy with the lowess() function my intent is smoothing confidence intervals for predicted y values in a linear model lm() setting since in the predict() function there exists an option for predicting NA values, I instead encounter problems when I fit a missing values x variable to the predicted terms.,impossible task!!! I've been banging on my head all morning, I can't find a solution also using a loess() I can't get the function to keep the x missing values with the result that it automaticaly drop 'em, so that the smoothing looks like a mess in the plot I attach my scripts, a function that would automatize the process explained above #d = a data.frame with some Na #c1,c2 =column indexes regr.plot<-function(d,c1,c2, method=c("spearman","pearson"),titolo,xtag,ytag,livello){ test<-cor.test(d[,c2],d[,c1],method="spearman") pval<-round(test$p.val,4) corr<-round(test$estimate,4) I<-predict(lm(d[,c1]~d[,c2]),newdata=d,interval="confidence",level=livello) plot(d[,c2],d[,c1],xlab=as.character(xtag),ylab=as.character(ytag)) title(as.character(titolo),font=3) abline(lm(d[,c1]~d[,c2]),lwd=2,col="red") #lines(loess(I[,2]~d[,c2])#,data=d,control =loess.control(surface = "direct") #,na.action=na.omit(d[,c2]) #,lty=2,lwd=2,col="blue") #lines(loess(I[,3]~d[,c2])#,data=d,control =loess.control(surface = "direct") #,na.action=na.omit(d[,c2]) #,lty=2,lwd=2,col="blue") lines(lowess(d[,c2],I[,2]),lty=2,lwd=2,col="blue") lines(lowess(d[,c2],I[,3]),lty=2,lwd=2,col="blue") text(locator(1),as.expression(substitute(r==corr),list(corr=corr)), cex=1.5,font=3) if(pval==0){ text(locator(1),expression("p<0.001"), cex=1.5,font=3) }else{ text(locator(1),as.expression(substitute(p==pval),list(pval=pval)), cex=1.5,font=3) } } #THANK U 4 THE EVENTUAL ADVISE.!! Federico, Student of Statistics at Milano-Bicocca [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] anomalies with the loess() function
Hello Masters, I run the loess() function to obtain local weighted regressions, given lowess() can't handle NAs, but I don't improve significantly my situation.., actually loess() performance leave me much puzzled I attach my easy experiment below #--SCRIPT-- #I explore the functionalities of lowess() & loess() #because I have encountered problems in execute local weighted regressions #with lowess() (in presence of NAs) & with loess() (always!!!) #I generate 2 fictious vectors a<-sample(c(sample(1:1000,100),rep(NA,50))) b<-sample(c(sample(1:1000,100),rep(NA,50))) #lm() has no problems..can handle the missing values plot(a,b) abline(lm(b~a),col="red",lwd=2) #loess return a plain mess like it would go dizzed with ordering or something. #Off course lowess() turns useless in presence of NAs, I don't even try it. lines(loess(b~a)) #I get rid off NAs and compare lowess() & loess() performance, expecting to #obtain the same result as both functions implement local weighted regressions a<-na.omit(a) b<-na.omit(b) #check out the evidence.something's wrong with loess()??? par(mfrow=c(1,2)) plot(a,b) lines(lowess(a,b),col="red")#if NAs are excluded lowess() runs regularly plot(a,b) lines(loess(b~a),col="red")#.but loess() keeps messing all over...!!??? [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] logistic regression with response 0,1
Dear Masters, first I'd like to wish u all a great 2011 and happy holydays by now, second (here it come the boring stuff) I have a question to which I hope u would answer: I run a logistic regression by glm(), on the following data type (y1=1,x1=x1); (y2=0,x2=x2);..(yn=0,xn=xn), where the response (y) is abinary outcome on 0,1 amd x is any explanatory variable (continuous or not) observed with the i-th outcome. This is indeed one of the most frequent case when challenged with binary responses, though I know R manages such responses slightly differently (a vector for the successes counts and one for the failures) and I'm not sure wheather my summary.glm gives me any senseful answer at all for the purpose I have tried to to emphasize the numbers so to obtain significant results y<-rbinom(2000,1,.7)#response for(i in 1:2000){ euro[i]<-if(y[i]==1){rnorm(1,300,20)#explanatory 1 }else{rnorm(1,50,12)} } for(i in 1:2000){ sex[i]<-if(y[i]==1){sample(c("m","f"),1,prob=c(.8,.2))#explanatory 2 }else{sample(c("m","f"),1,prob=c(.2,.8))} } m<-glm(y~euro+factor(sex),family=binomial) summary(m) My worries: - are the estimates correct? - degrees of freedom exponentiate dramatically (one per cell) , so may I risk to never obtain a significant result? I also take the chance to ask wheater u know any implemented method to plot logistic curves directly out of a glm() model I would like to thank u all by the way Federico Bonofiglio [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] graphics: 3D regression plane
Hello Masters, wishing you all a great 2011 I was also going to ask if anyone knows a quick and efficient way to plot a regression plane (z~x*y). I have tried the regr2.plot{HH} function but it is only an educational tool and has poor graphical properties. I also tried to run the following script on a fictitious longitudinal problem, with poor results set.seed(1234) id<-c(rep(1,3),rep(2,4),rep(3,2)) # subjects y<-rchisq(9,df=20) #response k<-rnorm(9,4,2) # x time<-as.Date(c("03/07/1981","15/11/1981","03/04/1983","08/12/1979", "30/12/1979","08/03/1980","12/08/1980","12/08/1973","28/03/1975"), format="%d/%m/%Y") fac<-c("m","m","m","f","f","f","f","m","m")# sex d1<-as.vector(by(time,factor(id),min)) t0<-as.Date(d1,origin=as.Date("1970-01-01"));t0 A<-data.frame(id=c("1","2","3"),t=t0) B<-data.frame(id=id,tempo=time) C<-merge(A,B);C rd<-as.vector(C$tempo-C$t);rd #time centered on sbj specific first occurrence mod<-lm(y~rd*k) newax<- expand.grid( days = giorni<-seq(min(rd),max(rd), length=100), expl= esplic<- seq(min(k), max(k), length=100) ) fit <- predict(mod,data.frame(rd=giorni,k=esplic)) graph <- persp(x=giorni, y=esplic,fit, expand=0.5, ticktype="detailed", theta=-45) #error : z argument not valid I would be grateful if someone would give me some suggestions. Thank u again and happy new year Federico Bonofiglio [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] piecewise regression
Hello everybody Quick question, if you'd like to throw a little tip: does anyone knows a function that runs piecewise regression models with coefficients estimation and inferences ? Thank you [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] piecewise regression
encer.gra...@structuremonitoring.com> wrote: > > > > > Hello: > > > > > > > > > I just found 58 help pages in 32 contributed packages containing > the > > > terms "piecewise" and "regression", as follows: > > > > > > > > > library(sos) > > > pr <- ???'piecewise regression' # for "piecewise" with "regression" > > > summary(pr) # 58 matches in 32 packages > > > pr # view the 58 links sorted by package in a web browser > > > > > > p.r <- ???"{piecewise regression}" # for "piecewise" followed immediate > by > > > "regression" > > > summary(p.r) # 15 matches in 13 packages > > > p.4 # to view in a web browser > > > > > > > > > Beyond this, the "structchange" package supports looking for > changes > > > that would be, e.g., linear in one segment then jump to a different > linear > > > regime. If you want the lines to connect, then you may want splines. > > > Probably the most common spline is cubic, but you can also get linear > > > splines. One of the most general spline packages is "fda" (functional > data > > > analysis), with multiple companion books by James Ramsay and others. > > > However, there are other spline packages. > > > > > > > > > If this does not answer your question, "PLEASE do read the posting > > > guide http://www.R-project.org/posting-guide.html and provide > commented, > > > minimal, self-contained, reproducible code", as it says in the > signature > > > block for each email on this list. > > > > > > > > > Hope this helps. > > > Spencer Graves > > > > > > > > > > > > On 1/14/2011 6:42 AM, Federico Bonofiglio wrote: > > > > > >> Hello everybody > > >> > > >> Quick question, if you'd like to throw a little tip: > > >> does anyone knows a function that runs piecewise regression models > with > > >> coefficients estimation and inferences ? > > >> > > >> Thank you > > >> > > >>[[alternative HTML version deleted]] > > >> > > >> __ > > >> R-help@r-project.org mailing list > > >> https://stat.ethz.ch/mailman/listinfo/r-help > > >> PLEASE do read the posting guide > > >> http://www.R-project.org/posting-guide.html > > >> and provide commented, minimal, self-contained, reproducible code. > > >> > > >> > > > __ > > > R-help@r-project.org mailing list > > > https://stat.ethz.ch/mailman/listinfo/r-help > > > PLEASE do read the posting guide > > > http://www.R-project.org/posting-guide.html > > > and provide commented, minimal, self-contained, reproducible code. > > > > > > > [[alternative HTML version deleted]] > > > > __ > > R-help@r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > > > Confidentiality Statement: > > This email message, including any attachments, is for th...{{dropped:6}} > > > > __ > > R-help@r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > > -- > Open WebMail Project (http://openwebmail.org) > > [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] TRADUCING lmer() syntax into lme()
Dear Rsociety, I'd like to kingly ask to anyone is willing to answer me how to implement a NON NESTED random effects structure in lme() In particular I've tried the following translation from lmer to lme, as suggested from some web example mod1<-lmer(y~x*z+(x*z|factorA1/factorB)+(x*z|factorA2/factorB)) # y,x,z continuous mod2<-lme(y~x*z, random= pdBlocked(list(pdIdent(~1|factorA1/factorB ),pdIdent(~1|factorA2/factorB In detail check how I've tried to state in mod1 that Iwant to evaluate randomness in the interaction x*z (i.e intercept, slope, interaction) grouped by by a general nesting structure that sets factorA1 and factorA2 as same level effects (hence non nested) and factorB as nested in both. I also must express my momentaneous sheer ignorange on the pdMat objects, thing that prabably is not helping me in the process Kindly Regards Federico Bonofiglio -- *Little u can do against ignorance,it will always disarm u: is the 2nd principle of thermodinamics made manifest, ...entropy in expansion.**But setting order is the real quest 4 truth, ..and the mission of a (temporally) wise dude. * [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] random interaction effect in lmer
Hi dears while modeling an interaction random effect in lmer i receive the instantaneous error message > ldlM4<-lmer(ldl~rt*cd4+age+rf+pharmac+factor(hcv)+ + hivdur+(rt:cd4|id),na.action=na.omit,REML=F) *Warning message: In mer_finalize(ans) : false convergence (8) * I think the matter lies in syntax, 'cause i sistematically receive the same message even when changing response... PS: the above model runs quite well in lme and only rarely I recive : * Errore in lme.formula(fixed = ldl ~ rt + cd4 + age + rf + pharmac + factor(hcv) + : nlminb problem, convergence error code = 1 message = iteration limit reached without convergence (9) * ..reason for which I switch into lmer Thank u in advance foe any hints...;) -- *Little u can do against ignorance,it will always disarm u: is the 2nd principle of thermodinamics made manifest, ...entropy in expansion.**But setting order is the real quest 4 truth, ..and the mission of a (temporally) wise dude. * [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] glht{multcomp} : use with lme {nlme}
Hi dears, I do > CHOL<-lme(chol~rt*cd4+sex+age+rf+nadir+pharmac+factor(hcv)+factor(hbs)+ haartd+hivdur+factor(arv), random= ~rt|id, na.action=na.omit) ...runs sweet,..then try a multicomparisons approach for the categorical rf > summary(glht(CHOL, linfct=mcp(rf="Tukey"))) * Error in model.frame.default(object, data, xlev = xlev) : l'oggetto non è una matrice * <-- Object isn't a matrix!!* Errore in factor_contrasts(model) : no model.matrix method for model found! Errore in summary(glht(CHOL, linfct = mcp(rf = "Tukey"))) : error in evaluating the argument 'object' in selecting a method for function 'summary'* data are ok...I have followed the examples on ?glht ..,there is still something I'm missing here? Thanks in advance if you wish to give a tiny answer. -- *Little u can do against ignorance,it will always disarm u: is the 2nd principle of thermodinamics made manifest, ...entropy in expansion.**But setting order is the real quest 4 truth, ..and the mission of a (temporally) wise dude. * [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] intervals {nlme} lower CI greater than upper CI !!!????
Hi folks... check this out.. > GLU<-lme(gluc~rt*cd4+sex+age+rf+nadir+pharmac+factor(hcv)+factor(hbs)+ + haartd+hivdur+factor(arv), + random= ~rt|id, na.action=na.omit) > intervals(GLU)$fixed lower est. upper (Intercept) 67.3467070345 7.362307e+01 7.989944e+01 rt *0.0148050160* 6.249304e-02 1.101811e-01 cd4 -0.0032133709 -6.566501e-04 1.900071e-03 sexM 0.9601467541 3.565606e+00 6.171066e+00 age *0.2436472425* 3.420025e-01 4.403578e-01 rfOTHER -5.0933264232 -1.198439e+00 2.696449e+00 rfETEROSESSUALE -2.1085096013 1.974725e+00 6.057960e+00 rfOMOSESSUALE -5.8156940466 -1.891870e+00 2.031953e+00 nadir -0.0063302516 8.604450e-04 8.051142e-03 pharmacNRTI -4.2586487617 -1.252531e+00 1.753586e+00 pharmacPI -3.3132784651 -8.829819e-01 1.547315e+00 factor(hcv)1-0.4338530808 2.374438e+00 5.182729e+00 factor(hbs)1-0.7670414150 3.274923e+00 7.316888e+00 haartd -0.0837790540 -4.655330e-02 -9.327549e-03 hivdur -0.1351044027 8.417971e-02 3.034638e-01 factor(arv)1-6.7807591712 -2.075808e+00 2.629143e+00 rt:cd4 -0.0001150701 -4.667572e-05 2.171870e-05 attr(,"label") [1] "Fixed effects:" ...check for instance on row 2 and 5 (bold numbers)isn't something really weird there ? lower interval greater than upper !! Federico -- *Little u can do against ignorance,it will always disarm u: is the 2nd principle of thermodinamics made manifest, ...entropy in expansion.**But setting order is the real quest 4 truth, ..and the mission of a (temporally) wise dude. * [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] PREDICTIONS from a PIECEWISE LINEAR (mixed) MODEL: THEY AIN'T LINEAR BETWEEN BREAK POINTS!!
I Dears, if that wouldn't take u too much effort I'd like to ask for a swift opinion: I have a alinear (mixed effect) model that I wish to run as a piecewise one. When I predict the values Iget quite some odds and disturbing results: The predicted stright line after the break point is not straight at all, instead behaves like if it was a hig order polynomial or something similar I attach the codes below, hoping someone can point me the mistake I sincerely express many thanks in advance .. Federico Bonofiglio frame3<-data.frame(id,chol,cd4,rt,sex,age, nadir,pharmac2,hbs,hcv,resp2, hivbl,switch) # I run a model specifing the function (t-t*)+ = 0 if t<=t* AND (t-t*)+ = t if t>t* with the following syntax, emulating MJ Crawley (The R book 2007, Wiley), # rt+rt*(rt>15) I run the model, is a mixed effect one... PCWISE<-lme(chol~rt+rt*(rt>15) +sex+age+cd4+cd4:rt+ cd4:rt*(rt>15)+nadir+hivbl+ pharmac2:rt+resp2:rt+ resp2:rt*(rt>15)+hbs+hcv+ switch,data=frame3, random= ~rt|id, na.action=na.omit) # I prepare a fictious data.frame for the predictions , in fact lme predict values inside the clusters , id in my case, they are patients.. I need instead a population result... nx<-seq( min(cd4,na.rm=T), max(cd4,na.rm=T), length.out=200) nt<-seq( min(rt,na.rm=T), max(rt,na.rm=T), length.out=200) newframe<-data.frame(#generate a fictious patient's profile id="48", cd4=nx, rt=nt, sex=rep("M",200), age=rep(mean(age),200), nadir=rep(mean(age),200), pharmac2=rep("PI.boost",200), resp2=rep("mild",200), hbs=rep(1,200), hcv=rep(1,200), hivbl=rep(mean(hivbl,na.rm=T),200), switch=rep(1,200)) yy<-predict(PCWISE,level=0,newdata=newframe) #predict the response # I plot the predicted values against the observed.. xyplot(chol~rt, groups=id, panel=function(x, y){panel.xyplot(x,y) panel.abline(v=15,col=2,lty=2) panel.lines(nt,yy,col=2,lwd=2) # plot the predicted values panel.axis("top",at=15,labels="15") }) # I get a quite puzzling result : predictions are not linear they act more like they were from a polynomial of high order!!! for instance after time 15, (rt>15), I get a line that is not straight, but flexes smoothly downward #Also if I try run the model separetely for time <= 15 and >15... prova2<-lme(chol~rt+sex+age+cd4+cd4:rt+ nadir+hivbl+ pharmac2:rt+resp2:rt+hbs+hcv+ switch,data=frame3, random= ~rt|id, subset=which(rt<=15), na.action=na.omit,correlation=corAR1()) prova3<-lme(chol~rt+sex+age+cd4+cd4:rt+ nadir+hivbl+ pharmac2:rt+resp2:rt+hbs+hcv+ switch,data=frame3, random= ~rt|id, subset=which(rt>15), na.action=na.omit,correlation=corAR1()) # ..I get quite completely different estimeates for time and all the time interactions (sometimes an order of magnitude!) -- *Little u can do against ignorance,it will always disarm u: is the 2nd principle of thermodinamics made manifest, ...entropy in expansion.**But setting order is the real quest 4 truth, ..and the mission of a (temporally) wise dude. * [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] strange PREDICTIONS from a PIECEWISE LINEAR (mixed) MODEL
Hi Dears, When I introduce an interaciton in a piecewise model I obtain some quite unusual results. If that would't take u such a problem I'd really appreciate an advise from you. I've reproduced an example below... Many thanks x<-rnorm(1000) y<-exp(-x)+rnorm(1000) plot(x,y) abline(v=-1,col=2,lty=2) mod<-lm(y~x+x*(x>-1)) summary(mod) yy<-predict(mod) lines(x[order(x)],yy[order(x)],col=2,lwd=2) #--lme #grouping factor, unbalanced g<-as.character(c(1:200)) id<-sample(g,size=1000,replace=T, prob=sample(0:1,200,rep=T)) table(id) #unbalanced mod2<-lme(y~x+x*(x>-1),random=~x|id, data=data.frame(x,y,id)) summary(mod2) newframe<-data.frame( #fictious id id="fictious", x) newframe[1:5,] #predictions yy2<-predict(mod2,level=0, newdata=newframe) lines(x[order(x)],yy2[order(x)],col="blue",lwd=2) # add variable in the model z<-rgamma(1000,4,6) mod3<-lme(y~x+x*(x>-1)+z ,random=~x|id, data=data.frame(x,y,z,id)) summary(mod3) #new id newframe2<-data.frame( #fictious id id="fictious", x, z) #predict yy3<-predict(mod3,level=0, newdata=newframe2) lines(x[order(x)],yy3[order(x)],col="green",lwd=2) # ADD INTERACTION z:x mod4<-lme(y~x+x*(x>-1)+ z+ z:x+ z:x*(x>-1) ,random=~x|id, data=data.frame(x,y,z,id)) #predict yy4<-predict(mod4,level=0, newdata=newframe2) lines(x[order(x)],yy4[order(x)],col="violet",lwd=2) #something bizarre #starts to happen #in the predicted values # they begin to jiggle around the straight line -- *Little u can do against ignorance,it will always disarm u: is the 2nd principle of thermodinamics made manifest, ...entropy in expansion.**But setting order is the real quest 4 truth, ..and the mission of a (temporally) wise dude. * [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.