[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.
[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 tt*
with the following syntax, emulating MJ Crawley (The R book 2007, Wiley),
# rt+rt*(rt15)
I run the model, is a mixed effect one...
PCWISE-lme(chol~rt+rt*(rt15)
+sex+age+cd4+cd4:rt+
cd4:rt*(rt15)+nadir+hivbl+
pharmac2:rt+resp2:rt+
resp2:rt*(rt15)+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, (rt15), 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(rt15),
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] 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] 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] 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.
Re: [R] piecewise regression
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] 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.
[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] 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] 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] 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(p0.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] 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] 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] 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.
