Re: [R] Win upgrade pb (virus)
Thanks to all that responded privately!
If you have AVG and needs the colorspace (or ggplot2) package do the
following from a suitably privileged account:
1) Open AVG gui and select Components-Anti virus from the menu
2) Click Advanced settings... and then go to Temporary disable AVG
protection Here you need half a million clicks to disable protection
for 10 min
3) Install colorspace package from R (I installed ggplot2 that installs
colorspace)
4) Go to the user library and note the path to the colorspace.dll . In
my case it is:
C:\Users\mvalle\R\win-library\2.14\colorspace\libs\i386\colorspace.dll
(SHIFT-click and then Copy as path is handy)
5) In the open AVG window click Manage Exceptions then Add File and
add the path above
6) On the main AVG gui there is a big button that reenables protection.
Voilà, no more problems.
Thanks again!
mario
On 11-Nov-11 11:34, Mario Valle wrote:
I just upgraded my Win7 32bits installation to 2.14.0 after
deinstalling 2.12.x
First thing I moved the win-library from 2.12 to 2.14 and executed a
update.packages(ask='graphics',checkBuilt=TRUE)
(Swiss mirror).
This aborts with the console message:
Error in if (any(diff)) { : missing value where TRUE/FALSE needed
And the antivirus (AVG) pops up complaining that colorspace.dll
contains the virus Win32/Heur
I cannot ignore the thread because update.packages has already crashed
when I can reach the antivirus dialogbox to push the ignore button.
To continue I had to remove the colorspace and ggplot2 packages.
Is it a known problem? What else can I do (except removing the above
packages or changing antivirus)?
Thanks!
mario
--
Ing. Mario Valle
Data Analysis and Visualization Group| http://www.cscs.ch/~mvalle
Swiss National Supercomputing Centre (CSCS) | Tel: +41 (91) 610.82.60
v. Cantonale Galleria 2, 6928 Manno, Switzerland | Fax: +41 (91) 610.82.82
__
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] Estimating IRT models by using nlme() function
Note that there are several IRT dedicated packages in R, such as the ltm
and eRm packages. For more info you may check the Psychometrics Task
View at: http://cran.r-project.org/web/views/Psychometrics.html
I hope it helps.
Best,
Dimitris
On 11/11/2011 9:38 PM, Cengiz Zopluoğlu wrote:
Hi,
I have a question about estimating IRT models by using nlme, not just rasch
model, but also other models.
Behavior Research Methods
http://www.springerlink.com/content/1554-351x/ Volume
37, Number 2http://www.springerlink.com/content/1554-351x/37/2/, 202-218,
DOI: 10.3758/BF03192688
Using SAS PROC NLMIXED to fit item response theory models (2005). Ching-Fan
Sheuhttps://mail.google.com/mail/u/0/html/compose/static_files/goog_851944013,
Cheng-Te
Chenhttps://mail.google.com/mail/u/0/html/compose/static_files/goog_851944013,
Ya-Hui Su and Wen-Chung Wang
An example of this is provided in the paper above. They use SAS PROC
NLMIXED, and estimate all possible binary and Polytomous IRT models. I want
to replicate what they did using R and nlme package. I am able to fit rasch
model, but have some error messages for 2PL IRT model. So, I could not go
beyond that.
This is really important for me, because any nonlinear model can be fitted
by using this approach. I would like to see how well nlme() package
recovers true parameters, not just regular IRT models but also some other
less used IRT models.
Here is the example. Let's say I have the following dataset. 1000 people
responds five items.
ID Item Response d1 d2 d3 d4 d5
1.1 110 1 0 0 0 0
1.2 120 0 1 0 0 0
1.3 130 0 0 1 0 0
1.4 141 0 0 0 1 0
1.5 150 0 0 0 0 1
2.1 211 1 0 0 0 0
2.2 220 0 1 0 0 0
2.3 230 0 0 1 0 0
2.4 241 0 0 0 1 0
2.5 250 0 0 0 0 1
3.1 310 1 0 0 0 0
3.2 320 0 1 0 0 0
3.3 331 0 0 1 0 0
3.4 341 0 0 0 1 0
3.5 350 0 0 0 0 1
..
..
1000.1 100011 1 0 0 0 0
1000.2 100020 0 1 0 0 0
1000.3 100031 0 0 1 0 0
1000.4 100041 0 0 0 1 0
1000.5 100050 0 0 0 0 1
The items are nested within students. Response is the actual dependent
variable. d1, d2, d3, d4, and d5 are the corresponding dummy codes for each
item. We treat persons as random effects, and item parameters as fixed
effects. To fit Rasch model, I run the following code:
d- read.table(data.txt, header=TRUE)
d1- d$d1
d2- d$d2
d3- d$d3
d4- d$d4
d5- d$d5
###
onePL- function(b1,b2,b3,b4,b5,theta) { #nonlinear model to fit
b= b1*d1+b2*d2+b3*d3+b4*d4+b5*d5
exp((theta-b))/(1+exp((theta-b)))
}
#
nlme(model=Response ~ onePL(b1,b2,b3,b4,b5,theta),
data = d,
fixed = b1+b2+b3+b4+b5 ~ 1,
random = theta ~ 1 | ID,
start=c(b1=0,b2=0,b3=0,b4=0,b5=0),
)
*OUTPUT *
Nonlinear mixed-effects model fit by maximum likelihood
Model: Response ~ onePL(b1, b2, b3, b4, b5, theta)
Data: d
Log-likelihood: -2597.344
Fixed: b1 + b2 + b3 + b4 + b5 ~ 1
b1 b2 b3 b4 b5
-1.1240371 1.5931634 0.2574785 -2.0993236 0.8881437
Random effects:
Formula: theta ~ 1 | ID
theta Residual
StdDev: 0.9765999 0.3780802
Number of Observations: 5000
Number of Groups: 1000
This make sense to me. Actually, the data was simulated data, and the b
parameters were close to true values used to generate data. Also the
standard deviation of random effects (theta or latent ability level in this
case) was close to 1 that was used to generate IRT person ability
parameters.
Then, I run the following code to estimate 2 PL IRT model. It was all same,
just add an additional a parameter for each item.
twoPL- function(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,theta) {
a= a1*d1+a2*d2+a3*d3+a4*d4+a5*d5
b= b1*d1+b2*d2+b3*d3+b4*d4+b5*d5
exp(a*(theta-b))/(1+exp(a*(theta-b)))
}
nlme(model=Response ~ twoPL(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,theta),
data = d,
fixed = a1+a2+a3+a4+a5+b1+b2+b3+b4+b5 ~ 1,
random = theta ~ 1 | ID,
start=c(a1=1,a2=1,a3=1,a4=1,a5=1,b1=0,b2=0,b3=0,b4=0,b5=0),
)
It gives the following error:
*Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1*
Is there anyone who have an idea what I am doing wrong? Is there any error
in the syntax? I never used nlme package before, so I might be missing some
components in the code. Or, is this just estimation problem and nlme() can
not fit this model for some reason?
I would appreciate any help.
Thanks
Re: [R] Help
You've been told how to do what you ask. But I'm not convinced that you really want to do what you asked. It might be better to do whatever you want with the data leaving it all in one object. There are many ways of doing that, the 'by' function is one of them. On 11/11/2011 20:24, Francesca wrote: Dear Contributors I would like to perform this operation using a loop, instead of repeating the same operation many times. The numbers from 1 to 4 related to different groups that are in the database and for which I have the same data. x-c(1,3,7) datiP1- datiP[datiP$city ==1,x]; datiP2- datiP[datiP$city ==2,x]; datiP3- datiP[datiP$city ==3,x] datiP4- datiP[datiP$city ==4,x]; -- Patrick Burns pbu...@pburns.seanet.com twitter: @portfolioprobe http://www.portfolioprobe.com/blog http://www.burns-stat.com (home of 'Some hints for the R beginner' and 'The R Inferno') __ 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] Estimating IRT models by using nlme() function
Thanks for your suggestion, Dimitris. I know (and use) all of these R packages for IRT (I am actually writer of one of them, EstCRM). I also use specialized commercial software such as BILOG and MULTILOG. My goal is to be able to estimate IRT model parameters in a nonlinear mixed effects regression framework. Recently, there has been some efforts to do that as I give one reference in m previous e-mail. There are two advantages of using non-linear mixed effects regression for IRT estimation. One of them is that it gives a kind of unified approach, because it is very flexible. You can change the nonlinear model however you like, and try to estimate very different models. Then, see which one fits your data better. So, you don't have to force yourself to fit regular and well know IRT models if you have reasons to fit another model. Second, it allows you to model higher clusters. Most of the time, we collect data from clusters, and the subjects within a cluster are not independent from each other. It's a well known concept in HLM framework, and measured by intra-class correlation coefficient. I have never seen in IRT literature about this issue, and we really don't know how cluster dependencies might affect IRT parameter estimation. So, nonlinear mixed effects regression approach might allow you to model a third (or may be a fourth) cluster above the students, and take the cluster dependencies into account while estimating IRT parameters. It also provides flexibilities for other IRT applications ( such as detecting DIF items). Anyway, my goal is not basically estimating IRT parameters. I can do it in so many different ways by using other softwares. My goal is to do it in a nonlinear mixed effects regression framework, and see how well it does comparing to the traditional estimation methods. I think ICC component is very important and its possible impacts were never adressed in IRT literature for parameter estimation. This is open to investigation. But, I need to figure out this nlme() thing first. I can use SAS, because people already provided all SAS syntax to estimate IRT parameters for all well-known models in nonlinear mized effect approach. But, I want to see it if nlme() can replicate what they did since I am an R fan. 2011/11/12 Dimitris Rizopoulos d.rizopou...@erasmusmc.nl Note that there are several IRT dedicated packages in R, such as the ltm and eRm packages. For more info you may check the Psychometrics Task View at: http://cran.r-project.org/web/**views/Psychometrics.htmlhttp://cran.r-project.org/web/views/Psychometrics.html I hope it helps. Best, Dimitris On 11/11/2011 9:38 PM, Cengiz ZopluoÄlu wrote: Hi, I have a question about estimating IRT models by using nlme, not just rasch model, but also other models. Behavior Research Methods http://www.springerlink.com/**content/1554-351x/http://www.springerlink.com/content/1554-351x/ Volume 37, Number 2http://www.springerlink.com/**content/1554-351x/37/2/http://www.springerlink.com/content/1554-351x/37/2/, 202-218, DOI: 10.3758/BF03192688 Using SAS PROC NLMIXED to fit item response theory models (2005). Ching-Fan Sheuhttps://mail.google.com/**mail/u/0/html/compose/static_** files/goog_851944013https://mail.google.com/mail/u/0/html/compose/static_files/goog_851944013 , Cheng-Te Chenhttps://mail.google.com/**mail/u/0/html/compose/static_** files/goog_851944013https://mail.google.com/mail/u/0/html/compose/static_files/goog_851944013 , Ya-Hui Su and Wen-Chung Wang An example of this is provided in the paper above. They use SAS PROC NLMIXED, and estimate all possible binary and Polytomous IRT models. I want to replicate what they did using R and nlme package. I am able to fit rasch model, but have some error messages for 2PL IRT model. So, I could not go beyond that. This is really important for me, because any nonlinear model can be fitted by using this approach. I would like to see how well nlme() package recovers true parameters, not just regular IRT models but also some other less used IRT models. Here is the example. Let's say I have the following dataset. 1000 people responds five items. ID Item Response d1 d2 d3 d4 d5 1.1 110 1 0 0 0 0 1.2 120 0 1 0 0 0 1.3 130 0 0 1 0 0 1.4 141 0 0 0 1 0 1.5 150 0 0 0 0 1 2.1 211 1 0 0 0 0 2.2 220 0 1 0 0 0 2.3 230 0 0 1 0 0 2.4 241 0 0 0 1 0 2.5 250 0 0 0 0 1 3.1 310 1 0 0 0 0 3.2 320 0 1 0 0 0 3.3 331 0 0 1 0 0 3.4 341 0 0 0 1 0 3.5 350 0 0 0 0 1 .. .. 1000.1 100011 1 0 0 0 0
Re: [R] need help
Supreet kaur [ksupre...@gmail.com]
how do I calculate the reliability between the two groups
using the ICCs?
Looking for the irr (inter-rater reliability package) might help; as the name
suggests, it includes a variety of measures of rater reliability.
S Ellison
***
This email and any attachments are confidential. Any use...{{dropped:8}}
__
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] plotting a bivariate quadratic curve
Hi, I would like to plot a x,y curve described by the equation : Ax^2 + Bx + Cy^2 + Dy + E == 0 (A,B,C,D,E are constants) This sounds like quite the common task but haven't been able to figure out how to do this. Could you please help ? I am new to R so probably missing something basic. Thanks __ 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] Second-order effect in Parametric Survival Analysis
Hi experts, http://r.789695.n4.nabble.com/file/n4034318/Parametric_survival_analysis_2nd-order_efffect.JPG Parametric_survival_analysis_2nd-order_efffect.JPG As we know a normal survival regression is the equation (1) Well, I'ld like to modify it to be 2nd-order interaction model as shown in equation(2) Question: Assume a and z is two covariates. x = dummy variable (1 or 0) z = factors (peoples' name) fit - survreg(Surv(time,censor)~x*z, data=sample, dist=exponential) I tried to apply survreg(), while I have few questions: 1) If */survreg(Surv(time,censor)~x*z, data=sample, dist=exponential)/* correct? 2) If the baseline hazard is the value excluded both x and z effects? 3) How can I get the value and plot the hazard with only x effect (but exclude effect z) Thanks Best, Ryusuke -- View this message in context: http://r.789695.n4.nabble.com/Second-order-effect-in-Parametric-Survival-Analysis-tp4034318p4034318.html Sent from the R help mailing list archive at Nabble.com. __ 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 using lmer
Hi can we perform logistic regression using lmer function. Please help me with this -- View this message in context: http://r.789695.n4.nabble.com/logistic-Regression-using-lmer-tp4034056p4034056.html Sent from the R help mailing list archive at Nabble.com. __ 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] Please Help
HiI want to construct a logliikelood function in RHere is the situationy=number
of particles emitted in 1 hr period~pois(30)p=probability of detection of
radiation particlesx=number of particles detected by a radiation
detector~pois(30p)where p~beta(a,1)I have to calculate the loglikehood for a
for the range a(2,50)I wish to simulate 100 random samples for each aHere is my
code:-m=481n=100x = c(15, 36, 29, 28, 37, 32, 25, 27, 31, 21, 25, 27, 28, 31,
28, 20, 34, 25, 20, 34, 15,21, 28, 24, 31, 19, 34, 29, 18, 25, 16, 19, 44, 26,
34, 31, 21, 28, 11, 31, 21, 34, 25, 25,30, 23, 21, 35, 36, 21, 27, 29, 30, 22,
25, 30, 24, 27, 28, 22, 36, 29, 33, 35, 30, 32, 27,26, 25, 27, 23, 21, 39, 33,
24, 21, 19, 34, 32, 28, 27, 28, 23, 20, 24, 29, 21, 22, 31, 28,27, 28, 29, 21,
30, 28, 31, 22, 29, 18)a=c(1:m)*0.1+1.9p-rbeta(n,a,1)loglik-numeric(m)for(i
in
1:n){p[i]-rbeta(n,a[i],1)x[i]-rpois(n,p[i]*30)loglik-function(x,p)logLik(a=a(i))-n*log(a)+(log(30)*sum(x))-sum(log(factorial(x)))-30*(!
sum(p[i]))+sum(x[i]*log(p[i]))+(a[j]-1)*sum(log(p[i]))}But somehow I am not
getting the answer please HELPP.S. the data for x is a random sample of number
of particles detected in 1 hr (given in the question)Really appreciate the
helpThanks
[[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] Subsetting data leads to funky plots
Hi, On Sat, Nov 12, 2011 at 2:04 AM, Vinny Moriarty vwmoria...@gmail.com wrote: I'm trying out a basic plot, but something about the way I subset my data leads to problems with the plot. Here is the first bit of my data set year,date,location,quadrat_juvenile,photo_location,photo_exists,genus,count,divers 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Acanthastrea,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Acropora,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Astreopora,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Cyphastrea,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Favia,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Fungia,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Gardinoseris,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Herpolitha,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Leptastrea,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Leptoseris,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Lobophyllia,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Millepora,0,HP+MEM 2005,2005-04-30, 1 Fringing Reef,1, 1 Fringing Reef Coral Transect Pole 1-2 Quadrat 1,t,Montastrea,0,HP+MEM dput() with a sample of your data would be useful, so that we ca try it. Among other things, the example you give above doesn't *have* the values that you're subsetting on, so we can't easily try your specific code. str() of your data would also be useful. Are these factors? Strings? The small reproducible example is important. I need to breakdown the data before I can plot it, so this is the stepwise code I use: Juv=read.csv(juvenile_density_20110506.csv) Juv$count[Juv$count0]- NA #Missing data is imputed as negative data, so first label it as NA JuvFor=subset(Juv,(location== 2 Outer 10 m) | (location== 1 Outer 10 m)) #Subset the data of interest Juv_Sum_by_quad=aggregate(count~year+location+quadrat_juvenile,data= JuvFor,sum) # Calculate sum of each quadrat Juv_Avg=aggregate(count~year+location,data=Juv_Sum_by_quad, mean) #Calculate yearly means So far so good. Is it? Did you look at Juv_Avg to make sure that it contains what you expect? I thought I could do this: plot(Juv_Sum1$year, Juv_Sum1$count,type=L) There's no type=L, although there is type=l Is Juv_Sum1 sorted by year? Without a reproducible example, I can't see which of the many things that could go wrong have gone wrong. R does exactly what you tell it, and it sounds like you may have given it data that weren't sorted. For that matter, there's no Juv_Sum1 in the code you gave, so I don't have the slightest idea what might be in it. But not only do I get odd separate lines as opposed to a time series line plot, but I the first tick mark on the x-axis lists the amount of rows in my data. I can fix it by saving the final object as a .csv and reloading it the data- but that defeats the whole purpose. Can anyone explain why I get this plot result? Looking through past posts I found that I can get away with just plotting the y-axis. And if I modify the x axis separately I can get the plot I'm looking for. plot(Juv_Sum1$count,xaxt=n) axis(1,at=1:length(Juv_Sum1$year),labels=Juv_Sum1$year) But I'm concerned about why my data subsets are giving me odd results in the first place, especially as I move on into more complicated bits of code. This list is generally quite helpful, but not telepathic. See these two lines from every message (below)? That's what we need. PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Sarah -- Sarah Goslee http://www.functionaldiversity.org __ 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] processing data of hierarchical river network
Hi, I have some hydrochemical data from observation stations of a river network. The data set provides information about the river hierarchy, i.e. the order of the rivers (main river is 1st order, contributing is 2nd order, and contributing to 2nd is 3rd order... etc). Now I would like to use the hierarchical information to: a) plot a dendrogram of the observation stations b) perform some simple statistics on the different branches What packages are needed to conduct that investigation? How the hierarchical information should be formated to be easily processed with R? Thank you, Torsten __ 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] Please Help
Two problems with your post:
a) you used a format that is not accepted by the mail server and
resulted in all of your linefeeds being removed, so it is now
basically unreadable. As the Posting Guide says, please use plain
text. (And please read the Posting Guide for yourself ... and do it
now before reading further.) All the free mailing services provide
mechanisms for constructing plain text emails.
b) it appears to be homework, and homework submission are discouraged.
You can attach a plea describing why you should be exempted from this
rule, but it might require identifying your self more completely and
explaining why you are not using the services of the educational
institution you are studying under.
On Nov 12, 2011, at 8:09 AM, avish kamdar wrote:
HiI want to construct a logliikelood function in RHere is the
situationy=number of particles emitted in 1 hr
period~pois(30)p=probability of detection of radiation
particlesx=number of particles detected by a radiation
detector~pois(30p)where p~beta(a,1)I have to calculate the
loglikehood for a for the range a(2,50)I wish to simulate 100 random
samples for each aHere is my code:-m=481n=100x = c(15, 36, 29, 28,
37, 32, 25, 27, 31, 21, 25, 27, 28, 31, 28, 20, 34, 25, 20, 34,
15,21, 28, 24, 31, 19, 34, 29, 18, 25, 16, 19, 44, 26, 34, 31, 21,
28, 11, 31, 21, 34, 25, 25,30, 23, 21, 35, 36, 21, 27, 29, 30, 22,
25, 30, 24, 27, 28, 22, 36, 29, 33, 35, 30, 32, 27,26, 25, 27, 23,
21, 39, 33, 24, 21, 19, 34, 32, 28, 27, 28, 23, 20, 24, 29, 21, 22,
31, 28,27, 28, 29, 21, 30, 28, 31, 22, 29, 18)a=c(1:m)*0.1+1.9p-
rbeta(n,a,1)loglik-numeric(m)for(i in 1:n){p[i]-rbeta(n,a[i],
1)x[i]-rpois(n,p[i]*30)loglik-function(x,p)logLik(a=a(i))-n*log(a)
+(log(30)*sum(x))-sum(log(factorial(x)))-30*(!
sum(p[i]))+sum(x[i]*log(p[i]))+(a[j]-1)*sum(log(p[i]))}But somehow I
am not getting the answer please HELPP.S. the data for x is a random
sample of number of particles detected in 1 hr (given in the
question)Really appreciate the helpThanks
[[alternative HTML version deleted]]
David Winsemius, MD
West Hartford, CT
__
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] State space model
Hi,
I'm trying to estimate the parameters of a state space model of the
following form
measurement eq:
z_t = a + b*y_t + eps_t
transition eq
y_t+h = (I -exp(-hL))theta + exp(-hL)y_t+ eta_{t+h}.
The problem is that the distribution of the innovations of the transition
equation depend on the previous value of the state variable.
To be exact: y_t|y_{t-1} ~N(mu, Q_t) where Q is a diagonal matrix with
elements equal to
Q_{i,t} =
sigma_i*(1-exp(-kappa_i*h)/kappa_i*(theta_i/2*(1-exp(kappa_i*h)+exp(-kappa_i*h)y_{t-1,i}
The fkf returns the filtered states variables so y_{t-1,i} is available. I
just can't figure out how to write my program in such a way that this
information is included and updated in the state space model for each
iteration in optim. Any suggestions on how to solve this are much
appreciated.
Thank you,
Kristian.
Below my program where the variance matrices are just identity matrices and
the data is just random numbers. I used the example from the FKF package as
framework.
library(FKF) #loading Fast Kalman Filter package
library(Matrix) # matrix exponential package
library(BB)
library(alabama)
x - matrix(abs(rnorm(400)), nrow=10, ncol=40)
m - 2 # m is the number of state variables
n - ncol(x) # is the length of the observed sample
d - nrow(x) # is the number of observed variables.
h - 1/52
## creating state space representation of 2-factor CIR model
CIR2ss - function(K_1, K_2, sigma_1, sigma_2, lambda_1, lambda_2, theta_1,
theta_2, delta_0, delta_1, delta_2) {
## defining auxilary parameters,
phi_11 - sqrt((K_1+lambda_1)^2+2*sigma_1^2*delta_1)
phi_21 - sqrt((K_2+lambda_2)^2+2*sigma_2^2*delta_2)
phi_12 - K_1+lambda_1+phi_11
phi_22 - K_2+lambda_2+phi_21
phi_13 - 2*K_1*theta_1/sigma_1^2
phi_23 - 2*K_2*theta_2/sigma_2^2
a - array(1, c(d,n))
phi_14 - numeric(d)
phi_24 - numeric(d)
b1 - numeric(d)
b2 - numeric(d)
for(t in 1:d){
phi_14[t] - 2*phi_11+phi_12*(exp(phi_11*t)-1)
phi_24[t] - 2*phi_21+phi_22*(exp(phi_21*t)-1)
a[t] - delta_0 - phi_13/(t)*log(2*phi_11*exp(phi_12*(t)/2)/phi_14[t])-
phi_23/(t)*log(2*phi_21*exp(phi_22*(t)/2)/phi_24[t])
b1[t] - (delta_1/(t))*2*(exp(phi_11*(t))-1)/phi_14[t]
b2[t] - (delta_2/(t))*2*(exp(phi_21*(t))-1)/phi_24[t]
}
b - array(c(b1,b2), c(d,m,n))
j - -matrix(c(K_1, 0, 0, K_2), c(2,2))*h
explh - as.matrix(expm(j))
Tt - array(explh, c(m,m,n)) #array giving the factor of the transition
equation
Zt - b #array giving the factor of the measurement equation
ct - a #matrix giving the intercept of the measurement equation
dt - as.matrix((diag(m)-explh)%*%c(theta_1, theta_2)) #matrix giving
the intercept of the transition equation
GGt - array(diag(d), c(d,d,1)) #array giving the variance of the
disturbances of the measurement equation
HHt - diag(m) #array giving the variance of the innovations of the
transition equation
a0 - c(0.5, 0.5) #vector giving the initial value/estimation of the
state variable
P0 - matrix(1e6, nrow = 2, ncol = 2) # matrix giving the variance of a0
return(list(a0 = a0, P0 = P0, ct = ct, dt = dt, Zt = Zt, Tt = Tt, GGt =
GGt,
HHt = HHt))
}
## Objective function passed to optim
objective - function(parm, yt) {
sp - CIR2ss(parm[K_1], parm[K_2], parm[sigma_1], parm[sigma_2],
parm[lambda_1], parm[lambda_2],
parm[theta_1], parm[theta_2], parm[delta_0],
parm[delta_1], parm[delta_2])
ans - fkf(a0 = sp$a0, P0 = sp$P0, dt = sp$dt, ct = sp$ct, Tt = sp$Tt,
Zt = sp$Zt, HHt = sp$HHt, GGt = sp$GGt, yt = yt)
return(-ans$logLik)
}
parm - c(K_1 = 0.6048, K_2 = 0.656, sigma_1 =0.1855, sigma_2 =0.5524,
lambda_1 =0.55, lambda_2 =0.009, theta_1 = 0.173, theta_2 = 0.12,
delta_0 =0.686, delta_1 =-0.003, delta_2=0.0025)
hin - function(parm, yt){
h- numeric(2)
h[1] - 2*parm[1]*parm[7]-parm[3]^2
h[2] - 2*parm[2]*parm[8]-parm[4]^2
h
}
##optimizing objective function
fit - auglag(par = parm, fn = objective, yt = x, control.outer =
list(method = CG), hin = hin)
print(fit)
[[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] logistic Regression using lmer
arunkumar akpbond007 at gmail.com writes: Hi can we perform logistic regression using lmer function. Please help me with this Yes. library(lme4) glmer([reponse]~[fixed effects (covariates)]+(1|[grouping variable]), data=[data frame], family=binomial) Further questions on this topic should probably go to the r-sig-mixed-models mailing list. __ 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] plotting a bivariate quadratic curve
On 11-11-12 5:50 AM, Yackov Lubarsky wrote: Hi, I would like to plot a x,y curve described by the equation : Ax^2 + Bx + Cy^2 + Dy + E == 0 (A,B,C,D,E are constants) This sounds like quite the common task but haven't been able to figure out how to do this. Could you please help ? I am new to R so probably missing something basic. R is not very good at that, but one way to do it is to draw a contour plot of the bivariate function, with a single contour at level 0. For example, A - B - C - D - 1 E - -1 x - y - seq(-3, 3, len=100) z - outer(x, y, function(x, y) A*x^2 + B*x + C*y^2 + D*y + E) contour(x, y, z, levels=0, drawLabels=FALSE) Another would be to solve for y as a function of x, and draw both branches. Duncan Murdoch __ 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] plotting a bivariate quadratic curve
On Sat, Nov 12, 2011 at 10:50 AM, Yackov Lubarsky yluba...@gmail.com wrote: Hi, I would like to plot a x,y curve described by the equation : Ax^2 + Bx + Cy^2 + Dy + E == 0 (A,B,C,D,E are constants) This sounds like quite the common task but haven't been able to figure out how to do this. Could you please help ? I am new to R so probably missing something basic. What you may be missing was all figured out by the ancient Greeks and they didnt even have R in their alphabet! You've got this: http://mathworld.wolfram.com/QuadraticCurve.html but with b=0. I'm not sure how many of the awkward cases (line pairs, imaginary everything) disappear with b=0, but you can compute the Delta, I, J, and K determinants to figure it out. Then set theta to seq(0,2*pi,len=100) and compute r from the polar equation form. Plot r*cos(theta), r*sin(theta)... Confession: its been a long time since I did maths like this (forgive me father, for I have not sin'd). Barry __ 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] Passing parameters
Hi Annie, The error is in the R script...you start reading with Args[1] but instead you should start with Args[5] etc. because (I believe) that R --slave --vanilla --args count as [1] [2][3][4] Therefore, Args[5] it's the first argument you use...try that and let me know!! I am trying something similar for one of my scripts Francesco -- View this message in context: http://r.789695.n4.nabble.com/Passing-parameters-tp3762480p4034548.html Sent from the R help mailing list archive at Nabble.com. __ 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] With an example - Re: rbind.data.frame drops attributes for factor variables
When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. See the sample example below to reproduce the problem. Please suggest How I can fix this. Thanks, Sammy a=c(Male, Male, Female, Male) b=c(1,2,3,4) c=c(great, bad, good, bad) dataset- data.frame (gender = a, count = b, answer = c) dataset gender count answer 1 Male 1 great 2 Male 2bad 3 Female 3 good 4 Male 4bad attributes(dataset$answer) $levels [1] bad good great $class [1] factor Now adding some custom attributes to column dataset$answer attributes(dataset$answer)-c(attributes(dataset$answer),list(newattr1=custom-attr1)) attributes(dataset$answer)-c(attributes(dataset$answer),list(newattr2=custom-attr2)) attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 However as soon as I add a row to this data frame (dataset) by rbind(), it loses the custom attributes (newattr1 and newattr2) I have just added newrow = c(gender=Female, count = 5, answer = great) dataset - rbind(dataset, newrow) attributes(dataset$answer) $levels [1] bad good great $class [1] factor the two custom attributes are dropped!! Any suggestion why this is happening. On Fri, Nov 11, 2011 at 11:44 AM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote: As the doctor says, if it hurts don't do that. A factor is a sequence of integers with a corresponding list of character strings. Factors in two separate vectors can and usually do map the same integer to different strings, and R cannot tell how you want that resolved. Convert these columns to character before combining them, and only convert to factor when you have all of your possibilities present (or you specify them in the creation of the factor vector). --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Sammy Zee szee2...@gmail.com wrote: Hi all, When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. I see the following post with same problem. However i did not see any reply to the following posting offering a solution. Could someone please help. http://r.789695.n4.nabble.com/rbind-data-frame-drops-attributes-for-factor-variables-td919575.html Thanks, Sammy [[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. [[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] changelog for MASS?
Does anyone know where I can find a changelog for MASS? It's difficult to know whether I should ask my company to update the package or not. We are usually required to show the changelog, but I can't find one. Thank you, Xu -- View this message in context: http://r.789695.n4.nabble.com/changelog-for-MASS-tp4034473p4034473.html Sent from the R help mailing list archive at Nabble.com. __ 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] R v2.13.2 - Cannot find Rcmd on path?
Dear all: I could not able to find rcmd using the following path: either R version (2.13.2 or 2.14.0); either rtools version (2.13 or 2.14). My os is winxp. The variables set in path are as follows. c:\Rtools\bin;c:\Rtools\MinGW\bin;C:\Perl\bin; C:\Perl\site\bin; c:\Rtools\bin; C:\Program Files\R\R-2.13.2\bin\i386; C:\Program Files\MiKTeX 2.9\miktex\bin;%SystemRoot%\system32;%SystemRoot%;%SystemRoot%\System32\Wbem;C:\Program Files\Common Files\Roxio Shared\DLLShared\;C:\Program Files\Common Files\Roxio Shared\10.0\DLLShared\;c:\Program Files\Microsoft SQL Server\90\Tools\binn\ After I tried setpath.bat file, it doesn't work either. Thank you, Jixiang Wu -- View this message in context: http://r.789695.n4.nabble.com/R-v2-13-2-Cannot-find-Rcmd-on-path-tp3927126p4034596.html Sent from the R help mailing list archive at Nabble.com. __ 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] list.dir() function
firstly - it looks like you define orgDir but then your code later refers to indexDir? that being said, I don't think that is relevant. You could always try upgrading to the latest version and see if the problem goes away (if you're not already). I can't replicate the issue here but we can't guess what version you're running of R. It's always helpful to include What build/operating system etc. Use the following command and copy and paste. into your reply please. R.version One option you've got perhaps is to indexDir into the directory and then try the same commands you've mentioned in your post. another option which might work would be based on these two lines. a1-file.info(list.files(indexDir)) mydir-row.names(a1[a1$isdir==TRUE,]) let me know if you can't get it. I don't know why your code doesn't work. -- View this message in context: http://r.789695.n4.nabble.com/list-dir-function-tp4033044p4035093.html Sent from the R help mailing list archive at Nabble.com. __ 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] changelog for MASS?
It's on CRAN, where you can also see the record of updates: http://cran.r-project.org/web/packages/MASS/index.html Specifically, http://cran.r-project.org/web/packages/MASS/NEWS Hope this helps On Sat, Nov 12, 2011 at 8:55 AM, Xu Wang xuwang...@gmail.com wrote: Does anyone know where I can find a changelog for MASS? It's difficult to know whether I should ask my company to update the package or not. We are usually required to show the changelog, but I can't find one. Thank you, Xu -- View this message in context: http://r.789695.n4.nabble.com/changelog-for-MASS-tp4034473p4034473.html Sent from the R help mailing list archive at Nabble.com. __ 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.
Re: [R] With an example - Re: rbind.data.frame drops attributes for factor variables
On Nov 12, 2011, at 2:47 PM, Sammy Zee wrote: When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. See the sample example below to reproduce the problem. Please suggest How I can fix this. Thanks, Sammy a=c(Male, Male, Female, Male) b=c(1,2,3,4) c=c(great, bad, good, bad) dataset- data.frame (gender = a, count = b, answer = c) dataset gender count answer 1 Male 1 great 2 Male 2bad 3 Female 3 good 4 Male 4bad attributes(dataset$answer) $levels [1] bad good great $class [1] factor Now adding some custom attributes to column dataset$answer attributes(dataset$answer)-c(attributes(dataset $answer),list(newattr1=custom-attr1)) attributes(dataset$answer)-c(attributes(dataset $answer),list(newattr2=custom-attr2)) If you look through the code of rbind.data.frame you see that column values are processed with the 'factor' function. attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 attributes(factor(dataset$answer)) $levels [1] bad good great $class [1] factor So I think you are out of luck. You will need to restore the special attributes yourself. -- David. attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 However as soon as I add a row to this data frame (dataset) by rbind(), it loses the custom attributes (newattr1 and newattr2) I have just added newrow = c(gender=Female, count = 5, answer = great) dataset - rbind(dataset, newrow) attributes(dataset$answer) $levels [1] bad good great $class [1] factor the two custom attributes are dropped!! Any suggestion why this is happening. On Fri, Nov 11, 2011 at 11:44 AM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote: As the doctor says, if it hurts don't do that. A factor is a sequence of integers with a corresponding list of character strings. Factors in two separate vectors can and usually do map the same integer to different strings, and R cannot tell how you want that resolved. Convert these columns to character before combining them, and only convert to factor when you have all of your possibilities present (or you specify them in the creation of the factor vector). --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Sammy Zee szee2...@gmail.com wrote: Hi all, When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. I see the following post with same problem. However i did not see any reply to the following posting offering a solution. Could someone please help. http://r.789695.n4.nabble.com/rbind-data-frame-drops-attributes-for-factor-variables-td919575.html Thanks, Sammy [[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. [[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. David Winsemius, MD West Hartford, CT __ 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] Second-order effect in Parametric Survival Analysis
On Nov 12, 2011, at 7:37 AM, ryusuke wrote: Hi experts, http://r.789695.n4.nabble.com/file/n4034318/Parametric_survival_analysis_2nd-order_efffect.JPG Parametric_survival_analysis_2nd-order_efffect.JPG As we know a normal survival regression is the equation (1) Well, I'ld like to modify it to be 2nd-order interaction model as shown in equation(2) Question: Assume a and z is two covariates. x = dummy variable (1 or 0) z = factors (peoples' name) fit - survreg(Surv(time,censor)~x*z, data=sample, dist=exponential) I tried to apply survreg(), while I have few questions: 1) If */survreg(Surv(time,censor)~x*z, data=sample, dist=exponential)/* correct? The formula interface for R would expand x*z to x + z + x:z (Which is not the formula in your Nabble-provided-jpg, but from your later questions is probably what you want anyway.) 2) If the baseline hazard is the value excluded both x and z effects? Maybe. You won't be excluding them so much as holding their values jointly at zero, which may or may not be the same thing. 3) How can I get the value and plot the hazard with only x effect (but exclude effect z) You will never be able to do so. If you have an interacting variable in a model, there will always be an effect of that covariate on predictions associated with any covariate with which it is interacting. You should be able to display or plot the x- effects (note the plural) that are estimated for chosen levels of z, however. To accomplish that you should construct an appropriate data.frame and offer it as the newdata argument to predict(fit) just as you would do with any properly constructed R/S regression package. There is a worked example on this posting from the Master: http://finzi.psych.upenn.edu/Rhelp10/2010-May/240458.html -- David Winsemius, MD West Hartford, CT __ 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] dev.new() within a loop
Hello,
I have a loop where I iterate performance data files within a folder, parse and
plot them in one shot (see below).
However, when executing plot_raw which invokes dev.new(..) all windows come out
blank whereas if I execute each file outside of a loop then I can see the plots
properly. What's wrong here?
Thanks in advance,
Best regards,
Giovanni
# given a directory name, it will iterate all files that match the given pattern
#basedir - /Users/bravegag/code/asl11/data/2k-r1-test-2011_data/
basedir - /Users/bravegag/code/asl11/data/nclients_2_128-2010_data/
pattern - paste(logs.*cl\\-.*mw\\-.*db\\-.*\\-client\\.dat,sep=)
all_files - dir(path=basedir, pattern=pattern)
throughput - NULL
response - NULL
#file_name - all_files[1]
# iterate all files
for (file_name in all_files) {
print(paste(processing, file_name, ...))
df - read.table(paste(basedir, file_name, sep=)) # read
the data as a data frame
names(df)[names(df)==V1] - Time
names(df)[names(df)==V2] - Partitioning
names(df)[names(df)==V3] - Workload
names(df)[names(df)==V4] - Runtime
# get rid of first and last n minutes
df - subset(df, df$Time warmup_cooldown_minutes)
df - subset(df, df$Time (max(df$Time) - warmup_cooldown_minutes))
#
=
# Throughput
#
=
if (decouple) {
dft - aggregate(x=df$Runtime, by=list(df$Time,df$Workload),
FUN=length)
names(dft)[names(dft)==Group.1] - Time
names(dft)[names(dft)==Group.2] - Workload
names(dft)[names(dft)==x] - Y
} else {
dft - aggregate(x=df$Runtime, by=list(df$Time), FUN=length)
names(dft)[names(dft)==Group.1] - Time
names(dft)[names(dft)==x] - Y
}
dft$se - 0
plot_raw(dft,connect=TRUE,label=file_name)
}
[[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] With an example - Re: rbind.data.frame drops attributes for factor variables
Thanks David. Besides rbind(), is there any other way to add a row to a data frame so that I do not lose the custom attributes. Thanks, Sammy On Sat, Nov 12, 2011 at 5:17 PM, David Winsemius dwinsem...@comcast.netwrote: On Nov 12, 2011, at 2:47 PM, Sammy Zee wrote: When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. See the sample example below to reproduce the problem. Please suggest How I can fix this. Thanks, Sammy a=c(Male, Male, Female, Male) b=c(1,2,3,4) c=c(great, bad, good, bad) dataset- data.frame (gender = a, count = b, answer = c) dataset gender count answer 1 Male 1 great 2 Male 2bad 3 Female 3 good 4 Male 4bad attributes(dataset$answer) $levels [1] bad good great $class [1] factor Now adding some custom attributes to column dataset$answer attributes(dataset$answer)-c(**attributes(dataset$answer),** list(newattr1=custom-attr1)) attributes(dataset$answer)-c(**attributes(dataset$answer),** list(newattr2=custom-attr2)) If you look through the code of rbind.data.frame you see that column values are processed with the 'factor' function. attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 attributes(factor(dataset$**answer)) $levels [1] bad good great $class [1] factor So I think you are out of luck. You will need to restore the special attributes yourself. -- David. attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 However as soon as I add a row to this data frame (dataset) by rbind(), it loses the custom attributes (newattr1 and newattr2) I have just added newrow = c(gender=Female, count = 5, answer = great) dataset - rbind(dataset, newrow) attributes(dataset$answer) $levels [1] bad good great $class [1] factor the two custom attributes are dropped!! Any suggestion why this is happening. On Fri, Nov 11, 2011 at 11:44 AM, Jeff Newmiller jdnew...@dcn.davis.ca.us**wrote: As the doctor says, if it hurts don't do that. A factor is a sequence of integers with a corresponding list of character strings. Factors in two separate vectors can and usually do map the same integer to different strings, and R cannot tell how you want that resolved. Convert these columns to character before combining them, and only convert to factor when you have all of your possibilities present (or you specify them in the creation of the factor vector). --**--** --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --**--** --- Sent from my phone. Please excuse my brevity. Sammy Zee szee2...@gmail.com wrote: Hi all, When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. I see the following post with same problem. However i did not see any reply to the following posting offering a solution. Could someone please help. http://r.789695.n4.nabble.com/**rbind-data-frame-drops-** attributes-for-factor-**variables-td919575.htmlhttp://r.789695.n4.nabble.com/rbind-data-frame-drops-attributes-for-factor-variables-td919575.html Thanks, Sammy [[alternative HTML version deleted]] __** R-help@r-project.org mailing list https://stat.ethz.ch/mailman/**listinfo/r-helphttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/**posting-guide.htmlhttp://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-helphttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/** posting-guide.html http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD West Hartford, CT [[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
Re: [R] efpFunctional construction (strucchange package)
Unfortunately, I couldn't find neither this source file, nor function wmax(). Was it any old version of strucchange? There are only files strucchange.R, strucchange.rdb and strucchange.rdx in directory \strucchange\R. Best, Julia -- View this message in context: http://r.789695.n4.nabble.com/efpFunctional-construction-strucchange-package-tp4023903p4035497.html Sent from the R help mailing list archive at Nabble.com. __ 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] Simulation over data repeatedly for four loops
Dear Contributors,
I am trying to perform a simulation over sample data,
but I need to reproduce the same simulation over 4 groups of data. My
ability with for loop is null, in particular related
to dimensions as I always get, no matter what I try,
number of items to replace is not a multiple of replacement length
This is what I intend to do: replicate this operation for
four times, where the index for the four groups is in the
part of the code: datiPc[[1]][,2].
I have to replicate the following code 4 times, where the
changing part is in the data from which I pick the sample,
the data that are stored in datiPc[[1]][,2].
If I had to use data for the four samples, I would substitute the 1 with a
j and replicate a loop four times, but it never worked.
My desired final outcome is a matrix with 1 observations for each
couple of extracted samples, i.e. 8 columns of 1 observations of means.
db-c()
# Estrazione dei campioni dai dati di PGG e TRUST
estr1 - c();
estr2 - c();
m1-c()
m2-c()
tmp1- data1[[1]][,2];
tmp2- data2[[2]][,2];
for(i in 1:100){
estr1-sample(tmp1, 1000, replace = TRUE)
estr2-sample(tmp2, 1000, replace = TRUE)
m1[i]-mean(estr1,na.rm=TRUE)
m2[i]-mean(estr2,na.rm=TRUE)
}
db-data.frame(cbind(m1,m2))
Thanks for any help you can provide.
Best Regards
--
Francesca
--
[[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] With an example - Re: rbind.data.frame drops attributes for factor variables
On Nov 12, 2011, at 6:40 PM, Sammy Zee wrote: Thanks David. Besides rbind(), is there any other way to add a row to a data frame so that I do not lose the custom attributes. I have already told you the method that I know of. You don't seem to have taken my poin that it is not a data.frame specific problem but rahter a facor problem. You are welcome to redefine `rbind.data.frame`. The R language is rather flexible in that manner. -- David. Thanks, Sammy On Sat, Nov 12, 2011 at 5:17 PM, David Winsemius dwinsem...@comcast.net wrote: On Nov 12, 2011, at 2:47 PM, Sammy Zee wrote: When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. See the sample example below to reproduce the problem. Please suggest How I can fix this. Thanks, Sammy a=c(Male, Male, Female, Male) b=c(1,2,3,4) c=c(great, bad, good, bad) dataset- data.frame (gender = a, count = b, answer = c) dataset gender count answer 1 Male 1 great 2 Male 2bad 3 Female 3 good 4 Male 4bad attributes(dataset$answer) $levels [1] bad good great $class [1] factor Now adding some custom attributes to column dataset$answer attributes(dataset$answer)-c(attributes(dataset $answer),list(newattr1=custom-attr1)) attributes(dataset$answer)-c(attributes(dataset $answer),list(newattr2=custom-attr2)) If you look through the code of rbind.data.frame you see that column values are processed with the 'factor' function. attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 attributes(factor(dataset$answer)) $levels [1] bad good great $class [1] factor So I think you are out of luck. You will need to restore the special attributes yourself. -- David. attributes(dataset$answer) $levels [1] bad good great $class [1] factor $newattr1 [1] custom-attr1 $newattr2 [1] custom-attr2 However as soon as I add a row to this data frame (dataset) by rbind(), it loses the custom attributes (newattr1 and newattr2) I have just added newrow = c(gender=Female, count = 5, answer = great) dataset - rbind(dataset, newrow) attributes(dataset$answer) $levels [1] bad good great $class [1] factor the two custom attributes are dropped!! Any suggestion why this is happening. On Fri, Nov 11, 2011 at 11:44 AM, Jeff Newmiller jdnew...@dcn.davis.ca.uswrote: As the doctor says, if it hurts don't do that. A factor is a sequence of integers with a corresponding list of character strings. Factors in two separate vectors can and usually do map the same integer to different strings, and R cannot tell how you want that resolved. Convert these columns to character before combining them, and only convert to factor when you have all of your possibilities present (or you specify them in the creation of the factor vector). --- Jeff NewmillerThe . . Go Live... Sammy Zee szee2...@gmail.com wrote: Hi all, When I use rbind() or rbind.data.frame() to add a row to an existing dataframe, it appears that attributes for the column of type factor are dropped. I see the following post with same problem. However i did not see any reply to the following posting offering a solution. Could someone please help. http://r.789695.n4.nabble.com/rbind-data-frame-drops-attributes-for-factor-variables-td919575.html Thanks, Sammy [[alternative HTML version deleted]] ___ David Winsemius, MD West Hartford, CT __ 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] dev.new() within a loop
On Nov 12, 2011, at 6:04 PM, Giovanni Azua wrote:
Hello,
I have a loop where I iterate performance data files within a
folder, parse and plot them in one shot (see below).
However, when executing plot_raw which invokes dev.new(..) all
windows come out blank whereas if I execute each file outside of a
loop then I can see the plots properly.
Perhaps ...(you did not say what package this plot_raw function comes
from) ... Read the FAQ about why lattice plot don't print. (It
applies to all grid based plotting functions.)
--
David
What's wrong here?
Thanks in advance,
Best regards,
Giovanni
# given a directory name, it will iterate all files that match the
given pattern
#basedir - /Users/bravegag/code/asl11/data/2k-r1-
test-2011_data/
basedir - /Users/bravegag/code/asl11/data/
nclients_2_128-2010_data/
pattern - paste(logs.*cl\\-.*mw\\-.*db\\-.*\\-client\\.dat,sep=)
all_files - dir(path=basedir, pattern=pattern)
throughput - NULL
response - NULL
#file_name - all_files[1]
# iterate all files
for (file_name in all_files) {
print(paste(processing, file_name, ...))
df - read.table(paste(basedir, file_name, sep=)) #
read the data as a data frame
names(df)[names(df)==V1] - Time
names(df)[names(df)==V2] - Partitioning
names(df)[names(df)==V3] - Workload
names(df)[names(df)==V4] - Runtime
# get rid of first and last n minutes
df - subset(df, df$Time warmup_cooldown_minutes)
df - subset(df, df$Time (max(df$Time) - warmup_cooldown_minutes))
#
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
==
# Throughput
#
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
==
if (decouple) {
dft - aggregate(x=df$Runtime, by=list(df$Time,df$Workload),
FUN=length)
names(dft)[names(dft)==Group.1] - Time
names(dft)[names(dft)==Group.2] - Workload
names(dft)[names(dft)==x] - Y
} else {
dft - aggregate(x=df$Runtime, by=list(df$Time), FUN=length)
names(dft)[names(dft)==Group.1] - Time
names(dft)[names(dft)==x] - Y
}
dft$se - 0
plot_raw(dft,connect=TRUE,label=file_name)
}
[[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.
David Winsemius, MD
West Hartford, CT
__
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] Myriam Saavedra M. Sc. Questions about maximun radius distance
Dear Mr. Baddeley I just graduated from a Masters in Applied Mathematics on Jun19th. My thesis was about spatial distribution /a nalysis of some trees in a part of the Congo Basic Forest. In my thesis I used your spatial package in R, and today I'm doing a more deeper study about how we choise the r distance in Function F(). I would like to be able to understand about value of rmaxdefault as: ripley - min(diff(W$xrange), diff(W$yrange))/4 rlarge - if (!missing(lambda)) sqrt(1000/(pi * lambda)) else Inf rmax - min(rlarge, ripley) For the ripley's calculation, I found your explanation in the internet but for the rlarge I couldn't find it. Could you explain why it is using the value inside of the sqrt (1000/(pi * lambda)). Many thanks in advance if you do have the time to answer, Myriam Saavedra Universidad San Francisco de Quito ... Adress Home 3808 Lillooet St. VANCOUVER BC CANADA, V5R-2E6 [[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] identify duplicate from more than one column
Hi all, I've searched everywhere to try to find out how to do this and have had no luck. I am trying to construct identifiers for couples in a dataset. Essentially, I want to identify couples using more than one column as identifiers. Take for instance: obs unithome zsex age 1 015029 18 11 053 2 015029 18 12 049 3 015029 01 11 038 4 015029 01 12 033 5 015029 02 11 036 6 015029 02 12 033 7 015029 03 11 023 8 015029 03 12 019 9 015029 04 12 045 10 015029 05 12 047 Where unit is the housing unit, home is household. Of course, there are more values for unit, although these first ten observations consist of the same unit (which could possibly be an apartment complex). Nonetheless, I want to construct an identifier for couples if unit, home match, but only if both male and female are within the same household. Taking the example data above, I want to see this: unithomez sex age couple 1 015029 18 11 053 1 2 015029 18 12 049 1 3 015029 01 11 038 2 4 015029 01 12 033 2 5 015029 02 11 036 3 6 015029 02 12 033 3 7 015029 03 11 023 4 8 015029 03 12 019 4 9 015029 04 12 045 0 10 015029 05 12 047 0 As you can see in the last two observations, there were no males identified within the same household, thus the last two observations would not contain couple identifiers, rather some other identifier (but the same one) so I can detect them and remove them later. I've tried using the duplicated function but was not very useful. Any help would be greatly appreciated!!! Thanks, Carlos -- View this message in context: http://r.789695.n4.nabble.com/identify-duplicate-from-more-than-one-column-tp4035888p4035888.html Sent from the R help mailing list archive at Nabble.com. __ 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] identify duplicate from more than one column
Hi Carlos,
Here is one option:
## read in your data
dat - read.table(textConnection(
obs unithome zsex age
1 015029 18 11 053
2 015029 18 12 049
3 015029 01 11 038
4 015029 01 12 033
5 015029 02 11 036
6 015029 02 12 033
7 015029 03 11 023
8 015029 03 12 019
9 015029 04 12 045
10 015029 05 12 047),
header = TRUE, stringsAsFactors = FALSE)
closeAllConnections()
## create a unique ID for matching unit and home
dat$mID - with(dat, paste(unit, home, sep = ''))
## somewhat messy way of creating a couple number
## for each mID, if there is more than 1 row, and more than 1 sex
## it creates a couple id, otherwise 0
i - 0L
dat$couple - with(dat, unlist(lapply(split(sex, mID), function(x) {
i - i + 1L
if (length(x) 1 length(unique(x)) 1) {
rep(i, length(x))
} else 0L
})))
## view results
dat
obs unit home z sex age mID couple
11 15029 18 1 1 53 1502918 1
22 15029 18 1 2 49 1502918 1
33 150291 1 1 38 150291 2
44 150291 1 2 33 150291 2
55 150292 1 1 36 150292 3
66 150292 1 2 33 150292 3
77 150293 1 1 23 150293 4
88 150293 1 2 19 150293 4
99 150294 1 2 45 150294 0
10 10 150295 1 2 47 150295 0
See these functions for more details:
?ave # where I got my idea
?split
?lapply
?`-`
Cheers,
Josh
On Sat, Nov 12, 2011 at 8:16 PM, jour4life jour4l...@gmail.com wrote:
Hi all,
I've searched everywhere to try to find out how to do this and have had no
luck. I am trying to construct identifiers for couples in a dataset.
Essentially, I want to identify couples using more than one column as
identifiers. Take for instance:
obs unit home z sex age
1 015029 18 1 1 053
2 015029 18 1 2 049
3 015029 01 1 1 038
4 015029 01 1 2 033
5 015029 02 1 1 036
6 015029 02 1 2 033
7 015029 03 1 1 023
8 015029 03 1 2 019
9 015029 04 1 2 045
10 015029 05 1 2 047
Where unit is the housing unit, home is household. Of course, there are more
values for unit, although these first ten observations consist of the same
unit (which could possibly be an apartment complex). Nonetheless, I want to
construct an identifier for couples if unit, home match, but only if both
male and female are within the same household. Taking the example data
above, I want to see this:
unit home z sex age couple
1 015029 18 1 1 053 1
2 015029 18 1 2 049 1
3 015029 01 1 1 038 2
4 015029 01 1 2 033 2
5 015029 02 1 1 036 3
6 015029 02 1 2 033 3
7 015029 03 1 1 023 4
8 015029 03 1 2 019 4
9 015029 04 1 2 045 0
10 015029 05 1 2 047 0
As you can see in the last two observations, there were no males identified
within the same household, thus the last two observations would not contain
couple identifiers, rather some other identifier (but the same one) so I can
detect them and remove them later. I've tried using the duplicated function
but was not very useful.
Any help would be greatly appreciated!!!
Thanks,
Carlos
--
View this message in context:
http://r.789695.n4.nabble.com/identify-duplicate-from-more-than-one-column-tp4035888p4035888.html
Sent from the R help mailing list archive at Nabble.com.
__
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.
--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, ATS Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/
__
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
Re: [R] Second-order effect in Parametric Survival Analysis
Thank you Dr. David. I try to summarize it. Assumes x and z are two covariates: x = dummy variable (1 or 0) z = factors (people name) x*z = x + z + x*z therefore this is not a 2nd-order interactions, it should be (for an exponential survival regression):- h(t|(X=x,Z=z)) = exp(Beta0 + XZBeta1) #--- I believe there is no 2nd-order interactions survival regression as I searched over www.rseek.org. While I tried to read through the codes of survreg(), I stuck (cannot understand) at survreg6.c survreg6.c apply C Language which involves Cholesky decomposition multi-matrix (first-order interactions) calculation. 1) chinv2.c 2) cholesky3.c 3) chsolve2.c (only solve the equations of first-order interactions) If someone gives some idea or suggestion on these? Thank you. Best, Ryusuke -- View this message in context: http://r.789695.n4.nabble.com/Second-order-effect-in-Parametric-Survival-Analysis-tp4034318p4036005.html Sent from the R help mailing list archive at Nabble.com. __ 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.
