[R] How/When to use :: and :::?
Hi Dear R users, I check help(Syntax) but I still don't know how to use :: and ::: and when I should use it. Can anyone give me a answer? Thanks in advance Chunhao Tu __ 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] Question about nlmer
Hi R user, I have a question about "nlmer". Can I specify the multilevel random effects in the the nlmer? For example, level 1 is "hospital" and level 2 is patient, so nlmer(response~SSfol(Dose,Time,lke,lka,lcl)~(lke+lka+lcl|hospital)~(lke+lka+lcl|patient%in%hospital), data=health, start=c(lke=-2.5,lka=0.5,lcl=-3),verb=1)) I can do it in the nlme but I still don't know how to get it in nlmer! many thanks in advance Chunhao __ 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] ubuntu repository for R
Hi All R+Linux users, I am very new to Linux and I also have the (dis)similar problem. I try to install R 2.7.0 on Ubuntu. But I only got the R 2.6.2 version. Could anyone tell me how to uninstall R 2.6.2 and install R2.7.0 on Ubuntu. many thanks in advance. Chunhao Quoting Graham Smith <[EMAIL PROTECTED]>: Maybe this is more a Ubuntu question,but... I am trying to get synaptic to work with R I have added http://cran.uk.r-project.org/bin/linux/ubuntu as the URL and hardy/ as the distribution. But when I reload the repositories, I get an error Failed to fetch http://cran.uk.r-project.org/bin/linux/ubuntu/hardy/Packages.gz 301 Moved Permanently Some index files failed to download, they have been ignored, or old ones used instead. I'm still struggling with Linux/Ubuntu so can anyone suggest what I might be doing wrong. Many thanks, Graham [[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.
Re: [R] For loop
Hi,
I guess this what your need. I assume your output is 10 by 11 matrix.
However, you still need to define your geneset function(?)
result<- list()
output<-matrix(NA, nrow=10, ncol=11)
for(i in 1:length(ncol(output)))
{
result[[i]] <- geneset(which(geneset %n% output[1,]))
}
Chunhao
Quoting Rajasekaramya <[EMAIL PROTECTED]>:
Hi,
Could you please let me know to use a list in a for loop here geneset is a
loop.I am trying to match the names of the list with 1st row of the output.
result<- list()
for(i in 1:length(output)
{
result[[i]] <- geneset(which(geneset %n% output[,1]))
}
Kindly help me out
--
View this message in context:
http://www.nabble.com/For-loop-tp18163665p18163665.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] NLME questions -- interpretation of results
Hi Jenny, I try your code but I did not get in converge in fm3 (see the below). For the first question, you could use fm1 to interpret the result without bothering fm2 and fm3. It means that R0 and lrc can be treated as pure fixed effects (Pinherir and Bates, 2000 Book). For the second question, your want to know "is AB/E different from the AB/S" The simplest way is to change your fixed statement: fixed = Asym+R0+lrc ~ dir %in% loc and specify the correct length of starting values. If I am wrong please correct me~ Hope this helpful. Chunhao Tu test<-read.table(file="C:\\Documents and Settings\\ado_cabgfaculty\\Desktop\\sun.txt", header=T) LAST<-groupedData(Y~X|loc/dir, data=test) fm1 <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST, + random = Asym ~1, + fixed = Asym+R0+lrc ~ 1, + start=c(Asym = 0.97, R0 = 1.14, lrc = -0.18)) fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1)) fm3 <- update(fm2, random = pdDiag(Asym+R0+lrc~ 1)) Error in nlme.formula(model = Y ~ SSlogis(X, Asym, R0, lrc), data = LAST, : Step halving factor reduced below minimum in PNLS step Quoting Jenny Sun <[EMAIL PROTECTED]>: My special thanks to Chunhao Tu for the suggestions about testing significance of two locations. I used logistic models to describe relationships between Y and X at two locations (A & B). And within each location, I have four groups (N,E,S,W)representing directions. So the test data can be arranged as: Y X dir loc 0.6295 0.8667596S A 0.7890 0.7324820S A 0.4735 0.9688875S A 0.7805 1.1125239S A 0.8640 0.9506174E A 0.9445 0.6582157E A 0.8455 0.5558860E A 0.9380 0.3304870E A 0.4010 1.1763090N A 0.2585 1.3202890N A 0.3750 1.1763090E A 0.3855 1.3202890E A 0.3020 1.1763090S A 0.2300 1.3202890S A 0.3155 1.1763090W A 0.8890 0.6915861W B 0.9185 0.6149019W B 0.9275 0.5289258W B 0.8365 0.9507088S B 0.7720 0.8842165N B 0.8615 0.8245123N B 0.9170 0.7559687W B 0.9590 0.6772720W B 0.9900 0.5872023W B 0.9940 0.4849064W B 0.7500 0.9560776W B The data is grouped using: LAST<-groupedData(Y~X|loc/dir, data=test) I then used logistic models to define the relationship between Y and X, and got fm1, fm2, and fm3 as follows: -- fm1 <- nlme(DIFN ~ SSlogis(SVA, Asym, R0, lrc),data = LAST,fixed = Asym + R0 + lrc ~ 1,random = Asym ~ 1,start =c(Asym = 1, R0 = 1, lrc = -5)) fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1)) fm3 <- update(fm2, random = pdDiag(Asym + R0 + lrc ~ 1)) anova(fm1,fm2,fm3) ANOVA showed: anova(fm1,fm2,fm3) Model df AIC BIC logLik Test L.Ratio p-value fm1 1 7 -1809.913 -1774.304 910.9564 fm2 2 9 -1805.774 -1758.295 910.8871 1 vs 2 0.1386696 0. fm3 3 12 -1801.822 -1742.473 910.9109 2 vs 3 0.0475543 0.9666 ** question: do the results show that fm1 could represent the results of fm2 and fm3? coef(fm1) Asym R0lrc AB/E 0.9148927 1.389432 -0.3009858 AB/N 0.8775250 1.389432 -0.3009858 AB/S 0.9247592 1.389432 -0.3009858 AB/W 0.8479180 1.389432 -0.3009858 BC/E 0.8791908 1.389432 -0.3009858 BC/N 0.8414229 1.389432 -0.3009858 BC/S 0.9169323 1.389432 -0.3009858 BC/W 0.8817838 1.389432 -0.3009858 ** question: how could I know if any of the models is significantly different from the other ones? (eg. AB/E is different from the AB/S)? summary(fm1) Nonlinear mixed-effects model fit by maximum likelihood Model: DIFN ~ SSlogis(SVA, Asym, R0, lrc) Data: LAST AIC BIC logLik -1809.913 -1774.304 910.9564 Random effects: Formula: Asym ~ 1 | loc Asym StdDev: 2.303402e-05 Formula: Asym ~ 1 | dir %in% loc Asym Residual StdDev: 0.03208693 0.1741559 Fixed effects: Asym + R0 + lrc ~ 1 Value Std.Error DF t-value p-value Asym 0.8855531 0.015375906 2783 57.59355 0 R01.3894322 0.009418047 2783 147.52869 0 lrc -0.3009858 0.012833066 2783 -23.45393 0 Correlation: Asym R0 R0 -0.440 lrc -0.452 0.150 Standardized Within-Group Residuals: Min Q1Med Q3Max -4.1326757 -0.6117037 0.1082112 0.6575250 3.3297270 Number of Observations: 2793 Number of Groups: loc dir %in% loc 28 I have marked all the codes and questions(**). Any answers and suggestions are appreciated. Have a good day! Jenny __ 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, mi
Re: [R] NLME questions -- interpretation of results
Hi Jenny, (I use the data you provide in the previous e-mail)
For the 1st question, let me assume you only want to compare loc: A vs. B
So you could specified your code like this:
fmAB <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST,
random = Asym ~1,
fixed = Asym+R0+lrc ~ loc,
start=c(0.97,0,
1.14,0,
-0.18,0))
summary(fmAB)
Nonlinear mixed-effects model fit by maximum likelihood
Model: Y ~ SSlogis(X, Asym, R0, lrc)
Data: LAST
AIC BIC logLik
-31.02303 -19.70017 24.51152
Random effects:
Formula: Asym ~ 1 | loc
Asym.(Intercept)
StdDev: 1.549779e-06
Formula: Asym ~ 1 | dir %in% loc
Asym.(Intercept) Residual
StdDev: 6.124237e-08 0.09426086
Fixed effects: Asym + R0 + lrc ~ loc
Value Std.Error DF t-value p-value
Asym.(Intercept) 0.9477404 0.1037337 14 9.136280 0.
Asym.locB 0.0982456 0.4175971 14 0.235264 0.8174
R0.(Intercept)1.1289387 0.0652189 14 17.310003 0.
R0.locB 0.1390656 0.3946717 14 0.352358 0.7298
lrc.(Intercept) -0.2110057 0.0656513 14 -3.214036 0.0062
lrc.locB -0.0820484 0.6826382 14 -0.120193 0.9060
Then you know Asym, R0, and lrc of loc B are not significant.
Moreover, you can test the joint fixed effect by
anova(fmAB)(Pinherio and Bate, 2000 Book, p 374)
for the 2nd question, How to get the fitted value for particular level?
Based on this example, let me assume you want to get the fitted value of A/N.
then you could write a small code like this:
FV<-data.frame(F.V=fitted(fmAB), group=summary(fmAB)$groups$dir)
A.N<-FV[is.element(FV$group, c("A/N")),]
F.V group
9 0.4209011 A/N
10 0.2726129 A/N
hope this is helpful~
Chunhao
Quoting Jenny Sun <[EMAIL PROTECTED]>:
Thank you for your reply Chunhao!
I attached only part of the test data and that is why you might not
be able to get convergence. Sorry.
I have a couple more questions:
For the second question you answered, how to specify the correct
length of starting values. I tried using the length of levels in
each of the parameters in the start list but found:
fm1 <- nlme(DIFN ~ SSlogis(SVA, Asym, R0, lrc),data = LAST,fixed =
Asym + R0 + lrc ~ dir %in% loc,random = Asym ~ 1,start =list(Asym =
c(1,1,1,1), R0 = c(1,1,1,1), lrc = c(-5,-2,-2,-2)))
Error in nlme.formula(DIFN ~ SSlogis(SVA, Asym, R0, lrc), data = LAST, :
start must have a component called "fixed"
I've got two loc levels (A,B) with four group levels(N,E,S,W); How
I am gonna define the list and the component called"fixed"?
My another question is about the fitted value of the model. If I
want to calculate adjusted R square, I have to get fitted(fm1).
WHich has values like this;
>fitted(fm1)[1:40]
AB/N AB/N AB/E AB/S AB/W AB/W
AB/W AB/W AB/W AB/W
0.6541876 0.7421748 0.8408251 0.5879220 0.4889387 0.6129576
0.5097593 0.6195679 0.5152567 0.5680860
AB/W AB/W AB/W AB/W AB/W AB/W
AB/N AB/N AB/N AB/E
0.4724423 0.8128148 0.7674529 0.7106698 0.6553155 0.6074771
0.5036201 0.5464105 0.6062978 0.6878438
AB/N AB/N AB/N AB/S AB/S AB/S
AB/S AB/S AB/S AB/S
0.7792725 0.8411961 0.7942503 0.7354845 0.5895700 0.6781973
0.6286886 0.5212052 0.8864748 0.8370021
AB/S AB/S AB/N AB/N AB/N AB/N
AB/E AB/E AB/E AB/E
0.7750731 0.7147024 0.6625288 0.5492599 0.5959280 0.6612426
0.7501786 0.8498928 0.6274681 0.7118615
My question is how to get the fitted values for specified group
levels (eg. values for AB/E)?
Thank all very much!
Jenny
Hi Jenny,
I try your code but I did not get in converge in fm3 (see the below).
For the first question, you could use fm1 to interpret the result
without bothering fm2 and fm3. It means that R0 and lrc can be treated
as pure fixed effects (Pinherir and Bates, 2000 Book).
For the second question, your want to know "is AB/E different from the AB/S"
The simplest way is to change your fixed statement:
fixed = Asym+R0+lrc ~ dir %in% loc
and specify the correct length of starting values.
If I am wrong please correct me~
Hope this helpful.
Chunhao Tu
test<-read.table(file="C:\\Documents and
Settings\\ado_cabgfaculty\\Desktop\\sun.txt", header=T)
LAST<-groupedData(Y~X|loc/dir, data=test)
fm1 <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST,
+ random = Asym ~1,
+ fixed = Asym+R0+lrc ~ 1,
+ start=c(Asym = 0.97, R0 = 1.14, lrc = -0.18))
fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1))
fm3 <- update(fm2, random = pdDiag(Asym+R0+lrc~ 1))
Error in nlme.formula(model = Y ~ SSlogis(X, Asym, R0, lrc), data = LAST, :
Step halving factor reduced below minimum in PNLS step
Quoting Jenny Sun <[EMAIL PROTECTED]>:
My special thanks to Chunhao Tu fo
Re: [R] Matrix set value problem
Hi, when you do the trunc the mx is not a real integer 1 so you must round up m<-matrix(data=NA, nrow=10,ncol=10) i<-1001 mx<-round(trunc(i/1000)) my<-round((i/1000-mx)*1000) m[mx,my]<-1 m [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,]1 NA NA NA NA NA NA NA NANA [2,] NA NA NA NA NA NA NA NA NANA [3,] NA NA NA NA NA NA NA NA NANA [4,] NA NA NA NA NA NA NA NA NANA [5,] NA NA NA NA NA NA NA NA NANA [6,] NA NA NA NA NA NA NA NA NANA [7,] NA NA NA NA NA NA NA NA NANA [8,] NA NA NA NA NA NA NA NA NANA [9,] NA NA NA NA NA NA NA NA NANA [10,] NA NA NA NA NA NA NA NA NANA Chunhao Tu Quoting Valentino Botta <[EMAIL PROTECTED]>: the following code: m<-matrix(nrow=10,ncol=10) i<-1001 mx<-trunc(i/1000) my<-(i/1000-mx)*1000 m[mx,my]<-1 does not assign the value at the matrix m[1,1] Any hints? Thanks v -- View this message in context: http://www.nabble.com/Matrix-set-value-problem-tp18253809p18253809.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] looking for alternative of 'if-else'
Why don't you try "switch" Let me assume that you want to calculate
a=3 and b =5 by using one of "MEAN", "SUM","MIN", and "MAX functions.
Then you could write your code: (If you want to calculate "MEAN")
example<-function(fun=c("MEAN", "SUM", "MIN", "MAX")){
+ fun<-match.arg(fun);a=3;b=5
+ switch( fun,
+MEAN=(a+b)/2,
+SUM=a+b,
+MIN=min(a,b),
+MAX=max(a,b))
+}
example("MEAN")
[1] 4
Hope this is helpful,
Chunhao Tu
Quoting Arun Kumar Saha <[EMAIL PROTECTED]>:
There is "if-else" loop if I have to choose 1 item from a 2-item list.
However if I have a list of 4 items (let say) then how i can choose a single
item without employing 'if-else' loop? I mean in VBA I can use
"select-case", is there any equivalent in R as well?
Regards,
[[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.
Re: [R] looking for alternative of 'if-else'
Hi,
I write the code for you. Hope this helps
SD <- 1
date1 <- seq(as.Date("2001-01-01"), as.Date("2002-12-1"), by = "day")
len1 <- length(date1)
set.seed(1) # to make it reproducible
data1 <- zoo(rnorm(len1), date1)
plot(data1)
# calculation for monthly or quarterly mean
calc = function(frequ = c("M.S", "Q.S")){
+frequ<-match.arg(frequ)
+switch(frequ,
+ M.S = aggregate(data1, as.yearmon, mean),
+ Q.S = aggregate(data1, as.yearqtr, mean))
+ }
Take a look what is the difference between my code and yours then you
will know why your code does not work.
Chunhao Tu
Quoting Arun Kumar Saha <[EMAIL PROTECTED]>:
I am not sure whether I could understand all things. Here is my code :
library(zoo)
# an reproducible example
SD <- 1
date1 <- seq(as.Date("2001-01-01"), as.Date("2002-12-1"), by = "day")
len1 <- length(date1)
set.seed(1) # to make it reproducible
data1 <- zoo(rnorm(len1), date1)
plot(data1)
# calculation for monthly or quarterly mean
calc = function(frequ = c("M S", "Q S"))
{
switch(frequ
"M S" = aggregate(data1, as.yearmon, mean)
"Q S" = aggregate(data1, as.yearqtr, mean))
}
However I am getting lot of errors while executing. Any further suggestion?
Regards,
On Sun, Jul 6, 2008 at 3:38 PM, <[EMAIL PROTECTED]> wrote:
Why don't you try "switch" Let me assume that you want to calculate a=3 and
b =5 by using one of "MEAN", "SUM","MIN", and "MAX functions. Then you could
write your code: (If you want to calculate "MEAN")
example<-function(fun=c("MEAN", "SUM", "MIN", "MAX")){
+ fun<-match.arg(fun);a=3;b=5
+ switch( fun,
+MEAN=(a+b)/2,
+SUM=a+b,
+MIN=min(a,b),
+MAX=max(a,b))
+}
example("MEAN")
[1] 4
Hope this is helpful,
Chunhao Tu
Quoting Arun Kumar Saha <[EMAIL PROTECTED]>:
There is "if-else" loop if I have to choose 1 item from a 2-item list.
However if I have a list of 4 items (let say) then how i can choose a
single
item without employing 'if-else' loop? I mean in VBA I can use
"select-case", is there any equivalent in R as well?
Regards,
[[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.
Re: [R] Method for checking automatically which distribtions fits a data
Hi, In my experience, I just plot the data set then figure it out. Maybe you could try this? I really wonder that there is such a R function exists. If yes, please let me know. Thanks Chunhao Tu Quoting Gundala Viswanath <[EMAIL PROTECTED]>: Hi, Suppose I have a vector of data. Is there a method in R to help us automatically suggest which distributions fits to that data (e.g. normal, gamma, multinomial etc) ? - Gundala Viswanath Jakarta - Indonesia __ 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] GLM, LMER, GEE interpretation
I agree with Ben. Theoretically, Laplace (lmer) provides a better
approximation.
Chunhao
Quoting Ben Bolker <[EMAIL PROTECTED]>:
Daniel Malter umd.edu> writes:
Hi, my dependent variable is a proportion ("prob.bind"), and the independent
variables are factors for group membership ("group") and a covariate
("capacity"). I am interested in the effects of group, capacity, and their
interaction. Each subject is observed on all (4) levels of capacity (I use
capacity as a covariate because the effect of this variable is normatively
linear). I fit three models, but I am observing differences between the
three.
The first model is a quasibinomial without any subject effects using glm.
The second is a random-effects model using lmer. The third model is a
generalized estimating equation using gee from the gee package in which I
cluster for the subject using an unstructured correlation matrix. The
results of the first and the third model almost coincide, but the second,
using lmer, shows an insginficant coefficient where I would expect a
significant one. The other 2 models show the coefficient significant. I do
not really have experience with gee. Therefore I apologize in advance for my
ignorant question whether one of lmer and gee is preferable over the other
in this setting?
[glm]
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.4274 0.4641 -7.386 1.10e-12 ***
capacity 0.9931 0.1281 7.754 9.55e-14 ***
group20.7242 0.6337 1.143 0.25392
group32.0264 0.6168 3.286 0.00112 **
capacity:group2 -0.1523 0.1764 -0.863 0.38864
capacity:group3 -0.3885 0.1742 -2.231 0.02633 *
[lmer]
Generalized linear mixed model fit using Laplace
Formula: prob.bind ~ capacity * group + (1 | subject)
Subset: c(combination == "gnl")
Family: quasibinomial(logit link)
[snip]
Fixed effects:
Estimate Std. Error t value
(Intercept) -3.8628 1.2701 -3.041
capacity 1.1219 0.1176 9.542
group20.9086 1.7905 0.507
group32.3700 1.7936 1.321
capacity:group2 -0.1745 0.1610 -1.083
capacity:group3 -0.3807 0.1622 -2.348
[gee]
Coefficients:
Estimate Naive S.E.Naive z Robust S.E. Robust z
(Intercept) -3.4798395 0.4910274 -7.0868545 0.4739913 -7.3415687
capacity 1.0149659 0.1366365 7.4282170 0.1284162 7.9037210
group2 0.7781014 0.6691731 1.1627806 0.7424769 1.0479807
group3 2.0720270 0.6527565 3.1742727 0.6234005 3.3237495
capacity:group2 -0.1750448 0.1877361 -0.9323982 0.2060484 -0.8495325
capacity:group3 -0.4021872 0.1865916 -2.1554413 0.1724780 -2.3318168
I assume you're talking about the differences in
the estimated standard errors of the group3 (and group2)
parameters (everything else looks pretty similar)?
All else being equal I would trust lmer slightly more
than gee (and the non-clustered glm is not reliable for
inference in this situation, since it ignores the clustering) --
but I'm pretty ignorant of gee, so take that with a grain of salt.
I would make the following suggestions --
1. consider whether it even makes sense to test the
significance of the group3 main effect in the presence
of the capacity:group3 interaction. Is the value capacity=0
somehow intrinsically interesting?
2. all of these standard error estimates are pretty crude/
rely on large-sample assumptions (how big is your data set?);
unfortunately more sophisticated estimates of uncertainty
are currently unavailable for GLMMs in lmer. I would try
your problem again with glmmML, just to check that it gives
similar answers to lmer.
3. if you need more advice, consider asking this on r-sig-mixed
instead ...
Ben Bolker
__
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] extracting index list when using tapply()
Hi, How about using "subset"? x1<-tapply(subset(years, length(area)>20), function(x) length(unique(x))) I hope this works Chunhao Quoting hesicaia <[EMAIL PROTECTED]>: Hello, The quick version of my question is how can I extract a matrix instead of a vector using tapply()? I would like to be able to access both the results of tapply() and also the index variables. In case further explanation would help: I am analyzing a large (3million rows x 9 columns) spatial/temporal dataset and am attempting to calculate the number of unique years containing any data within each geographic area (10 degree cells in this case). I can do this, but I also want to extract a subset vector of the index variable (area). My script to calculate the number of unique years containing any data for each area is: x<-tapply(years, area, function(x) length(unique(x))) Now, I want to extract the vector of areas where the number of unique years containing any data is >20, but tapply() only returns a vector of unique years and I was a matrix. I could use a looping function to do this, but tapply() is much faster with large datasets and so I would like to use it if possible. Any help is appreciated. Thanks. -- View this message in context: http://www.nabble.com/extracting-index-list-when-using-tapply%28%29-tp18345794p18345794.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.
[R] Fwd: Re: extracting index list when using tapply()
The following message is provided by Erik Please provide the reproducible code to do this. Generate a sample data set using the random data generating functions and show us what you'd like, we can then more easily help. [EMAIL PROTECTED] wrote: Hi, How about using "subset"? x1<-tapply(subset(years, length(area)>20), function(x) length(unique(x))) I hope this works Chunhao Quoting hesicaia <[EMAIL PROTECTED]>: Hello, The quick version of my question is how can I extract a matrix instead of a vector using tapply()? I would like to be able to access both the results of tapply() and also the index variables. In case further explanation would help: I am analyzing a large (3million rows x 9 columns) spatial/temporal dataset and am attempting to calculate the number of unique years containing any data within each geographic area (10 degree cells in this case). I can do this, but I also want to extract a subset vector of the index variable (area). My script to calculate the number of unique years containing any data for each area is: x<-tapply(years, area, function(x) length(unique(x))) Now, I want to extract the vector of areas where the number of unique years containing any data is >20, but tapply() only returns a vector of unique years and I was a matrix. I could use a looping function to do this, but tapply() is much faster with large datasets and so I would like to use it if possible. Any help is appreciated. Thanks. -- View this message in context: http://www.nabble.com/extracting-index-list-when-using-tapply%28%29-tp18345794p18345794.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. - End forwarded message - __ 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] Loglikelihood for x factorial?
Hi Rers, I have a silly question. I don't know how to express the loglikelihood function of 1/(x!) where x=x1,x2,xn in R. Could anyone give me a hint? Thank you in advance. Chunhao Tu __ 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] What's the T-Value in fisher.test
Hi Dave, As I know there is no T value for Fisher's exact test (no matter you use SAS, R or other packages) Fisher's exact test is for the categorical data analysis. Fisher's exact test for testing the null of independence of rows and columns in a contingency table with fixed marginals so your result is fail to reject the independence. (If I am wrong please correct me) Hope this helps. Chunhao Tu Quoting DaveFrisch <[EMAIL PROTECTED]>: I do not understand how to interpret this to find the T Value for the data. Is there a way to figure this out, or another function that will provide this for me using Fisher's Exact Test? The outcome of my data is listed below. data: DATA p-value = 0.1698 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.6026805 79.8309210 sample estimates: odds ratio 5.430473 -- View this message in context: http://www.nabble.com/What%27s-the-T-Value-in-fisher.test-tp18386075p18386075.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] rounding
I got round(0.55,1) [1] 0.6 round(2.55,1) [1] 2.5 Your answers are also correct. Please check ?round Quoting "Korn, Ed (NIH/NCI) [E]" <[EMAIL PROTECTED]>: Hi, Round(0.55,1)=0.5 Round(2.55,1)=2.6 Can this be right? Thanks, Ed [[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.
Re: [R] How to output 0.405852702657738 as 0.405853 ?
round(0.405852702657738,6)
[1] 0.405853
Super
Chunhao Tu
Quoting Daren Tan <[EMAIL PROTECTED]>:
I am confused by options("digits") and options("scipen"), which
should be used to output 0.405852702657738 as 0.405853 in
write.table ?
_
[[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.
Re: [R] create multi windows plot at the same time
Hi, you only need to add x11() in front of plot(); such as x11() Plot(v1) x11() Plot(v2) .. then R will open multi windows for you Chunhao Quoting Alessandro <[EMAIL PROTECTED]>: Ciao All I'd like to generate 4 windows plot to compare different relationship at the same time in R. Each plot code is : Plot(v1) Plot(v2) Plot(v4) Plot(v4) But I did not find a right code in R. Ale [[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.
Re: [R] R: create multi windows plot at the same time
par(mfrow=c(2,2)) is trying to open a 2*2 figures in one window, so you could see your v1,..v4 in one windows at the same time if you want they plot separately then you could use x11() for each plot() Chunhao Quoting Alessandro <[EMAIL PROTECTED]>: Hi I have a problem with par(mfrow=c(2,2)) code, because I dont understand this: par(mfrow=c(2,2)) plot(v1) plot(v2) plot(v3) plot(v4) it opens a v1 window and after close that window to open next v2 plot Da: stephen sefick [mailto:[EMAIL PROTECTED] Inviato: venerdì 25 luglio 2008 12.13 A: [EMAIL PROTECTED] Cc: Alessandro; r-help@r-project.org Oggetto: Re: [R] create multi windows plot at the same time or quartz() or windows() you could also use par(mfrow=c(2,2)) On Fri, Jul 25, 2008 at 3:05 PM, <[EMAIL PROTECTED]> wrote: Hi, you only need to add x11() in front of plot(); such as x11() Plot(v1) x11() Plot(v2) .. then R will open multi windows for you Chunhao Quoting Alessandro <[EMAIL PROTECTED]>: Ciao All I'd like to generate 4 windows plot to compare different relationship at the same time in R. Each plot code is : Plot(v1) Plot(v2) Plot(v4) Plot(v4) But I did not find a right code in R. Ale [[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. -- Let's not spend our time and resources thinking about things that are so little or so large that all they really do for us is puff us up and make us feel like gods. We are mammals, and have not exhausted the annoying little problems of being mammals. -K. Mullis __ 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] About clustering techniques
Hi Paco, I got the same problem with you before. Thus, I just impute the missing values For example: newdata<-as.matrix(impute(olddata, fun="random")) then I believe that you could analyze your data. Hopefully it helps. Chunhao Quoting pacomet <[EMAIL PROTECTED]>: Hello R users It's some time I am playing with a dataset to do some cluster analysis. The data set consists of 14 columns being geographical coordinates and monthly temperatures in annual files latitutde - longitude - temperature 1 -. - temperature 12 I have some missing values in some cases, maybe there are 8 monthly valid values at some points with four non valid. I don't want to supress the whole row with 8 good/4 bad values as I wanna try annual and monthy analysis. I first tried kmeans but found a problem with missing values. When trying without omitting missing values kmeans gives an error and when excluding invalid data too many values are excluded in some years of the data series. Now I have been reading about pam, pamk and clara, I think they can handle missing values. But can't find out the way to perform the analysis with these functions. As I'm not an statistics nor an R expert the fpc or cluster package documentation is not enough for me. If you know about a website or a tutorial explaining the way to use that functions, with examples to check if possible, please post them. Any other help or suggestion is greatly appreciated. Thanks in advance Paco -- _ El ponent la mou, el llevant la plou Usuari Linux registrat: 363952 --- Fotos: http://picasaweb.google.es/pacomet [[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.
[R] Repeated Measure ANOVA-Old question
Hi R users,
I google the website and I found that there are three ways to perform
repeated measure ANOVA: aov, lme and lmer.
http://www.mail-archive.com/[EMAIL PROTECTED]/msg58502.html
But the questions are which one is good to use and how to do post-hoc test?
I use the example that is provided in the above link and I try
tt<-aov(p.pa~group*time+Error(subject/time),data=P.PA)
TukeyHSD(tt)
Error in UseMethod("TukeyHSD") : no applicable method for "TukeyHSD"
Any suggestions?
Appreciate in advance~
Chunhao
__
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] Re peated Measure ANOVA-Old question
Tank you (Anna & Mark) very much for this super information
Let me confirm that IF I want to perform the RM-ANOVA
I should use "lmer" and perform the post-hoc test by using "glht", right?
Because I am not so familiar with "lmer" so I have two more questions.
Let me use that following example again,
1. Is the random effect statement correct? if no what will be?
2. What did I do wrong in glht?
require(lme4);require(nlme)
group <-rep(rep(1:2, c(5,5)), 3);time <-rep(1:3, rep(10,3));subject <-
rep(1:10,3)
p.pa <- c(92, 44, 49, 52, 41, 34, 32, 65, 47, 58, 94, 82, 48, 60, 47,
46, 41, 73, 60, 69, 95, 53, 44, 66, 62, 46, 53, 73, 84, 79)
P.PA <- data.frame(subject, group, time, p.pa)
P.PA$group=factor(P.PA$group);P.PA$time=factor(P.PA$time)
P.PA$subject=factor(P.PA$subject)
tlmer<-lmer (p.pa ~ time * group + (time*group | subject), data=P.PA )
summary(glht(tlmer,linfct=mcp(group="Tukey")))
Error in UseMethod("fixef") : no applicable method for "fixef"
In addition: Warning message
Thank you very much
Chunhao
__
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] Mixed effects model where nested factor is not the repeated across treatments lme???
Hi Miki,
I just got the same problem with you couple hours ago.
Rusers (Anna, and Mark {thank you guys}) provide me a vary valuable
information.
link to following address.
http://www.nabble.com/Tukey-HSD-(or-other-post-hoc-tests)-following-repeated-measures-ANOVA-td17508294.html#a17559307
for the A vs. B, A vs. C
You could install and download the multcomp package and perform the
post hoc test
such as
summary(glht(lmel,linfct=mcp(treatment="Tukey")))
hopefully it helps
Chunhao
Quoting M Ensbey <[EMAIL PROTECTED]>:
Hi,
I have searched the archives and can't quite confirm the answer to this.
I appreciate your time...
I have 4 treatments (fixed) and I would like to know if there is a
significant difference in metal volume (metal) between the treatments.
The experiment has 5 blocks (random) in each treatment and no block is
repeated across treatments. Within each plot there are varying numbers
of replicates (random) (some plots have 4 individuals in them some have
14 and a range in between). NOTE the plots in one treatment are not
replicated in the others.
So I end up with a data.frame with 4 treatments repeated down one column
(treatment=A, B, C, D), 20 plots repeated down the next (block= 1 to 20)
and records for metal volume (metal- 124 of these)
I have made treatment and block a factor. But haven't grouped them (do I
need to and how if so)
The main question is in 3 parts:
1. is this the correct formula to use for this situation:
lme1<-lme(metal~treatment,data=data,random=~1|block) (or is lme even the
right thing to use here?)
I get:
summary(lme1)
Linear mixed-effects model fit by REML
Data: data
AIC BIClogLik
365.8327 382.5576 -176.9163
Random effects:
Formula: ~1 | block
(Intercept) Residual
StdDev: 0.4306096 0.9450976
Fixed effects: Cu ~ Treatment
Value Std.Error DF t-value p-value
(Intercept) 5.587839 0.2632831 104 21.223688 0. ***
TreatmentB -0.970384 0.3729675 16 -2.601792 0.0193 ***
TreatmentC -1.449250 0.3656351 16 -3.963651 0.0011 ***
TreatmentD -1.319564 0.3633837 16 -3.631323 0.0022 ***
Correlation:
(Intr) TrtmAN TrtmCH
TreatmentB -0.706
TreatmentC -0.720 0.508
TreatmentD -0.725 0.511 0.522
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.85762206 -0.68568460 -0.09004478 0.56237152 3.20650288
Number of Observations: 124
Number of Groups: 20
2. if so how can I get p values for comparisons between every
group... ie is A different from B, is A different from C, is A different
from D, is B different from C, is B different from D etc... is there a
way to get all of these instead of just "is A different from B, is A
different from C, is A different from D" which summary seems to give?
3. last of all what is the best way to print out all the residuals
for lme... I can get qqplot(lme1) is there a pre-programmed call for
multiple diagnostic plots like in some other functions...
Thankyou so Much for your time
It is much appreciated
;-)
Miki
[[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.
Re: [R] How can I create a vector of vector?
Hi Jim, x<-as.vector(list(c(31,2,3), c(1,2,5,34,5656),c(211,3243,5,343,3))) x[1] [[1]] [1] 31 2 3 x[2] [[1]] [1]125 34 5656 x[3] [[1]] [1] 211 32435 3433 Hope this helps Chunhao Quoting Jim <[EMAIL PROTECTED]>: Dear All, How can I create a vector of vector in R? for example, I would input the data as follows: x[1] <- c(31,2,3) x[2] <- c(1,2,5,34,5656) x[3] <- c(211,3243,5,343,3) Jim Liu __ 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.
[R] correlation with boot
Hi R users,
I have one simple question but I really don't know why I can't get it work.
The data was adopted from Efron's An introduction to the bootstrap book.
nlaw
LSAT GPA
[1,] 576 3.39
[2,] 635 3.30
[3,] 558 2.81
[4,] 578 3.03
[5,] 666 3.44
[6,] 580 3.07
[7,] 555 3.00
[8,] 661 3.43
[9,] 651 3.36
[10,] 605 3.13
[11,] 653 3.12
[12,] 575 2.74
[13,] 545 2.76
[14,] 572 2.88
[15,] 594 2.96
nlaw<-as.matrix(law[,-1])
corr.t<-function(data){
+return(boot::corr(data))
+}
law.boot<-boot(data=nlaw,statistic=corr.t,R=49)
Error in statistic(data, original, ...) : unused argument(s) (1:15)
many thanks for the help
Chunhao
__
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] correlation with boot
Hi Professor Ripley,
(Thank you very much)^infinite for your response. I still have few questions
1. R shows that I have two column in nlaw data.
ncol(nlaw)
[1] 2
corr(nlaw)
[1] 0.7763745
2. Why can't I use boot::corr? I check corr in boot package and it
seems to be used for computing correlation coefficient,right?
3. The last question, cor is for computed the correlation variance and
covariance matrix. Would you give a explanation for
cor(data[ind,])[1,2]
I do not get it.
Appreciate,
Chunhao
Quoting Prof Brian Ripley <[EMAIL PROTECTED]>:
Please look at the help for boot, in particular what it says the form
of 'statistic' should be. Possibly what you want is
corr.t <- function(data, ind) cor(data[ind,])[1,2]
nlaw appears to be a single-column matrix: that will not work either.
On Mon, 1 Sep 2008, [EMAIL PROTECTED] wrote:
Hi R users,
I have one simple question but I really don't know why I can't get it work.
The data was adopted from Efron's An introduction to the bootstrap book.
nlaw
LSAT GPA
[1,] 576 3.39
[2,] 635 3.30
[3,] 558 2.81
[4,] 578 3.03
[5,] 666 3.44
[6,] 580 3.07
[7,] 555 3.00
[8,] 661 3.43
[9,] 651 3.36
[10,] 605 3.13
[11,] 653 3.12
[12,] 575 2.74
[13,] 545 2.76
[14,] 572 2.88
[15,] 594 2.96
nlaw<-as.matrix(law[,-1])
corr.t<-function(data){
+return(boot::corr(data))
+}
law.boot<-boot(data=nlaw,statistic=corr.t,R=49)
Error in statistic(data, original, ...) : unused argument(s) (1:15)
many thanks for the help
Chunhao
__
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.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
__
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 for par(mfrow)
Hi, R users
I have a question with par(mfrow). I try to histograms and qqplots
form boot output for 5 statistics but par(mfrow=c(5,2)) or
par(mfrow=c(5,1)) does not work. R still display each figure
separately. What did I do wrong? (I check ?par)
c1.boot<-boot(c1data,c1.fun,R=999)
c1.boot
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = c1data, statistic = c1.fun, R = 999)
Bootstrap Statistics :
originalbiasstd. error
t1* 0.30447696 0.0014228306 0.01026153
t2* 0.05967183 -0.0014505285 0.01274977
t3* 1.11852318 -0.0009339349 0.31586852
t4* 4.12856813 -0.0310221125 0.39602210
t5* -1.29958196 0.0222876889 0.67625299
par(mfrow=c(5,2))
for (i in 1:5) {plot(c1.boot, index=i, col=i)}
Many Thanks
Chunhao
__
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 use try function with boot
Hi R users, Is is possible for me to use the try function with boot? I would to do the bootstraping with a nonlinear model(it works well when R < 1000). But it does not work very well (when R is large) thus I try to use "try" to resolve. I put the try function in two cases: case1: put the try in front of the boot c1.try<-try(boot(c1data, statistic = c1.fun, R=3999),silent=T) c1.try [1] "Error in nls(formula = density ~ nmf(time, alpha, delta, psi, tau, gamma), : \n Convergence failure: false convergence (8)\n" attr(,"class") [1] "try-error" case2: put the try in front of the nls c1.nmf<-try(nls(density~nmf(time, alpha, delta, psi, tau, gamma), +algorithm="port",data=c1, +lower=c(alpha=0.1, delta=0, psi=0.2, tau=1, gamma=-5), +upper=c(alpha=0.6, delta=0.1, psi=3, tau=7, gamma=2), +start=c(alpha=0.35, delta=0, psi=0.99, tau=4.5, gamma=-1.2)),silent=T) c1.try<-boot(c1data, statistic = c1.fun, R=3999) Error in nls(formula = density ~ nmf(time, alpha, delta, psi, tau, gamma), : Convergence failure: iteration limit reached without convergence Any suggestion will be helpful. Many thanks in advance Chunhao __ 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] Help use try function with boot
Hi Jinsong,
I try to put the "try" in the c1.fun but it still does not work. Do
you mean this?
c1.fun<-try(function(data,i){
+ d<-data
+ d$density<-d$fitted+d$res[i]
+ coef(update(c1.nmf,data=d))
+ } ,silent=T)
c1.try<-boot(c1data, statistic = c1.fun, R=5000)
Error in nls(formula = density ~ nmf(time, alpha, delta, psi, tau, gamma), :
Convergence failure: false convergence (8)
Thanks
Chunhao
Quoting Jinsong Zhao <[EMAIL PROTECTED]>:
[EMAIL PROTECTED] wrote:
Hi R users,
Is is possible for me to use the try function with boot? I would to
do the bootstraping with a nonlinear model(it works well when R <
1000). But it does not work very well (when R is large) thus I try
to use "try" to resolve. I put the try function in two cases:
case1: put the try in front of the boot
c1.try<-try(boot(c1data, statistic = c1.fun, R=3999),silent=T)
c1.try
[1] "Error in nls(formula = density ~ nmf(time, alpha, delta, psi,
tau, gamma), : \n Convergence failure: false convergence (8)\n"
attr(,"class")
[1] "try-error"
case2: put the try in front of the nls
c1.nmf<-try(nls(density~nmf(time, alpha, delta, psi, tau, gamma),
+algorithm="port",data=c1,
+lower=c(alpha=0.1, delta=0, psi=0.2, tau=1, gamma=-5),
+upper=c(alpha=0.6, delta=0.1, psi=3, tau=7, gamma=2),
+start=c(alpha=0.35, delta=0, psi=0.99, tau=4.5,
gamma=-1.2)),silent=T)
c1.try<-boot(c1data, statistic = c1.fun, R=3999)
Error in nls(formula = density ~ nmf(time, alpha, delta, psi, tau,
gamma), :
Convergence failure: iteration limit reached without convergence
Any suggestion will be helpful.
Many thanks in advance
Chunhao
put try() in your statistic function. in your case, it may be inside
the c1.fun function.
HTH,
Jinsong
__
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 parametric boot
Hi R users
Is there any example for nonlinear parametric boot? I google it but I
can't find it. I am interested in the parameter estimators of a
nonlinear model. But I really don't know how to code it in the
"ran.gen" statement (data set from ?nls)
fm1 <- nls(weight ~ Asym/(1+exp((xmid-Time)/scal)), data = ChickWeight,
+start=c(Asym=337, xmid=16, scal=8))
fm1.fun<-function(data){coef(update(fm1,data=data))}
ran.sim<-function(data,mle){out<-rnorm(n=nrow(data,mle));out}
fm1.boot<-boot(ChickWeight, statistic = fm1.fun, R=99, sim="parametric",
+ran.gen=ran.sim, mle=coef(fm1))
Error in nrow(data, mle) :
unused argument(s) (c(337.605336871528, 16.0688379710354, 8.00747460385483))
Any suggestion would be very helpful.
many thanks in advance
Chunhao
__
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] S.O.S "try" doesnot work in boot?
First thanks for Jinsong's suggestions
I would like to do a bootstrap in a nonlinear model. But it fails to
converge in most of time. (it did converge if I just use nls without
boot). Thus, I use "try" function to resolve my problem. This
following code is from Jinsong's suggestion.
h1a.nls<-nls(density~nmf(time, alpha, delta, psi, tau, gamma),data=h1a,
start=c(alpha=0.3, delta=0.08869, psi=1.26523, tau=3.93919,
gamma=-1.41927))
h1a.data<-data.frame(h1a,res=resid(h1a.nls),fitted=fitted(h1a.nls))
h1a.fun<-function(data,i){
d<-data
d$density<-d$fitted+d$res[i]
try(update(h1a.nls,data=d),silent=T)
if(!inherits(h1a.nls,"try-error")) h1a.coef<-coef(h1a.nls)
else h1a.coef<-NA
h1a.coef
}
h1a.boot<-boot(h1a.data, statistic = h1a.fun, R=1000)
h1a.boot
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = h1a.data, statistic = h1a.fun, R = 1000)
Bootstrap Statistics :
original biasstd. error
t1* 0.27892590 0 0
t2* 0.08869433 0 0
t3* 1.26523275 0 0
t4* 3.93919567 0 0
t5* -1.41926966 0 0
all of the values of each column in h1a.boot$t are the same.
Is anyone know to how I can solve this problem?
Appreciate in advance
Chunhao
__
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] Bug in "is" ?
Hi R users Is there anything wrong in "is" function? (R 2.7.2) I believe that everyone will agree that "7" is an integer, right? but why R shows 7 is not an integer is.integer(7) [1] FALSE is(7,"integer") [1] FALSE is(as.integer(7), "integer") [1] TRUE Thank you very much in advance Chunhao __ 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] Bug in "is" ?
This is really bothering me! In the Dr. Venables and Dr. Ripley's book "S Programming" Page 105 shows that c(is(10,"integer"),is(10.5,"integer")) [1] T F But I try this in R 2.7.2 it shows c(is(10,"integer"),is(10.5,"integer")) [1] FALSE FALSE Does anyone know what is going on here? Appreciate, Chunhao Quoting Yihui Xie <[EMAIL PROTECTED]>: Yes, everyone will agree "7" is an integer, but I don't think computers will agree too :-) R thinks it's a double-precision number, except when you explicitly specify it as an integer (say, as.integer()). class(7) [1] "numeric" is.double(7) [1] TRUE Regards, Yihui -- Yihui Xie <[EMAIL PROTECTED]> Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086 Mobile: +86-15810805877 Homepage: http://www.yihui.name School of Statistics, Room 1037, Mingde Main Building, Renmin University of China, Beijing, 100872, China On Wed, Sep 24, 2008 at 12:40 PM, <[EMAIL PROTECTED]> wrote: Hi R users Is there anything wrong in "is" function? (R 2.7.2) I believe that everyone will agree that "7" is an integer, right? but why R shows 7 is not an integer is.integer(7) [1] FALSE is(7,"integer") [1] FALSE is(as.integer(7), "integer") [1] TRUE Thank you very much in advance Chunhao __ 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] Bug in "is" ?
Thank you for all of you. Intuitively, 7 is an integer for people who
live in this planet. It is just very difficult for me to believe that
R does not think 7 is an integer but 7L is.
is.integer(7) # R 2.7.2
[1] FALSE
Thus, based on Martin's comments, I try it again on the S-PLUS 8.0 and
it shows
is.integer(7) # S-PLUS 8.0
[1] T
Hopefully, someday and someone will fix it therefore, R users don't
need to use as.integer(7) to tell R that 7 is an integer.
Thanks again
Chunhao
Quoting Martin Maechler <[EMAIL PROTECTED]>:
"KJ" == Keith Jewell <[EMAIL PROTECTED]>
on Wed, 24 Sep 2008 09:46:08 +0100 writes:
KJ> "7" is an integer, but it's also a real.
KJ> In R '?is' and '?is.integer' are clear that you're testing
the class(es) of
KJ> objects, not their values.
KJ> I can't comment on the relationship with "S Programming"
I can:
In S, and S-plus upto version 3.4,
numeric constants such as '7' where "double" as they are in R.
Then in S-plus 5.1, they became "integer",
and there were tools so users could change all(!!) their S
scripts to use '7.' instead of '7' in all places where numeric
constants were seen, in order to keep behavior back compatible.
R never made such a step (backwards ;-), and never will,
notably since in R we had introduced the explicit long (= long
integer) constants, using the 'L' suffix,
i.e., 7L is "integer"
7 is "double"
Note however that for both, is.numeric(.) is fulfilled and
class(.) and mode(.) return "numeric".
Only typeof(.), storage.mode(.) or str(.)
(or functions building on these) tell you the difference.
Martin Maechler, ETH Zurich and R core team
[And, yes, if you think further and are wondering:
If we'd design things from scratch, we would only have S4
classes and "double" would be a proper class and
"numeric" would be the class union of {"integer", "double"}
]
KJ> <[EMAIL PROTECTED]> wrote in message
KJ> news:[EMAIL PROTECTED]
>> This is really bothering me! In the Dr. Venables and Dr.
Ripley's book "S
>> Programming" Page 105
>> shows that
>>> c(is(10,"integer"),is(10.5,"integer"))
>> [1] T F
>>
>> But I try this in R 2.7.2 it shows
>>> c(is(10,"integer"),is(10.5,"integer"))
>> [1] FALSE FALSE
>> Does anyone know what is going on here?
>>
>> Appreciate,
>> Chunhao
>>
>> Quoting Yihui Xie <[EMAIL PROTECTED]>:
>>
>>> Yes, everyone will agree "7" is an integer, but I don't think
>>> computers will agree too :-) R thinks it's a double-precision number,
>>> except when you explicitly specify it as an integer (say,
>>> as.integer()).
>>>
class(7)
>>> [1] "numeric"
>>>
is.double(7)
>>> [1] TRUE
>>>
>>> Regards,
>>> Yihui
>>> --
>>> Yihui Xie <[EMAIL PROTECTED]>
>>> Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
>>> Mobile: +86-15810805877
>>> Homepage: http://www.yihui.name
>>> School of Statistics, Room 1037, Mingde Main Building,
>>> Renmin University of China, Beijing, 100872, China
>>>
>>>
>>>
>>> On Wed, Sep 24, 2008 at 12:40 PM, <[EMAIL PROTECTED]> wrote:
Hi R users
Is there anything wrong in "is" function? (R 2.7.2)
I believe that everyone will agree that "7" is an integer,
right? but
why R
shows 7 is not an integer
> is.integer(7)
[1] FALSE
>
> is(7,"integer")
[1] FALSE
>
> is(as.integer(7), "integer")
[1] TRUE
Thank you very much in advance
Chunhao
KJ> __
KJ> R-help@r-project.org mailing list
KJ> https://stat.ethz.ch/mailman/listinfo/r-help
KJ> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
KJ> 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.
__
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] Bug in "is" ?
Hi Keith,
No doubt, 7.0 is integer in math. But if people can write 7 why people
need to write 7.0 (I do not see any reason to do this). My point is
that R maybe can do something like S-plus. No point to argue. don't
you think so?
Thanks
Chunhao
Quoting Keith Jewell <[EMAIL PROTECTED]>:
Have you tried is.integer(7.0) in S-Plus? (I have)
Do you think 7.0 is integer?
IMHO in R there is nothing to be fixed (in this regard) except your
understanding.
This is a computer language, not English; intuition isn't reliable, so we
have help pages.
is.integer(x) is not intended to indicate whether the value of x is a whole
number, it indicates whether x has class "integer".
All objects of class "integer" have whole number values, but not all objects
with whole number values have class "integer".
If you want to know whether a value is a whole number you could try (but
there may be a better way, and beware of computer precision)
x == as.integer(x)
If you want a value to be stored in an object of class "integer" you'd
better say so (using as.integer or L or ...), else how is R to know what you
want? As Martin has pointed out, the system could "guess" based on the
presence or absence of a decimal point; I share his opinion that this would
be a "bad thing".
Nuff said.
Keith J
<[EMAIL PROTECTED]> wrote in message
news:[EMAIL PROTECTED]
Thank you for all of you. Intuitively, 7 is an integer for people who
live in this planet. It is just very difficult for me to believe that R
does not think 7 is an integer but 7L is.
is.integer(7) # R 2.7.2
[1] FALSE
Thus, based on Martin's comments, I try it again on the S-PLUS 8.0 and it
shows
is.integer(7) # S-PLUS 8.0
[1] T
Hopefully, someday and someone will fix it therefore, R users don't need
to use as.integer(7) to tell R that 7 is an integer.
Thanks again
Chunhao
Quoting Martin Maechler <[EMAIL PROTECTED]>:
"KJ" == Keith Jewell <[EMAIL PROTECTED]>
on Wed, 24 Sep 2008 09:46:08 +0100 writes:
KJ> "7" is an integer, but it's also a real.
KJ> In R '?is' and '?is.integer' are clear that you're testing the
class(es) of
KJ> objects, not their values.
KJ> I can't comment on the relationship with "S Programming"
I can:
In S, and S-plus upto version 3.4,
numeric constants such as '7' where "double" as they are in R.
Then in S-plus 5.1, they became "integer",
and there were tools so users could change all(!!) their S
scripts to use '7.' instead of '7' in all places where numeric
constants were seen, in order to keep behavior back compatible.
R never made such a step (backwards ;-), and never will,
notably since in R we had introduced the explicit long (= long
integer) constants, using the 'L' suffix,
i.e., 7L is "integer"
7 is "double"
Note however that for both, is.numeric(.) is fulfilled and
class(.) and mode(.) return "numeric".
Only typeof(.), storage.mode(.) or str(.)
(or functions building on these) tell you the difference.
Martin Maechler, ETH Zurich and R core team
[And, yes, if you think further and are wondering:
If we'd design things from scratch, we would only have S4
classes and "double" would be a proper class and
"numeric" would be the class union of {"integer", "double"}
]
KJ> <[EMAIL PROTECTED]> wrote in message
KJ> news:[EMAIL PROTECTED]
>> This is really bothering me! In the Dr. Venables and Dr.
Ripley's book "S
>> Programming" Page 105
>> shows that
>>> c(is(10,"integer"),is(10.5,"integer"))
>> [1] T F
>>
>> But I try this in R 2.7.2 it shows
>>> c(is(10,"integer"),is(10.5,"integer"))
>> [1] FALSE FALSE
>> Does anyone know what is going on here?
>>
>> Appreciate,
>> Chunhao
>>
>> Quoting Yihui Xie <[EMAIL PROTECTED]>:
>>
>>> Yes, everyone will agree "7" is an integer, but I don't think
>>> computers will agree too :-) R thinks it's a double-precision
number,
>>> except when you explicitly specify it as an integer (say,
>>> as.integer()).
>>>
class(7)
>>> [1] "numeric"
>>>
is.double(7)
>>> [1] TRUE
>>>
>>> Regards,
>>> Yihui
>>> --
>>> Yihui Xie <[EMAIL PROTECTED]>
>>> Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086
>>> Mobile: +86-15810805877
>>> Homepage: http://www.yihui.name
>>> School of Statistics, Room 1037, Mingde Main Building,
>>> Renmin University of China, Beijing, 100872, China
>>>
>>>
>>>
>>> On Wed, Sep 24, 2008 at 12:40 PM, <[EMAIL PROTECTED]> wrote:
Hi R users
Is there anything wrong in "is" function? (R 2.7.2)
I believe that everyone will agree that "7" is an integer,
right? but
why R
shows 7 is not an integer
> is.integer(7)
[1] FALSE
>
> is(7,"integer")
[1] FALSE
>
> is(as.integer(7), "integer")
[1] TRUE
Thank you very much in adva
Re: [R] goodness of fit test
Hi Edna, You could use lapply or sapply to perform the multiple goodness of fit tests at the same time. Chunhao Quoting Edna Bell <[EMAIL PROTECTED]>: Dear R Gurus; Is there an automated process for goodness of fit tests, please? I know there is prop.test for one at a time, but I was wondering about this, please. Thanks, Edna Bell __ 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.
[R] How to write a package based on nlme
Dear R Helpers,
I try to write a small package that based on nlme however my code does
not work.
R always appears this message:
Error in eval(expr, envir, enclos) : object "y" not found
where y is the response variable. Please help me out!
This is my code:
require(nlme)
AMPmixed<-function(y, x, S1=c("asymptotic","logistic"), random,data,
S2=c("asymptotic","logistic"), start)
{
logist.logist<-function(x,alpha,delta,psi.l,tau.l,gamma,h){
delta+(alpha-delta+gamma*(x-(h-1))/exp(x))/(1+exp(-(x-tau.l)/psi.l))}
logist.asymp<-function(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h){
delta+(alpha-delta)/(1+exp(-(x-tau.l)/psi.l))+(gamma*(x-(h-1))/exp(x))*exp(-exp(1/lpsi.a)*x)}
asymp.asymp<-function(x,alpha,delta,lpsi.a,gamma,h){
delta+(alpha-delta)*exp(-exp(1/lpsi.a)*x)+(gamma*(x-(h-1))/exp(x))*exp(-exp(1/lpsi.a)*x)}
asymp.logist<-function(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h){
delta+(alpha-delta)*exp(-exp(1/lpsi.a)*x)+(gamma*(x-(h-1))/exp(x))/(1+exp(-(x-tau.l)/psi.l))}
(logistic.logistic<-function(y, x, data, start, random){
nlme.out<-nlme(y~logist.logist(x,alpha,delta,psi.l,tau.l,gamma,h),
data=data, start=start, na.action)
list(nlme.out=summary(nlme.out))
})
(logistic.asymptotic<-function(y, x, data, start, random){
nlme.out<-nlme(y~logist.asymp(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h),
data=data, start=start, na.action)
list(nlme.out=summary(nlme.out))
})
(asymptotic.logistic<-function(y, x, data, start,random){
xy<-sortedXyData(x,y,data=data)
if(nrow(xy)<8){stop("Too few factor levels to fit an
asymptotic.logistic")}
nlme.out<-nlme(y~asymp.logist(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h),
data=data, start=start, na.action)
list(nlme.out=summary(nlme.out))
})
(asymptotic.asymptotic<-function(y, x, data, start, random){
nlme.out<-nlme(y~asymp.asymp(x,alpha,delta,lpsi.a,gamma,h), data=data,
start=start,
fixed=alpha+delta+lpsi.a+gamma+h~1,random=random)
list(nlme.out=summary(nlme.out))
})
if(S1=="logistic" && S2=="logistic")
{(AMPmixed=logistic.logistic(y, x, data, start, random))}
else if(S1=="logistic" &&
S2=="asymptotic"){(AMPmixed=logistic.asymptotic(y, x, data, start,
random))}
else if(S1=="asymptotic" &&
S2=="logistic"){(AMPmixed=asymptotic.logistic(y, x, data, start,
random))}
else if(S1=="asymptotic" &&
S2=="asymptotic"){(AMPmixed=asymptotic.asymptotic(y, x, data, start,
random))}
}
con rep biomass
10.00 1 1.126
20.32 1 1.096
31.00 1 1.163
43.20 1 0.985
5 10.00 1 0.716
6 32.00 1 0.560
7 100.00 1 0.375
80.00 2 0.833
90.32 2 1.106
10 1.00 2 1.336
11 3.20 2 0.754
12 10.00 2 0.683
13 32.00 2 0.488
14 100.00 2 0.344
iso<-read.table(file="E:\\Hormesis\\data\\isobutanol.txt", header=T)
aa<-groupedData(biomass~con|rep, data=iso)
van2<-AMPmixed(y=biomass, x=con, S1="asymptotic", S2="asymptotic", data=aa,
+ start=c(alpha= 0.7954, delta= 0.3231, lpsi.a=-0.2738,
gamma= 1.0366, h=0.8429))
Error in eval(expr, envir, enclos) : object "y" not found
Thank you very much~~
Chunhao
__
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] Need help for nlme
Hi everyone,
I try to write a module based on nlme however R always shows me the
error message
Error in eval(expr, envir, enclos) : object "y" not found. Does anyone
know how to solve this? There is no problem in nls at all.
require(nlme)
AMPmixed<-function(y, x, S1=c("asymptotic","logistic"),
S2=c("asymptotic","logistic"), data, start,random)
{
logist.logist<-function(x,alpha,delta,psi.l,tau.l,gamma,h){
delta+(alpha-delta+gamma*(x-(h-1))/exp(x))/(1+exp(-(x-tau.l)/psi.l))}
logist.asymp<-function(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h){
delta+(alpha-delta)/(1+exp(-(x-tau.l)/psi.l))+(gamma*(x-(h-1))/exp(x))*exp(-exp(1/lpsi.a)*x)}
asymp.asymp<-function(x,alpha,delta,lpsi.a,gamma,h){
delta+(alpha-delta)*exp(-exp(1/lpsi.a)*x)+(gamma*(x-(h-1))/exp(x))*exp(-exp(1/lpsi.a)*x)}
asymp.logist<-function(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h){
delta+(alpha-delta)*exp(-exp(1/lpsi.a)*x)+(gamma*(x-(h-1))/exp(x))/(1+exp(-(x-tau.l)/psi.l))}
(logistic.logistic<-function(y, x, data, start, random){
nlme.out<-nlme(y~logist.logist(x,alpha,delta,psi.l,tau.l,gamma,h),
data=data, start=start,
fixed=alpha+delta+psi.l+tau.l+gamma+h~1,
random=random)
list(nlme.out=summary(nlme.out))
})
(logistic.asymptotic<-function(y, x, data, start, random){
nlme.out<-nlme(y~logist.asymp(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h),
data=data, start=start,
fixed=alpha+delta+psi.l+tau.l+lpsi.a+gamma+h~1, random=random)
list(nlme.out=summary(nlme.out))
})
(asymptotic.logistic<-function(y, x, data, start,random){
nlme.out<-nlme(y~asymp.logist(x,alpha,delta,psi.l,tau.l,lpsi.a,gamma,h),
data=data, start=start,
fixed=alpha+delta+psi.l+tau.l+lpsi.a+gamma+h~1, random=random)
list(nlme.out=summary(nlme.out))
})
(asymptotic.asymptotic<-function(y, x, data, start, random){
nlme.out<-nlme(y~asymp.asymp(x,alpha,delta,lpsi.a,gamma,h), data=data,
start=start,
fixed=alpha+delta+lpsi.a+gamma+h~1,random=random)
list(nlme.out=summary(nlme.out))
})
if(S1=="logistic" && S2=="logistic")
{(AMPmixed=logistic.logistic(y, x, data, start, random))}
else if(S1=="logistic" &&
S2=="asymptotic"){(AMPmixed=logistic.asymptotic(y, x, data, start,
random))}
else if(S1=="asymptotic" &&
S2=="logistic"){(AMPmixed=asymptotic.logistic(y, x, data, start,
random))}
else if(S1=="asymptotic" &&
S2=="asymptotic"){(AMPmixed=asymptotic.asymptotic(y, x, data, start,
random))}
}
#
con rep biomass
10.00 1 1.126
20.32 1 1.096
31.00 1 1.163
43.20 1 0.985
5 10.00 1 0.716
6 32.00 1 0.560
7 100.00 1 0.375
80.00 2 0.833
90.32 2 1.106
10 1.00 2 1.336
11 3.20 2 0.754
12 10.00 2 0.683
13 32.00 2 0.488
14 100.00 2 0.344
iso<-read.table(file="E:\\Hormesis\\data\\isobutanol.txt", header=T)
aa<-groupedData(biomass~con|rep, data=iso)
van2<-AMPmixed(y=biomass, x=con, S1="asymptotic", S2="asymptotic", data=aa,
random=pdDiag(alpha+delta+lpsi.a+gamma+h~1),
start=c(alpha= 0.7954, delta= 0.3231, lpsi.a=-0.2738,
gamma= 1.0366, h=0.8429))
van2
Thank you very much in advance.
Chunhao
__
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] How to make R running faster
Hi everyone, I run the R loops on window XP and vista. Both are Intel core 2 Duo 2.2 GHz with 2 GB ram and XP is significantly faster than vista. Dose anyone know how speed up R loops in vista? Thank you in advance. Chunhao Tu __ 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] test for multivariate normality?
Yes. You could install "mvnormtest" Package and perform the multivariate normality test. By using mshapiro.test I wish this is helpful! Chunhao Tu Quoting HongSheng Liao <[EMAIL PROTECTED]>: My stat textbook tells me that using Shapiro-Wilk test for each variable one by one is not equal to a test for multivariate normality as a whole. Does R have a function of testing for multivariate normality? Thanks. Hongsheng (Hank) Liao, Ph.D. Lab Manager Center for Quantitative Fisheries Ecology 800 West 46th Street Old Dominion University Norfolk, Virginia 23508 Phone:757.683.4571 Fax:757.683.5293 __ 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] How do you exit a function in R?
You might try to use "on.exit" or "stop"?
# on.exit
if (nAssetPositions != nAssetPrices) {
on.exit(cat("Different number of assets! "\n"))
}
# stop
if (nAssetPositions != nAssetPrices) {
stop("Different number of assets!")
}
You could find these in "S programming" by W.N Venables and B.D. Ripley
Chunhao
Quoting Bill Cunliffe <[EMAIL PROTECTED]>:
For example, based on a certain condition, I may want to exit my code early:
# Are there the same number of assets in "prices" and
"positions"?
if (nAssetPositions != nAssetPrices) {
cat("Different number of assets! \n\n")
}
I have searched, but not found, a way of forcing a function to exit. Any
help greatly appreciated.
[[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.
Re: [R] Thank you
I totally agree both of you. This is a super place to mature the R. I learn a lot from this R heaven! Chunhao Quoting Esmail Bonakdarian <[EMAIL PROTECTED]>: Tubin wrote: In the past few weeks I have had to give myself a crash course in R, in order to accomplish some necessary tasks for my job. During that time, I've found this forum to be helpful time and time again - usually I find the answer to my problem by searching the archives; once or twice I've posted questions and been given near-immediate solutions. I just wanted to thank the forum participants for all of their contributions over the years! Hear hear! .. I've benefited greatly from reading the postings here and some of the members have been very generous with their knowledge too! Esmail __ 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] optimization setup
Hi,
I think Ray has answered this question in the previous e-mail.
Because optim can only use one single parameter thus you can not have
the parameters: theta, theta1 and x at the same time.
such as:
f1<-function(theta)
{theta[1]+theta[2]}
f2<-function(theta)
{f1(theta)*3}
f3<-function(theta)
{f1(theta)-f2(theta)}
optim(par=c(1,1), f3)
Chunhao
Quoting threshold <[EMAIL PROTECTED]>:
Hi, thanks for your replay the previous post of mine. Let me ask you one more
question similar to the previous one, based on a naive example:
f1<-function(theta, theta1)
{theta[1]+theta[2]+theta1[1]}
f2<-function(theta, x, theta1)
{f1(theta, theta1)*exp(x)*theta1[2]}
function to be optimized with respect to theta1:
f3<-function(theta1)
{f1(theta, theta1)-f2(theta,x,theta1)}
optim(par=c(1,1), f3)
If I do so I get:
Error in f1(theta, theta1) : argument "theta1" is missing, with no default
What I did wrong? I guess I miss something basic...
Thanks in advance, robert
Ray Brownrigg-2 wrote:
On Wed, 23 Apr 2008, threshold wrote:
Hi, here comes my problem, say I have the following functions (example
case) #
function1 <- function (x, theta)
{a <- theta[1] ( 1 - exp(-theta[2]) ) * theta[3] )
b <- x * theta[1] / theta[3]^2
return( list( a = a, b = b )) }
#---
function2<-function (x, theta)
{P <- function1(x, theta)
c <- P$a * x * exp(-theta[2])
d <- P$b * exp(x)
q <- theta[1] / theta[3]
res <- c + d + q; res}
# Function to be optimized
function3 <- function(theta1,theta2,theta3) {
n <- length(data)
-sum( function2(x = data[2:n], theta = c(theta1, theta2, theta3) ))}
# 'data' is my input ts class object
#--
Then I want to maximize function3 with respect to theta(s) (given some
starting values)
fit<-optim(par=c(theta1=1, theta2=1.2, theta3=.2), fn=function3)
I get the following:
Error in function1(x, theta) :
argument "theta2" is missing, with no default
Where I made a mistake? I will appreciate any help ...
r
Your function to be optimised must be a function of a single parameter
(which
may be a vector).
So you should have something like:
# Function to be optimized
function3 <- function(thetas) {
n <- length(data)
-sum( function2(x = data[2:n], theta = thetas))
}
[Not tested, your provided code has typos].
Ray
__
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.
--
View this message in context:
http://www.nabble.com/optimization-setup-tp16825410p17543578.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] optimization setup
I just wonder that if you know theta and x. Why can't you specify
these in your functions. such as
f1<-function(theta1){1+2+theta1[1]} (Let me assume theta[1]=1, theta[2]=2)
f2<-function(theta1){f1(theta1)*exp(3)*theta1[2]} (again, I assume x=3)
f3<-function(theta1){f1(theta1)-f2(theta1)}
optim(par=c(1.1, 2.1), f3)
If you know the exact values of theta you should specify in advance
otherwise optim will try to get the estimated theta[1] and theata[2].
In this case, you might need to use loop for your x vector.
If I am wrong please correct me!
Thanks
Chunhao
Quoting Paul Smith <[EMAIL PROTECTED]>:
On Fri, May 30, 2008 at 12:04 PM, threshold <[EMAIL PROTECTED]> wrote:
Thanks for all your replies and sorry for a negligence in my examples. They
are very simplified to reflect very roughly the structure of the case I deal
with, which is too complicated to be quoted here.
What I deal with is:
1) 'theta' which is vector with length 2, being known to me (eg.
theta=c(1,2))
2) 'x' which represents my empirical data (I know as well),
3) 'theta1' is a vector with length 2 I don't know, which suppose to
minimize my objective function f3 (see below)
Here come the structure of the functions:
f1<-function(theta, theta1)
{theta[1]+theta[2]+theta1[1]}
f2<-function(theta, x, theta1)
{f1(theta, theta1)*exp(x)*theta1[2]}
function to be optimized with respect to theta1:
f3<-function(theta1)
{f1(theta, theta1)-f2(theta,x,theta1)}
Again, I know vector 'theta' and 'x', and I look for (vector) 'theta1' which
minimize f3
Question are:
1) whether, given the case I deal with, the functions I provided are
specified correctly. If not what is the correct form?
2) how should I write my optim function with starting values for theta1
equal to 1.1 and 2.1. What I did was:
optim(par=c(1.1, 2.1), f3) but it did not work well with error message:
Error in f1(theta, theta1) : argument "theta1" is missing, with no default
If x is a vector, then your function f2 is not a real function; f2
returns a vector (and not a scalar). Optim does not optimize vectorial
functions.
Paul
__
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] function
I wrote a function and I hope this helps
b<-c(1,2,3,4);lambda<-c(0.3,0.4,0.5,0.6)
m<-c(0.10,0.11,0.12,0.13);t<-c(7,8,9,10)
N=4
Y<-matrix(data=NA,nrow=4, ncol=1)
for (i in 1:N){
y=b*(1-exp(lambda*(t^m)))
Y<-matrix(data=y, nrow=4, ncol=1)
Y<-as.numeric(Y)
sim<-data.frame(t, Y);plot(sim,xlab="age", ylab="response", col="red")
}
Chunhao
Quoting admin <[EMAIL PROTECTED]>:
hi i have to make my function
y=b*(1-exp(lambda*(t^m))
and my data are in excel file with names:
Year /b/lambda/m
and t is an age in another file calls pcb
could u help me please to calculate all y for all years please
i can make this step by step but i want an program for having all y and plot
them in the same windows with legend could u help me please thank u
[[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.
Re: [R] How to solve a non-linear system of equations using R
Hi Jorge,
Have you tried to use "systemfit" package. In this package, this is a
function call " nlsystemfit ". This might help.
Chunhao
Quoting Jorge Ivan Velez <[EMAIL PROTECTED]>:
Dear R-list members,
I've had a hard time trying to solve a non-linear system (nls) of equations
which structure for the equation i, i=1,...,4, is as follows:
f_i(d_1,d_2,d_3,d_4)-k_i(l,m,s) = 0(1)
In the expression above, both f_i and k_i are known functions and l, m and s
are known constants. I would like to estimate the vector d=(d_1,d_2,d_3,d_4)
which is solution of (1). Functions in R to estimate f_i-k_i are at the end
of this message.
Any help/suggestions/comments would be greatly appreciated.
Thanks in advance,
Jorge
# --
# Constants
# --
l=1
m=0.4795
s=0.4795
# --
# Functions to estimate f_i-k_i
# --
f1=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]
res=2*d1+2*sqrt(2)*d1*d2+2*sqrt(3)*d2*d3+4*d3*d4-l*m*(1+d1^2+d2^2+d3^2+d4^2)
res
}
f2=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]
res=2*sqrt(2)*d2+2*d1^2+2*sqrt(6)*d1*d3+4*d2^2+4*sqrt(3)*d2*d4+6*d3^2+8*d4^2-l*(m^2+m^3*s^(-1))*(1+d1^2+d2^2+d3^2+d4^2)
res
}
f3=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]
res=6*d1+12*sqrt(2)*d1*d2+18*sqrt(3)*d2*d3+48*d3*d4+2*sqrt(6)*d3+4*sqrt(6)*d1*d4-l*(m^3+3*m^4*s^(-1)+3*m^6*s^(-2))*(1+d1^2+d2^2+d3^2+d4^2)
res
}
f4=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]
res=12*sqrt(2)*d2+12*d1^2+36*d2^2+72*d3^2+120*d4^2+20*sqrt(6)*d1*d3+56*sqrt(3)*d2*d4+4*sqrt(6)*d4-l*((m^4+6*m^6*s^(-1)+15*m^6*s^(-2)+15*m^7*s^(-3))-3*l^(2)*(m^2+m^3*s^(-1))^2)*(1+d1^2+d2^2+d3^2+d4^2)
res
}
[[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.
Re: [R] splitting data frame based on a criteria
Or ?gsummary Quoting Moshe Olshansky <[EMAIL PROTECTED]>: Also check ?aggregate --- On Wed, 4/6/08, Ingmar Visser <[EMAIL PROTECTED]> wrote: From: Ingmar Visser <[EMAIL PROTECTED]> Subject: Re: [R] splitting data frame based on a criteria To: "Marvin Lists" <[EMAIL PROTECTED]> Cc: r-help@r-project.org Received: Wednesday, 4 June, 2008, 5:11 AM ?by may be helpful here eg if dat is your data.frame and yf is a factor (created using ifelse) use by(dat,yf,mean) to compute the means for each level of yf hth, Ingmar On Jun 3, 2008, at 8:37 PM, Marvin Lists wrote: > Hi, > I have a data frame that I want to split into two based on the > values of a > variable in it. > > The variable Y has numeric values ranging between 0 through 70. I > want to > plot the frequencies of another variable X in two different cases: > - When Y = 0 and > - When Y > 0 > > How does one go about doing this? > > In general, I want to do several analyses with this data frame that > are a > variation of the above situation, i.e. they require splitting the > data into > different age, gender etc. and then calculating separate means, > correlations > and so on for the different groups into which the data frame would > split. > > I am struggling with the correct syntax for achieving this. > > Reading through the documentation suggests that tapply and split > may be the > functions to use for my purposes but the examples in the > documentations > didn't help me understand how I could achieve this. > > I would appreciate any suggestions and help. > > Thanks, > Marvin > >[[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. __ 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] sorting the data~
id<-c(1,1,1,1,3,3,3,7,7,7,7,11,11,11,2,2,2,4,4,4,4,8,8,8)
sort(id)
[1] 1 1 1 1 2 2 2 3 3 3 4 4 4 4 7 7 7 7 8 8 8 11 11 11
Quoting Manli Yan <[EMAIL PROTECTED]>:
no,the id is variable of a table,such as:
treatment id age response
low 1 50 20
low 1 60 30
high5 50 30
high5 60 40
...
I want to rearranage the table according the id (increasing),since id is not
strictly from 1~n,it is in increasing order but sometime jump through many
number like 1 1 5 5,I like them to be 1 1 2 2~
2008/6/4 Erik Iverson <[EMAIL PROTECTED]>:
Are these the ranks of the data?
help.search("rank")
Manli Yan wrote:
id<-c(1,1,1,1,3,3,3,7,7,7,7,11,11,11)
how to sort this kind of data to
id:(1,1,1,1,2,2,2,3,3,3,3,4,4,4.)
thanks~
[[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-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] mean
TABLE<-matrix(data=c(12,13,14,15,24,10),byrow=T,nrow=2,ncol=3) TABLE [,1] [,2] [,3] [1,] 12 13 14 [2,] 15 24 10 apply(TABLE,1,mean) [1] 13.0 16.3 Chunhao Quoting Marco Chiapello <[EMAIL PROTECTED]>: Hi, I have a simple question. If I have a table and I want to have the mean for each row, how can I do?! Es: c1 c2 c3 mean 1 12 13 14 ?? 2 15 24 10 ?? ... Thanks, Marco __ 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.
