Re: [R] Change the colour of lines in the xyplot?
On 7/12/2010 2:52 AM, Jian Zhang wrote: I want to add the legend for my figure using auto.key in the xyplot function (library(lattice)). But I don't know how to change the colors of the lines in the legend part. I tried to add col=1:8 in the auto.key, but it changes the colors of the text of my Legend, not the lines. How can I change the colors of the lines? xyplot(estimated~year,data=res,type=l,group=sp, auto.key= list(points = FALSE, lines=TRUE, corner =c(0.1,0.6),size=5, col=1:8)) xyplot(estimated ~ year, data=res, type=l, group=sp, par.settings = list(superpose.line=list(col=1:8)), auto.key=list(points=FALSE, lines=TRUE, corner=c(0.1,0.6), size=5)) [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 with multiple imputation
On 6/30/2010 1:14 AM, Daniel Chen wrote:
Hi,
I am a long time SPSS user but new to R, so please bear with me if my
questions seem to be too basic for you guys.
I am trying to figure out how to analyze survey data using logistic
regression with multiple imputation.
I have a survey data of about 200,000 cases and I am trying to predict the
odds ratio of a dependent variable using 6 categorical independent variables
(dummy-coded). Approximatively 10% of the cases (~20,000) have missing data
in one or more of the independent variables. The percentage of missing
ranges from 0.01% to 10% for the independent variables.
My current thinking is to conduct a logistic regression with multiple
imputation, but I don't know how to do it in R. I searched the web but
couldn't find instructions or examples on how to do this. Since SPSS is
hopeless with missing data, I have to learn to do this in R. I am new to R,
so I would really appreciate if someone can show me some examples or tell me
where to find resources.
Here is an example using the Amelia package to generate imputations
and the mitools and mix packages to make the pooled inferences.
titanic -
read.table(http://lib.stat.cmu.edu/S/Harrell/data/ascii/titanic.txt;,
sep=',', header=TRUE)
set.seed(4321)
titanic$sex[sample(nrow(titanic), 10)] - NA
titanic$pclass[sample(nrow(titanic), 10)] - NA
titanic$survived[sample(nrow(titanic), 10)] - NA
library(Amelia) # generate multiple imputations
library(mitools) # for MIextract()
library(mix) # for mi.inference()
titanic.amelia - amelia(subset(titanic,
select=c('survived','pclass','sex','age')),
m=10, noms=c('survived','pclass','sex'),
emburn=c(500,500))
allimplogreg - lapply(titanic.amelia$imputations,
function(x){glm(survived ~ pclass + sex + age, family=binomial, data = x)})
mice.betas.glm - MIextract(allimplogreg, fun=function(x){coef(x)})
mice.se.glm - MIextract(allimplogreg, fun=function(x){sqrt(diag(vcov(x)))})
as.data.frame(mi.inference(mice.betas.glm, mice.se.glm))
# Or using only mitools for pooled inference
betas - MIextract(allimplogreg, fun=coef)
vars - MIextract(allimplogreg, fun=vcov)
summary(MIcombine(betas,vars))
Thank you!
Daniel
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] If-error-resume-next
On 6/24/2010 10:30 AM, Christofer Bogaso wrote: In VBA there is a syntax if error resume next which sometimes acts as very beneficial on ignoring error. Is there any similar kind of function/syntax here in R as well? ?try Thanks for your attention [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Power calculation
On 6/10/2010 8:26 AM, Samuel Okoye wrote: Hello, Is there any R function which does power calculation for unbalanced groups (n1 neq n2)? Since power.t.test has n Number of observations (per group). Many thanks, Samuel See pwr.t2n.test() in the pwr package. http://finzi.psych.upenn.edu/R/library/pwr/html/pwr.t2n.test.html You might have found it by doing the following in R: RSiteSearch(power analysis, restrict='function') [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] switching elements of a vector
On 5/24/2010 6:02 AM, speretti wrote: Hi, I would like to receive help for the following matter: If I'm dealing with a numeric vectors containing increasing elements. i.e. a-c(1,2,2,2,2,3,3,3,4,4,4,5,5,6,7,7,7) There exist an efficient way to obtain an vector that indicates the position of the changing element of a? In this case it would be something like: index-c(1,6,9,12,14,15) a - c(1,2,2,2,2,3,3,3,4,4,4,5,5,6,7,7,7) rle(a) Run Length Encoding lengths: int [1:7] 1 4 3 3 2 1 3 values : num [1:7] 1 2 3 4 5 6 7 cumsum(head(rle(a)$lengths, -1)) + 1 [1] 2 6 9 12 14 15 ?rle usually I'm used cycles to obtain boolean vectors of the same length of a indicating the changing elements ...later I've muliplied them for their numeric sequence and after that I've selected elements different from zero ...it is quite long... can you find an easier solution? Thank you for you help -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Getting dates in an SPSS file in right format.
On 5/18/2010 12:38 PM, Praveen Surendran wrote: Dear all, I am trying to read an SPSS file into a data frame in R using method read.spss(), sample - read.spss(file.name,to.data.frame=TRUE) But dates in the data.frame 'sample' are coming as integers and not in the actual date format given in the SPSS file. Appreciate if anyone can help me to solve this problem. Date variables in SPSS contain the number of seconds since October 14, 1582. You might try something like this: sample$MYDATE - as.Date(as.POSIXct(sample$MYDATE, origin=1582-10-14, tz=GMT)) Kind Regards, Praveen Surendran 2G, Complex and Adaptive Systems Laboratory (UCD CASL) School of Medicine and Medical Sciences University College Dublin Belfield, Dublin 4 Ireland. Office : +353-(0)1716 5334 Mobile : +353-(0)8793 13071 [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 ls() only functions or anything else but functions?
On 5/13/2010 10:02 AM, John Edwards wrote: Hello, How ls() only functions or only data objects (basically anything other than functions) such as data.frame, numeric ...? c(ls.str(mode = function)) ls()[!(ls() %in% c(ls.str(mode=function)))] ?ls.str John [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 ls() only functions or anything else but functions?
On 5/13/2010 11:04 AM, Chuck Cleland wrote: On 5/13/2010 10:02 AM, John Edwards wrote: Hello, How ls() only functions or only data objects (basically anything other than functions) such as data.frame, numeric ...? c(ls.str(mode = function)) ls()[!(ls() %in% c(ls.str(mode=function)))] ?ls.str Or better ... c(lsf.str()) ls()[!(ls() %in% c(lsf.str()))] John [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] vectorize a power analysis?
On 5/12/2010 3:34 PM, Jack Siegrist wrote: We are doing a power analysis by generating noisy data sets according to a model, fitting the model to the data, and extracting a p-value. What is the best way to do this many times? We are just using for loops and it is too slow because we are repeating the analysis for many parameterizations. I can think of several ways to do this: for loop sapply using the plyr package using the lme4 package You don't mention replicate(), which I would consider. For example: replicate(10, t.test(rnorm(20))$p.value) [1] 0.2622419 0.1538739 0.9759340 [4] 0.7413474 0.1541895 0.4321595 [7] 0.5800549 0.7329299 0.9625038 [10] 0.1315875 If you write a function that does the data generating, model fitting, and p-value extraction, then replicate can run that function many times. I don't know how the timing compares, but I like the simplicity and readability of replicate(). hope this helps, Chuck Someone told me that the apply functions are barely any faster than for loops, so what is the best way, in general, to approach this type of problem in R-style? Could someone point to a discussion of the comparative time efficiencies of these and other appropriate methods? I'm not looking for specific code, just sort of a list of the common approaches with information about efficiency. Thanks, Jack -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Determining Index of Last Element in Vector
On 4/25/2010 2:10 PM, Alan Lue wrote: Hi, Is there a way to specify the last element of a vector, similar to end in MATLAB? v[end] would be MATLAB for v(length(v)) in R. While `v(length(v))' does yield the last element, that approach fails in the following, rep(v, each=2)[-c(1,length(v))] which is meant to duplicate all elements of `v' except for the first and last. (I.e., if `v - 1:4', then we want '1 2 2 3 3 4'.) v - 1:4 rep(v, c(1, rep(2, length(v) - 2), 1)) [1] 1 2 2 3 3 4 So the question is, is there a better way specify the last element of a vector? If not, is there a better way to duplicate all elements of a vector except for the first and last? (I know you can achieve this using two lines, but I'm writing because I want to do it using one.) Alan -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] substract start from the end of the vector
with(df, substr(DESCRIPTION, start=1, stop=nchar(DESCRIPTION) - 10)) ?nchar On 4/23/2010 11:57 AM, arnaud Gaboury wrote: Dear group, Here is my df : df - structure(list(DESCRIPTION = c(PRM HGH GD ALUMINIUM USD 09/07/10 , PRM HGH GD ALUMINIUM USD 09/07/10 , PRIMARY NICKEL USD 04/06/10 ), CREATED.DATE = structure(c(18361, 18361, 18325), class = Date), QUANITY = c(-1L, 1L, 1L), CLOSING.PRICE = c(2,415.90, 2,415.90, 25,755.71)), .Names = c(DESCRIPTION, CREATED.DATE, QUANITY, CLOSING.PRICE), row.names = c(NA, 3L), class = data.frame) df DESCRIPTION CREATED.DATE QUANITY CLOSING.PRICE 1 PRM HGH GD ALUMINIUM USD 09/07/102020-04-09 -1 2,415.90 2 PRM HGH GD ALUMINIUM USD 09/07/102020-04-091 2,415.90 3 PRIMARY NICKEL USD 04/06/102020-03-041 25,755.71 In the DESCRIPTION column, I want to get rid of the date (09/07/10 .). I know the function substr(x, start, stop), but in my case, I need to indicate I want to start from the end of the vector, and I have no idea how to pass this argument. TY for any help *** Arnaud Gaboury Mobile: +41 79 392 79 56 BBM: 255B488F *** [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 make a boxplot with exclusion of certain groups
On 4/19/2010 12:21 PM, josef.kar...@phila.gov wrote:
This seems like a simple thing, but I have been stuck for some time. My
data has 2 columns. Column 1 is the value, and column 2 is the Site where
data was collected. Column 2 contains 7 different Sites (i.e. groups). I
am only interested in showing 3 groups on a single boxplot.
I have tried various methods of subsetting the data, in order to only have
the 3 groups in my subset. However even after doing this, all 7 groups
carry forward, so that when I make a boxplot of my subsetted data, all 7
groups still appear in the x-axis labels; all 7 groups also appear in the
boxplot summary (i.e. the values returned with boxplot (…plot=FALSE) ) .
Even if I delete the unwanted groups from the ‘levels’ of Column 2, they
still appear on the plot, and in the boxplot summary statistics.
There are various tricks I can do with the boxplot summary statistics to
correct for this, but they get complicated when I want to change the
algorithm for calculating outliers and their corresponding groups. Rather
than do all these tricks, it seems much simpler to fully exclude the
unwanted groups from the beginning. But this doesn’t appear to work
Any ideas?
library(gdata) # for drop.levels()
DF - data.frame(site = rep(LETTERS[1:7], each=20), y = runif(7*20))
boxplot(y ~ drop.levels(site), data=subset(DF, site %in% c('A','D','F'),
drop=TRUE))
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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 delete columns with NA values?
On 4/14/2010 10:56 AM, muting wrote: Hi everyone: I have a dataset: tm1 col1 col2 [1,]1 NA [2,]11 [3,]22 [4,]11 [5,]22 [6,]1 NA I need to delete entire column 2 that has NA in it(not all of them are NAs), and the result I want is tm1 col1 [1,]1 [2,]1 [3,]2 [4,]1 [5,]2 [6,]1 what should I do? subset(tm1, select=colMeans(is.na(tm1)) == 0) OR tm1[,colMeans(is.na(tm1)) == 0] I search a lot, all I found is how to delete column with all NA values.. Thanks a lot muting -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Finding common an unique elements in character vectors
On 3/29/2010 1:53 PM, Thomas Jensen wrote:
Dear R-list,
I have a problem which I think is quite basic, but so far google has not
helped me.
I have two vectors like this:
vector_1 - c(Belgium, Spain, Greece, Ireland, Luxembourg, Netherlands,
Portugal)
vector_2 - c(Denmark, Luxembourg)
I would like to find the elements in vector_1 that are not in vector_2
so that i get a vector with these countries: Belgium, Spain, Greece,
Ireland, Netherlands, Portugal.
vector_1 - c('Belgium', 'Spain', 'Greece', 'Ireland', 'Luxembourg',
'Netherlands', 'Portugal')
vector_2 - c('Denmark', 'Luxembourg')
setdiff(vector_1, vector_2)
?setdiff
Thanks a lot,
Thomas
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Stacking matrices
On 3/15/2010 1:37 PM, jlu...@ria.buffalo.edu wrote: What is an easy way to stack a matrix multiple times? E.g. I have a 6x6 matrix that I need to stack vertically 154 times to get a 6*154 by 6 matrix. I would rather not rbind(X,X,...,X) matrices.--Joe X - matrix(runif(36), ncol=6) matrix(rep(t(X), 154), ncol=6, byrow=TRUE) [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Randomly sampling subsets of dataframe variable
On 3/12/2010 3:06 PM, Hosack, Michael wrote: Fellow R users, I am stumped on what would seem to be something fairly simple. I have a dataframe that has a variable named 'WEEK' that takes the numbers 1:26 (26 week time-period) with each number repeated five times consecutively (once for each weekday, Monday through Friday). Ex. 123.2626262626. I would like to randomly extract two weekdays per five day week for each of 26 weeks and store this data as a separate dataframe. I have been unable to get the sample function to work properly. I have also tried using the runif function to assign random numbers to each row of my dataframe, sort the dataframe first by week number then by random number value, and finally select the first two elements from each week subset (26 weeks total, giving 52 randomly selected values). I can't figure out how to select the first two elements. My goal is to randomly select two weekdays per week (without replacement) for each of 26 consecutive weeks. Any advice would be greatly appreciated. DF - data.frame(WEEK = rep(1:26, each=5), DAY = rep(1:5, 26), X = runif(5*26)) DF2 - data.frame(DAY = c(replicate(26, sample(5, 2, replace=FALSE))), WEEK = rep(1:26, each=2)) new.DF - merge(DF, DF2, all=FALSE) Thank you, Mike __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] TukeyHSD model thing
On 3/6/2010 4:38 PM, casperyc wrote: Hi, I am trying to reproduce a tukey test in R == x=c(145,40,40,120,180, 140,155,90,160,95, 195,150,205,110,160, 45,40,195,65,145, 195,230,115,235,225, 120,55,50,80,45 ) y2=c( rep(as.character(1),5), rep(as.character(2),5), rep(as.character(3),5), rep(as.character(4),5), rep(as.character(5),5), rep(as.character(6),5) ) crd2=data.frame(x,y2) model1=aov(x~y2,data=crd2) TukeyHSD(model1) == The result in the 'diff' of the means are the same as those did using SAS, (which is in my tutorial sheet, i got a MAC, so cant use SAS) however, the 95% confiden limits are quite different. === 2-1 23 -73.817441 119.817441 0.975518208 3-1 59 -37.817441 155.817441 0.435116628 4-1 -7 -103.817441 89.817441 0.12318 5-1 95 -1.817441 191.817441 0.056613465 6-1 -35 -131.817441 61.817441 0.869242006 3-2 36 -60.817441 132.817441 0.855531189 4-2 -30 -126.817441 66.817441 0.926612938 5-2 72 -24.817441 168.817441 0.232896275 6-2 -58 -154.817441 38.817441 0.453535553 4-3 -66 -162.817441 30.817441 0.316718292 5-3 36 -60.817441 132.817441 0.855531189 6-3 -94 -190.817441 2.817441 0.060579795 5-4 1025.182559 198.817441 0.034819938 6-4 -28 -124.817441 68.817441 0.944203446 6-5 -130 -226.817441 -33.182559 0.004294761 === in the SAS output, it's (in slightly different order, you can just check one of the set) === 5 - 3 36.00 -28.63 100.63 5 - 2 72.00 7.37 136.63 *** 5 - 1 95.00 30.37 159.63 *** 5 - 4 102.00 37.37 166.63 *** 5 - 6 130.00 65.37 194.63 *** 3 - 5 -36.00 -100.63 28.63 3 - 2 36.00 -28.63 100.63 3 - 1 59.00 -5.63 123.63 3 - 4 66.00 1.37 130.63 *** 3 - 6 94.00 29.37 158.63 *** 2 - 5 -72.00 -136.63 -7.37 *** 2 - 3 -36.00 -100.63 28.63 2 - 1 23.00 -41.63 87.63 2 - 4 30.00 -34.63 94.63 2 - 6 58.00 -6.63 122.63 1 - 5 -95.00 -159.63 -30.37 *** 1 - 3 -59.00 -123.63 5.63 1 - 2 -23.00 -87.63 41.63 1 - 4 7.00 -57.63 71.63 1 - 6 35.00 -29.63 99.63 4 - 5 -102.00 -166.63 -37.37 *** 4 - 3 -66.00 -130.63 -1.37 *** 4 - 2 -30.00 -94.63 34.63 4 - 1 -7.00 -71.63 57.63 4 - 6 28.00 -36.63 92.63 6 - 5 -130.00 -194.63 -65.37 *** 6 - 3 -94.00 -158.63 -29.37 *** 6 - 2 -58.00 -122.63 6.63 6 - 1 -35.00 -99.63 29.63 6 - 4 -28.00 -92.63 36.63 === say, betweet treatment 5 and 3 R 5-3 36 -60.817441 132.817441 SAS 5-3 36 -28.63 100.63 i am wondering if i have done something wrong in R? It looks like the confidence intervals for your SAS output are not family-wise intervals. For example, you can get intervals that match your SAS output when you don't adjust for multiple comparisons. confint(lm(x ~ as.factor(y2), data=crd2)) 2.5 %97.5 % (Intercept) 59.302005 150.69800 as.factor(y2)2 -41.626725 87.62672 as.factor(y2)3 -5.626725 123.62672 as.factor(y2)4 -71.626725 57.62672 as.factor(y2)5 30.373275 159.62672 as.factor(y2)6 -99.626725 29.62672 full documents are in the attachment, can someone suggest to me the relevent R codes that does the same sort of thing? (tukeyHSD,fisherLSD,and anova table ) Thanks! casper http://n4.nabble.com/file/n1583109/SAS.pdf SAS.pdf http://n4.nabble.com/file/n1583109/R.pdf R.pdf http://n4.nabble.com/file/n1583109/ws1.pdf ws1.pdf http://n4.nabble.com/file/n1583109/ws1sols.pdf ws1sols.pdf -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] sem package and growth curves
On 3/2/2010 1:43 AM, Daniel Nordlund wrote:
I have been working through the book Applied longitudinal data analysis:
modeling change and event occurrence by Judith D. Singer and John B.
Willett. I have been working examples using SAS and also using it as an
opportunity for learning to use R for statistical analysis.
I ran into some difficulties in chapter 8 which deals with using structural
equation modeling. I have tried to use the sem package to replicate the
problem solutions in chapter 8. I am more familiar with RAM specifications
than I am with structural equations (though I am a novice at both). The
solutions I have tried seem to be very sensitive to starting values
(especially with more complex models). I don't know if this is just my lack
of knowledge in this area, or something else.
Has anyone worked out solutions to the Singer and Willett examples for
Chapter 8 that they would be willing to share? I would also be interested in
other simple examples using sem and RAM specifications. If anyone is
interested, I would also be willing to share the R code I have written for
other chapters in the Singer and Willett book.
Hi Dan,
See below for my code for Models A-D in Chapter 8. As you point out,
I find that this only works when good starting values are given. I took
the starting values from the results given for another program (Mplus)
at the UCLA site for this text:
http://www.ats.ucla.edu/stat/examples/alda.htm
I greatly appreciate John Fox's hard work on the sem package, but
since good starting values will generally not be available to applied
users I think the package is not as useful for these types of models as
it could be. If anyone has approaches to specifying the models that are
less sensitive to starting values, or ways for less sophisticated users
to generate good starting values, please share.
Chuck
# Begin Code for Models A-D, Chapter 8, Singer Willett (2003)
alc2 -
read.table(http://www.ats.ucla.edu/stat/mplus/examples/alda/alcohol2.txt;,
sep=\t, header=FALSE)
names(alc2) - c('ID','FEMALE','ALC1','ALC2','ALC3','PEER1','PEER2','PEER3')
alc2$UNIT - 1
library(sem)
alc2.modA.raw - raw.moments(subset(alc2,
select=c('ALC1','ALC2','ALC3','UNIT')))
modA - specify.model()
I - ALC1, NA, 1
I - ALC2, NA, 1
I - ALC3, NA, 1
S - ALC1, NA, 0
S - ALC2, NA, 0.75
S - ALC3, NA, 1.75
UNIT - I, Mi, 0.226
UNIT - S, Ms, 0.036
I - I, Vi, NA
S - S, Vs, NA
I - S, Cis, NA
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
sem.modA - sem(modA, alc2.modA.raw, 1122, fixed.x=UNIT, raw=TRUE)
summary(sem.modA)
alc2.modB.raw - raw.moments(subset(alc2,
select=c('FEMALE','ALC1','ALC2','ALC3','UNIT')))
modB - specify.model()
I - ALC1, NA, 1
I - ALC2, NA, 1
I - ALC3, NA, 1
S - ALC1, NA, 0
S - ALC2, NA, 0.75
S - ALC3, NA, 1.75
FEMALE - I, B1, NA
FEMALE - S, B2, NA
UNIT - I, Mi, 0.226
UNIT - S, Ms, 0.036
I - I, Vi, NA
S - S, Vs, NA
I - S, Cis, NA
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
sem.modB - sem(modB, alc2.modB.raw, 1122, fixed.x=c(FEMALE,UNIT),
raw=TRUE)
summary(sem.modB)
alc2.modC.raw - raw.moments(subset(alc2,
select=c('FEMALE','ALC1','ALC2','ALC3','UNIT')))
modC - specify.model()
I - ALC1, NA, 1
I - ALC2, NA, 1
I - ALC3, NA, 1
S - ALC1, NA, 0
S - ALC2, NA, 0.75
S - ALC3, NA, 1.75
FEMALE - I, B1, NA
FEMALE - S, NA, 0
UNIT - I, Mi, 0.226
UNIT - S, Ms, 0.036
I - I, Vi, NA
S - S, Vs, NA
I - S, Cis, NA
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
sem.modC - sem(modC, alc2.modC.raw, 1122, fixed.x=c(FEMALE,UNIT),
raw=TRUE)
summary(sem.modC)
alc2.modD.raw - raw.moments(subset(alc2,
select=c('PEER1','PEER2','PEER3','ALC1','ALC2','ALC3','UNIT')))
modD - specify.model()
Ia - ALC1, NA, 1
Ia - ALC2, NA, 1
Ia - ALC3, NA, 1
Sa - ALC1, NA, 0
Sa - ALC2, NA, 0.75
Sa - ALC3, NA, 1.75
UNIT - Ia, Mia, 0.226
UNIT - Sa, Msa, 0.036
Ip - PEER1, NA, 1
Ip - PEER2, NA, 1
Ip - PEER3, NA, 1
Sp - PEER1, NA, 0
Sp - PEER2, NA, 0.75
Sp - PEER3, NA, 1.75
Ip - Ia, B1, 0.799
Sp - Ia, B2, 0.080
Ip - Sa, B3, -0.143
Sp - Sa, B4, 0.577
UNIT - Ip, Mip, 0.226
UNIT - Sp, Msp, 0.036
Ia - Ia, Via, 0.042
Sa - Sa, Vsa, 0.009
Ia - Sa, Cisa, -0.006
Ip - Ip, Vip, 0.070
Sp - Sp, Vsp, 0.028
Ip - Sp, Cisp, 0.001
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
PEER1 - PEER1, Vd4, 0.106
PEER2 - PEER2, Vd5, 0.171
PEER3 - PEER3, Vd6, 0.129
ALC1 - PEER1, Cd1, 0.011
ALC2 - PEER2, Cd2, 0.034
ALC3 - PEER3, Cd3, 0.037
sem.modD - sem(modD, alc2.modD.raw, 1122, fixed.x=c(UNIT), raw=TRUE)
summary(sem.modD)
Thanks,
Dan
Daniel Nordlund
Bothell, WA USA
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th
Re: [R] sem package and growth curves
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
Behalf Of Chuck Cleland
Sent: March-03-10 9:03 AM
To: Daniel Nordlund
Cc: 'r-help'
Subject: Re: [R] sem package and growth curves
On 3/2/2010 1:43 AM, Daniel Nordlund wrote:
I have been working through the book Applied longitudinal data
analysis:
modeling change and event occurrence by Judith D. Singer and John B.
Willett. I have been working examples using SAS and also using it as an
opportunity for learning to use R for statistical analysis.
I ran into some difficulties in chapter 8 which deals with using
structural
equation modeling. I have tried to use the sem package to replicate the
problem solutions in chapter 8. I am more familiar with RAM
specifications
than I am with structural equations (though I am a novice at both). The
solutions I have tried seem to be very sensitive to starting values
(especially with more complex models). I don't know if this is just my
lack
of knowledge in this area, or something else.
Has anyone worked out solutions to the Singer and Willett examples for
Chapter 8 that they would be willing to share? I would also be interested
in
other simple examples using sem and RAM specifications. If anyone is
interested, I would also be willing to share the R code I have written for
other chapters in the Singer and Willett book.
Hi Dan,
See below for my code for Models A-D in Chapter 8. As you point out,
I find that this only works when good starting values are given. I took
the starting values from the results given for another program (Mplus)
at the UCLA site for this text:
http://www.ats.ucla.edu/stat/examples/alda.htm
I greatly appreciate John Fox's hard work on the sem package, but
since good starting values will generally not be available to applied
users I think the package is not as useful for these types of models as
it could be. If anyone has approaches to specifying the models that are
less sensitive to starting values, or ways for less sophisticated users
to generate good starting values, please share.
Chuck
# Begin Code for Models A-D, Chapter 8, Singer Willett (2003)
alc2 -
read.table(http://www.ats.ucla.edu/stat/mplus/examples/alda/alcohol2.txt;,
sep=\t, header=FALSE)
names(alc2) -
c('ID','FEMALE','ALC1','ALC2','ALC3','PEER1','PEER2','PEER3')
alc2$UNIT - 1
library(sem)
alc2.modA.raw - raw.moments(subset(alc2,
select=c('ALC1','ALC2','ALC3','UNIT')))
modA - specify.model()
I - ALC1, NA, 1
I - ALC2, NA, 1
I - ALC3, NA, 1
S - ALC1, NA, 0
S - ALC2, NA, 0.75
S - ALC3, NA, 1.75
UNIT - I, Mi, 0.226
UNIT - S, Ms, 0.036
I - I, Vi, NA
S - S, Vs, NA
I - S, Cis, NA
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
sem.modA - sem(modA, alc2.modA.raw, 1122, fixed.x=UNIT, raw=TRUE)
summary(sem.modA)
alc2.modB.raw - raw.moments(subset(alc2,
select=c('FEMALE','ALC1','ALC2','ALC3','UNIT')))
modB - specify.model()
I - ALC1, NA, 1
I - ALC2, NA, 1
I - ALC3, NA, 1
S - ALC1, NA, 0
S - ALC2, NA, 0.75
S - ALC3, NA, 1.75
FEMALE - I, B1, NA
FEMALE - S, B2, NA
UNIT - I, Mi, 0.226
UNIT - S, Ms, 0.036
I - I, Vi, NA
S - S, Vs, NA
I - S, Cis, NA
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
sem.modB - sem(modB, alc2.modB.raw, 1122, fixed.x=c(FEMALE,UNIT),
raw=TRUE)
summary(sem.modB)
alc2.modC.raw - raw.moments(subset(alc2,
select=c('FEMALE','ALC1','ALC2','ALC3','UNIT')))
modC - specify.model()
I - ALC1, NA, 1
I - ALC2, NA, 1
I - ALC3, NA, 1
S - ALC1, NA, 0
S - ALC2, NA, 0.75
S - ALC3, NA, 1.75
FEMALE - I, B1, NA
FEMALE - S, NA, 0
UNIT - I, Mi, 0.226
UNIT - S, Ms, 0.036
I - I, Vi, NA
S - S, Vs, NA
I - S, Cis, NA
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
sem.modC - sem(modC, alc2.modC.raw, 1122, fixed.x=c(FEMALE,UNIT),
raw=TRUE)
summary(sem.modC)
alc2.modD.raw - raw.moments(subset(alc2,
select=c('PEER1','PEER2','PEER3','ALC1','ALC2','ALC3','UNIT')))
modD - specify.model()
Ia - ALC1, NA, 1
Ia - ALC2, NA, 1
Ia - ALC3, NA, 1
Sa - ALC1, NA, 0
Sa - ALC2, NA, 0.75
Sa - ALC3, NA, 1.75
UNIT - Ia, Mia, 0.226
UNIT - Sa, Msa, 0.036
Ip - PEER1, NA, 1
Ip - PEER2, NA, 1
Ip - PEER3, NA, 1
Sp - PEER1, NA, 0
Sp - PEER2, NA, 0.75
Sp - PEER3, NA, 1.75
Ip - Ia, B1, 0.799
Sp - Ia, B2, 0.080
Ip - Sa, B3, -0.143
Sp - Sa, B4, 0.577
UNIT - Ip, Mip, 0.226
UNIT - Sp, Msp, 0.036
Ia - Ia, Via, 0.042
Sa - Sa, Vsa, 0.009
Ia - Sa, Cisa, -0.006
Ip - Ip, Vip, 0.070
Sp - Sp, Vsp, 0.028
Ip - Sp, Cisp, 0.001
ALC1 - ALC1, Vd1, 0.048
ALC2 - ALC2, Vd2, 0.076
ALC3 - ALC3, Vd3, 0.077
PEER1 - PEER1, Vd4, 0.106
PEER2 - PEER2, Vd5, 0.171
PEER3 - PEER3, Vd6, 0.129
ALC1 - PEER1, Cd1, 0.011
ALC2 - PEER2, Cd2, 0.034
ALC3 - PEER3, Cd3, 0.037
sem.modD - sem(modD, alc2.modD.raw, 1122
Re: [R] setting the steps for x axis labels on plot
On 3/1/2010 6:16 AM, Gonzalo Garcia-Perate wrote: Hello, I'm new to R, I've been working with it for the last 2 weeks. I am plotting some data and not getting the labels on the x axis I am expecting on my plot. my code reads #hours in the day h - c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23) #hp is a data frame with a pivot table of 25 columns (label and data for 24 hours) plot(h, as.matrix(hp[1,2:25]), main = hp[1,1], type=l, xlim=c(0,23), ylim=c(0,1800), xlab=hour, ylab=visitor activity) the result you can see here: http://www.flickr.com/photos/gonzillaaa/4397357491/ I was hoping the steps on the x axis would be every hour or at least every two hours I am getting unevenly spaced steps (0, 5, 10, 20) instead. apologies if this is too obvious a question, I've looked around on documentation etc. but found no reference to this. You can suppress the x axis in the call to plot() and then add it with a call to axis(). For, example plot(0:23, runif(24, min=0, max=1800), type=l, xlim=c(0,23), ylim=c(0,1800), xlab=hour, ylab=visitor activity, xaxt='n') axis(side=1, at=0:23, cex.axis=.5) ?axis Many thanks, Gonzalo. __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Subset and plot
On 2/2/2010 2:51 PM, Marlin Keith Cox wrote: Here is a runable program. When I plot Day and Wgt, it graphs all the data points. All I need is daily.sub1 plotted. I also need each Tanks to have its own col or pch. When I run it with the line with pch, it gives me nothing. DF - data.frame(Trial=rep(c(1,2),each=12), Tanks=rep(c(a3,a4,c4,h4),each=3,2), Day=rep(c(1:12),2), Wgt=c(1:24)) with(subset(DF, Trial==2 Tanks %in% c(a4,'c4','h4')), plot(Day, Wgt, pch=as.numeric(Tanks))) library(lattice) trellis.device(new=TRUE, color=FALSE) xyplot(Wgt ~ Day, groups=Tanks, data=subset(DF, Trial==2 Tanks %in% c(a4,'c4','h4')), par.settings=list(superpose.symbol=list(pch=c(2,19,21 rm(list=ls()) Trial-rep(c(1,2),each=12) Tanks=rep(c(a3,a4,c4,h4),each=3,2) Day=rep(c(1:12),2) Wgt=c(1:24) daily-cbind(Trial, Tanks, Day, Wgt) daily daily.sub-subset(daily, subset=Trial==2 Tanks==a4|Trial==2 Tanks==c4|Trial==2 Tanks==h4) daily.sub1-as.data.frame(daily.sub) attach(daily.sub1) daily.sub1 x11() plot(Day, Wgt) #plot(Day, Wgt, pch=c(2,19,21)[Tanks]) detach(daily.sub1) -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Using auto.key with two variable plots
On 1/30/2010 3:45 PM, Jonathan Greenberg wrote:
Rhelpers:
Having a problem solving this. I have an xyplot call that looks like
this:
print(xyplot(temp_species_EAM_Pred_Pop$x+temp_species_NULL_Pred_Pop$x~temp_species_EAM_Pred_Pop$Action,main=current_species,
xlab=Action,ylab=Predicted Pop,
xlim=c(xmin,xmax),ylim=c(ymin,ymax),
auto.key=list(corner=c(1,1
This is just a scatterplot with two response variables sharing the same
predictor variable (temp_species_EAM_Pred_Pop$Action). Right now, the
key has the words temp_species_EAM_Pred_Pop$x and
temp_species_NULL_Pred_Pop$x next to their symbols. I would like to
rename these in the key, say EAM and NULL -- using the group=
command doesn't work (since these aren't really different groups). What
is the right way to rename these variables in the key? Is using
auto.key the right approach?
Did you try setting the text parameter of the auto.key list?
Something like this:
print(xyplot(temp_species_EAM_Pred_Pop$x + temp_species_NULL_Pred_Pop$x
~ temp_species_EAM_Pred_Pop$Action,
main=current_species,
xlab=Action,ylab=Predicted Pop,
xlim=c(xmin,xmax),ylim=c(ymin,ymax),
auto.key=list(text=c('EAM','NULL'), corner=c(1,1
Thanks!
--j
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] extract R-squared and P-value from lm results
On 1/29/2010 9:04 AM, wenjun zheng wrote: Hi, R Users I find a problem in extracting the R-squared and P-value from the lm results described below (in Italic), *Residual standard error: 2.25 on 17 degrees of freedom* *Multiple R-squared: 0.001069, Adjusted R-squared: -0.05769 * *F-statistic: 0.01819 on 1 and 17 DF, p-value: 0.8943 * * * Any suggestions will be appreciated. Thanks. ?summary.lm In particular, see the 'Value' section which describes the components of the list returned when an lm() object is summarized. Notice the r.squared and coefficients components of the returned list. Wenjun [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Selective Plot Color
On 1/27/2010 10:00 AM, Richardson, Patrick wrote:
Is there a way (with a simple plot) to select all observations greater than a
certain value and plot them with a different color than the rest of the
observations in the plot? (i.e. for all observations greater than 10, I want
to plot them in red, but the rest of the observations remain black. Where can
I find how to do this?
Patrick
X - runif(30, min=5, max=15)
Y - rnorm(30)
plot(X, Y, col=ifelse(X 10, 'red', 'black'), pch=16)
The information transmitted is intended only for the per...{{dropped:10}}
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Extract R-squared from summary of lm
On 1/22/2010 10:10 AM, Trafim Vanishek wrote: Dear all, I cannot find to explicitly get the R-squared or adjusted R-squared from summary(lm()) fm1 - lm(Sepal.Length ~ Sepal.Width, data=iris) summary(fm1) str(summary(fm1)) summary(fm1)$r.squared ?str Thanks a lot! [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Estimation of S.E. based on bootstrapping (functions with two or more arguments)
On 1/21/2010 7:45 AM, Tomas Zelinsky wrote:
Hi all,
I need to estimate S.E. of a certain indicator. The function to compute the
value of indicator contains two arguments. Can anybody tell me how to do
it?
Example:
We have data:
a - c(1:10)
b - c(11:20)
data - data.frame(a, b)
Function to compute value of the indicator:
indicator - function(X, Y) sum(X)/(sum(Y)*2)
Next I need to do the bootstrapping and estimate mean value of indicator
and
its standard error.
If the function (indicator in my case) contained only one argument, there
would not be a problem, the code would look like:
resamples - lapply(1:1000, function(i) sample(data, replace = T))
r.indicator - sapply(resamples, indicator)
mean(r.indicator)
sqrt(var(r.indicator))
But in case of function with two arguments it doesn�t work. I tried to
do it
like:
resamples - lapply(1:1000, function(i) data[sample(1:nrow(data), replace =
TRUE),])
r.indicator - sapply(resamples, indicator)
but it didn't work.
Can anybody help?
How about using boot() in package boot? Using your example:
a - c(1:10)
b - c(11:20)
DF - data.frame(a,b)
library(boot)
boot(DF,
statistic = function(d, ind){sum(d$a[ind])/sum(d$b[ind]*2)},
R = 1000)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = DF, statistic = function(d, ind) {
sum(d$a[ind])/sum(d$b[ind] * 2)
}, R = 1000)
Bootstrap Statistics :
original biasstd. error
t1* 0.1774194 -0.001594390 0.01902264
Thanks a lot.
Tomas
__ Informacia od ESET NOD32 Antivirus, verzia databazy 4792
(20100121) __
Tuto spravu preveril ESET NOD32 Antivirus.
[1]http://www.eset.sk
References
1. http://www.eset.sk/
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] standardizing one variable by dividing each value by the mean - but within levels of a factor
On 1/20/2010 5:37 PM, Dimitri Liakhovitski wrote:
Hello!
I have a data frame with a factor and a numeric variable:
x-data.frame(factor=c(b,b,d,d,e,e),values=c(1,2,10,20,100,200))
For each level of factor - I would like to divide each value of
values by the mean of values that corresponds to the level of
factor
In other words, I would like to get a new variable that is equal to:
1/1.5
2/1.5
10/15
20/15
100/150
200/150
I realize I could do it through tapply starting with:
factor.level.means-tapply(x$values,x$factor,mean) ... etc.
But it seems clunky to me.
Is there a more elegant way of doing it?
with(x, ave(x=values, factor, FUN=function(x){x/mean(x)}))
[1] 0.667 1.333 0.667 1.333 0.667 1.333
?ave
Thanks a lot!
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Exchange NAs for mean
On 12/17/2009 9:31 AM, Joel Fürstenberg-Hägg wrote:
Hi all,
I'm have a matrix (X) with observations as rows and parameters as columns.
I'm trying to exchange all missing values in a column by the column mean
using the code below, but so far, nothing happens with the NAs... Can anyone
see where the problem is?
N-nrow(X) # Calculate number of rows = 108
p-ncol(X) # Calculate number of columns = 88
# Replace by columnwise mean
for (i in colnames(X)) # Do for all columns in the matrix
{
for (j in rownames(X)) # Go through all rows
{
if(is.na(X[j,i])) # Search for missing value in the given position
{
X[j,i]=mean(X[1:p, i]) # Change missing value to the mean of the
column
}
}
}
mymat - matrix(runif(50), ncol=5)
mymat[c(3,4,9),c(1,2,5)] - NA
mymat
[,1] [,2] [,3] [,4] [,5]
[1,] 0.05863667 0.23545794 0.25549983 0.96275422 0.1015316
[2,] 0.66107818 0.15846239 0.05112992 0.09081990 0.6097318
[3,] NA NA 0.86474629 0.73186676NA
[4,] NA NA 0.26226776 0.31987242NA
[5,] 0.78472732 0.09072242 0.24557669 0.57100857 0.1568413
[6,] 0.42431343 0.37551338 0.86085073 0.62782597 0.5090823
[7,] 0.44609972 0.90125504 0.52285650 0.41298482 0.3192929
[8,] 0.27676684 0.17200162 0.70648140 0.86983249 0.2035595
[9,] NA NA 0.34956846 0.07070571NA
[10,] 0.01482263 0.20074897 0.41553916 0.82367719 0.9674044
apply(mymat, 2, function(x){x - replace(x, is.na(x), mean(x,
na.rm=TRUE))})
[,1] [,2] [,3] [,4] [,5]
[1,] 0.05863667 0.23545794 0.25549983 0.96275422 0.1015316
[2,] 0.66107818 0.15846239 0.05112992 0.09081990 0.6097318
[3,] 0.38092069 0.30488025 0.86474629 0.73186676 0.4096348
[4,] 0.38092069 0.30488025 0.26226776 0.31987242 0.4096348
[5,] 0.78472732 0.09072242 0.24557669 0.57100857 0.1568413
[6,] 0.42431343 0.37551338 0.86085073 0.62782597 0.5090823
[7,] 0.44609972 0.90125504 0.52285650 0.41298482 0.3192929
[8,] 0.27676684 0.17200162 0.70648140 0.86983249 0.2035595
[9,] 0.38092069 0.30488025 0.34956846 0.07070571 0.4096348
[10,] 0.01482263 0.20074897 0.41553916 0.82367719 0.9674044
All the best,
Joel
_
Hitta hetaste singlarna på MSN Dejting!
http://dejting.se.msn.com/channel/index.aspx?trackingid=1002952
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] writing 'output.csv' file
On 12/4/2009 5:12 AM, Maithili Shiva wrote: Dear R helpers Suppose M - c(1:10) # length(M) = 10 N - c(25:50) # length(N) = 26 I wish to have an outut file giving M and N. So I have tried write.csv(data.frame(M, N), 'output.csv', row.names = FALSE) but I get the following error message Error in data.frame(M, N) : arguments imply differing number of rows: 10, 26 How do I modify my write.csv command to get my output in a single (csv) file irrespective of lengths. The first argument to write.csv() is preferably a matrix or data frame. If it is something else, write.csv() will try to make it a data frame. You may want to create the data frame like this: data.frame(M = c(M, rep(NA,length(N) - length(M))), N=N) M N 1 1 25 2 2 26 3 3 27 4 4 28 5 5 29 6 6 30 7 7 31 8 8 32 9 9 33 10 10 34 11 NA 35 12 NA 36 13 NA 37 14 NA 38 15 NA 39 16 NA 40 17 NA 41 18 NA 42 19 NA 43 20 NA 44 21 NA 45 22 NA 46 23 NA 47 24 NA 48 25 NA 49 26 NA 50 Plese Guide Thanks in advance Maithili The INTERNET now has a personality. YOURS! See your Yahoo! Homepage. [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Finding cases in one subset that are closet to another subset
On 12/2/2009 3:01 PM, Peter Flom wrote: Good afternoon Running R2.10.0 on Windows I have a data frame that includes (among much else) a factor (In_2006) and a continuous variable (math_3_4). I would like to find the 2 cases for In_2006 = 0 that are closest to each case where In_2006 = 1. My data looks like In_2006 math_3_4 0 55.1 1 51.6 1 18.1 1 26.6 1 14.1 1 9.6 1 48.9 1 12.9 0 63.0 0 51.8 etc. for several hundred rows. I would like a new data frame that has all the cases where In_2006 = 1, and those cases of In_2006 that are closest to those cases Hi Peter: How about using one of the various matching packages (MatchIt, optmatch, Matching)? For example, something like this: DF - data.frame(X = rbinom(200, 1, .1), Y = runif(200)) library(MatchIt) DF.match - matchit(X ~ Y, data=DF, method='optimal', ratio=2) DF[c(rownames(DF.match$match.matrix), c(DF.match$match.matrix)),] hope this helps, Chuck Thanks in advance Peter Peter L. Flom, PhD Statistical Consultant Website: www DOT peterflomconsulting DOT com Writing; http://www.associatedcontent.com/user/582880/peter_flom.html Twitter: @peterflom __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 z-standardize for subgroups?
On 11/29/2009 4:23 PM, John Kane wrote:
http://finzi.psych.upenn.edu/R/library/QuantPsyc/html/Make.Z.html
Make.Z in the QuantPsych package may already do it.
For a single variable, you could use ave() and scale() together like this:
with(iris, ave(Sepal.Width, Species, FUN = scale))
To scale more than one variable in a concise call, consider something
along these lines:
apply(iris[,1:4], 2, function(x){ave(x, iris$Species, FUN = scale)})
hope this helps,
Chuck Cleland
--- On Sun, 11/29/09, Karsten Wolf w...@uni-bremen.de wrote:
From: Karsten Wolf w...@uni-bremen.de
Subject: [R] How to z-standardize for subgroups?
To: r-help@r-project.org
Received: Sunday, November 29, 2009, 10:41 AM
Hi folks,
I have a dataframe df.vars with the follwing structure:
var1 var2 var3 group
Group is a factor.
Now I want to standardize the vars 1-3 (actually - there
are many more) by class, so I define
z.mean.sd - function(data){
return.values - (data -
mean(data)) / (sd(data))
return(return.values)
}
now I can call for each var
z.var1 - by(df.vars$var1, group, z.mean.sd)
which gives me the standardised data for each subgroup in a
list with the subgroups
z.var1 - unlist(z.var1)
then gives me the z-standardised data for var1 in one
vector. Great!
Now I would like to do this for the whole dataframe, but
probably I am not thinking vectorwise enough.
z.df.vars - by(df.vars, group, z.mean.sd)
does not work. I banged my head on other solutions trying
out sapply and tapply, but did not succeed. Do I need to
loop and put everything together by hand? But I want to keep
the columnnames in the vector…
-karsten
-
Karsten D. Wolf
Didactical Design of Interactive
Learning Environments
Universität Bremen - Fachbereich 12
web: http://www.ifeb.uni-bremen.de/wolf/
__
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.
__
Do You Yahoo!?
Tired of spam?
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] matching and extracting data
On 11/27/2009 10:25 AM, ram basnet wrote: Dear all, I have querry on how to extract the data by matching between two data set where one has the same elements multiple times? For example, I have two matrix X and Y. X [,1][,2] [,3] 1 A 5 P 2 B 6 P 3 C 7 P 4 D 5 Q 5 E 6 Q 6 F 7 Q 7 G 5 R 8 H 6 R 9 I 7 S 10 J 5 S 11 K6 T 12 L 7 T and Y [,1] 1 P 2 Q 3 R 4 S Now, I want to select and extract all the data of P, Q, R and S elements of column 3 of X matrix by matching with column 1 of Y matrix like below: [,1] [,2] [,3] 1 A 5 P 2 B 6 P 3 C 7 P 4 D 5 Q 5 E 6 Q 6 F 7 Q 7 G 5 R 8 H 6 R 9 I 7 S 10 J 5 S I guess, the answer might be simple but i am not getting way to figure out. And, i have to select subset from very huge data set. So, i need some kinds of automated procedure. If some one can help me, it will be great subset(X, X[,3] %in% Y[,1]) [,1] [,2] [,3] [1,] A 5 P [2,] B 6 P [3,] C 7 P [4,] D 5 Q [5,] E 6 Q [6,] F 7 Q [7,] G 5 R [8,] H 6 R [9,] I 7 S [10,] J 5 S Thanks in advance. Sincerely, Ram Kumar Basnet [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] interaction effects in a Linear model
On 11/23/2009 1:06 AM, ashok varma wrote: hello, i am fitting a Linear model using R. I have to fit the model considering all the interaction effects of order 1 of the independent variables. But I have 9 variables. So, it will be difficult for me to write all the 36 combinations in the model. Can anyone please help how to get this done in much smarter way?? This will include all main effects and two-way interactions: lm(Y ~ (A + B + C + D + E + F + G + H + I)^2, data=mydata) thanks a lot, Ashok Varma. [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] read.table (again)
On 11/4/2009 5:03 AM, Sybille Wendel (Udata) wrote: Dear R commnuity, Thanks a lot for your help. I want to read in tables, the problem is that the table is composed in a difficult way. In ariginal it looks like this: 669 736 842101610481029114711811166124312081128117611221026 9581024 992 685 720 829 925 995 96010241057116611501104106410711092 983 908 989 904 924 896 882 897 909 933 928 907 916 902 546 734 784 868 970 954 99210061012101610481045 9821057 989 914 956 960 948 947 949 939 916 950 562 598 771 827 941 922 905 877 860 862 931 952 954 977 960 901 978 975 970 950 973 953 958 931 912 558 593 725 786 884 866 838 797 815 802 809 822 853 833 863 852 903 9861015 9841019 993 955 932 918 874 518 580 600 764 834 804 814 824 849 831 939 757 769 790 809 892 89810251028104410371018 985 932 478 559 585 696 812 812 811 867 854 811 814 818 793 760 814 976 957104510551067106410501005 948 465 528 567 619 703 828 824 830 851 824 823 873 826 787 787 9791048105710621083108910841027 944 That is, that every 4 digits there is a new number, but when the number is 999, R thinks of course that the number consists of more than 4 digits. So, R can't read in the table. Is there a way I can tell R, that every 4 digits, a new number begins? ?read.fwf I treif different things with sep. data - read.table (RA19930101.txt,header=F,sep=) Thanks for your help again in this topic. Sybille -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Johnson-Neyman procedure (ANCOVA)
On 10/30/2009 9:20 PM, Stropharia wrote: Dear R users, Does anyone know of a package that implements the Johnson-Neyman procedure - for testing differences among groups when the regression slopes are heterogeneous (in an ANCOVA model)? I did not get any matches for functions when searching using R Site Search. If it has not been implemented in R, does anyone know of another statistical package that has this procedure? I found very old scripts available for Mathematica, but have had difficulty getting them to work. http://www.comm.ohio-state.edu/ahayes/SPSS%20programs/modprobe.htm Thanks in advance. best, Steve -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 is the difference between anova {stats} and aov?
On 10/27/2009 1:06 PM, Peng Yu wrote:
There are anova {stats} and aov in R. It seems that anova takes an
object returned by a model fitting function. But I don't see any
examples for anova. Can somebody give me a simple example on anova?
What is the difference between anova and aov?
http://finzi.psych.upenn.edu/Rhelp08/2009-January/184619.html
fm1 - aov(breaks ~ tension + wool, data=warpbreaks)
fm2 - aov(breaks ~ tension * wool, data=warpbreaks)
anova(fm1)
anova(fm2)
anova(fm1, fm2)
See ?anova.lm and ?anova.glm for more anova() examples.
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Bonferroni with unequal sample sizes
On 10/23/2009 1:30 PM, slrei...@vims.edu wrote: Hello- I have run an ANOVA on 4 treatments with unequal sample sizes (n=9,7,10 and 10). I want to determine where my sig. differences are between treatments using a Bonferroni test, and have run the code: pairwise.t.test(Wk16, Treatment, p.adf=bonf) I receive an error message stating that my arguments are of unequal length: Error in tapply(x, g, mean, na.rm = TRUE) : arguments must have same length Is there a way to run this test even with unequal sample sizes? Unequal sample sizes are not causing the reported error. The problem is that 'Wk16' and 'Treatment' do not have the same number of observations. If the group sizes are 9, 7, 10, and 10, then 'Wk16' and 'Treatment' should each contain 36 observations. Thanks. Stephanie Reiner Andrews Hall 410 Shipping: Rt. 1208 Greate Road Gloucester Point, VA 23062 work: 804-684-7869 cell: 443-286-4795 __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Simple question, I think
On 10/22/2009 12:37 PM, David Kaplan wrote: Greetings, I am recoding a dummy variable (coded 1,0) so that 0 = 2. I am using the line sciach$dummyba[sciach$ba==0] - 2 I notice that it creates a new column dummyba, with 0 coded as 2 but with 1's now coded as NA. Is there a simple way around this in the line I am using, or do I need to have an additional line sciach$dummyba[sciach$ba==1] - 1 You might do the recoding like this: with(sciach, ifelse(ba == 0, 2, ifelse(ba == 1, 1, NA))) or like this: sciach$ba * -1 + 2 Thanks in advance. David __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 average subgroups in a dataframe? (not sure how to apply aggregate(..))
On 10/21/2009 7:03 AM, Tony Breyal wrote:
Dear all,
Lets say I have the following data frame:
set.seed(1)
col1 - c(rep('happy',9), rep('sad', 9))
col2 - rep(c(rep('alpha', 3), rep('beta', 3), rep('gamma', 3)),2)
dates - as.Date(rep(c('2009-10-13', '2009-10-14', '2009-10-15'),6))
score=rnorm(18, 10, 3)
df1-data.frame(col1=col1, col2=col2, Date=dates, score=score)
col1 col2 Date score
1 happy alpha 2009-10-13 8.120639
2 happy alpha 2009-10-14 10.550930
3 happy alpha 2009-10-15 7.493114
4 happy beta 2009-10-13 14.785842
5 happy beta 2009-10-14 10.988523
6 happy beta 2009-10-15 7.538595
7 happy gamma 2009-10-13 11.462287
8 happy gamma 2009-10-14 12.214974
9 happy gamma 2009-10-15 11.727344
10 sad alpha 2009-10-13 9.083835
11 sad alpha 2009-10-14 14.535344
12 sad alpha 2009-10-15 11.169530
13 sad beta 2009-10-13 8.136278
14 sad beta 2009-10-14 3.355900
15 sad beta 2009-10-15 13.374793
16 sad gamma 2009-10-13 9.865199
17 sad gamma 2009-10-14 9.951429
18 sad gamma 2009-10-15 12.831509
Is it possible to get the following, whereby I am averaging the values
within each group of values in col2:
col1 col2 Date score Average
1 happy alpha 13/10/2009 8.120639 8.721561
2 happy alpha 14/10/2009 10.550930 8.721561
3 happy alpha 15/10/2009 7.493114 8.721561
4 happy beta 13/10/2009 14.785842 11.104320
5 happy beta 14/10/2009 10.988523 11.104320
6 happy beta 15/10/2009 7.538595 11.104320
7 happy gamma 13/10/2009 11.462287 11.801535
8 happy gamma 14/10/2009 12.214974 11.801535
9 happy gamma 15/10/2009 11.727344 11.801535
10 sad alpha 13/10/2009 9.083835 11.596236
11 sad alpha 14/10/2009 14.535344 11.596236
12 sad alpha 15/10/2009 11.169530 11.596236
13 sad beta 13/10/2009 8.136278 8.288990
14 sad beta 14/10/2009 3.355900 8.288990
15 sad beta 15/10/2009 13.374793 8.288990
16 sad gamma 13/10/2009 9.865199 10.882712
17 sad gamma 14/10/2009 9.951429 10.882712
18 sad gamma 15/10/2009 12.831509 10.882712
My feeling is that I should be using the ?aggregate is some fashion
but can't see how to do it. Or possibly there's another function i
should be using?
?ave
For example, try something like this:
transform(df1, Average = ave(score, col1, col2))
Thanks in advance,
Tony
O/S: Windows Vista Ultimate
sessionInfo()
R version 2.9.2 (2009-08-24)
i386-pc-mingw32
locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United Kingdom.
1252;LC_MONETARY=English_United Kingdom.
1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods
base
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Import SPSS file to R
On 10/19/2009 8:06 AM, Suman Kundu wrote:
Hello,
In R, How to read SPSS file and access the data item?
RSiteSearch('SPSS', restrict='function')
Thank you.
Regards,
Suman Kundu
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] modifying model coefficients
On 10/18/2009 11:26 PM, tdm wrote: I have build a model but want to then manipulate the coefficients in some way. I can extract the coefficients and do the changes I need, but how do I then put these new coefficients back in the model so I can use the predict function? my_model - lm(x ~ . , data=my_data) my_scores - predict(my_model, my_data) my_coeffs - coef(my_model) ## here we manipulate my_coeffs ## and then want to set the my_model ## coefficients to the new values so we ## predict using the new values my_model.coefs - my_coeffs ?? how is this done? ?? so that this will work with the new coefficients my_scores_new - predict(my_model, my_data) Any code snippets would be appreciated very much. Have you considered setting the coefficients to the values in my_coeffs using offset()? fm1 - lm(Sepal.Length ~ offset(0.90*Sepal.Width) + offset(0.50*Petal.Length) + offset(-0.50*Petal.Width), data = iris) summary(fm1) predict(fm1) all.equal(as.vector(predict(fm1)), c(as.matrix(cbind(1, iris[,2:4])) %*% c(1.8124, 0.9, 0.5, -0.5))) [1] TRUE ?offset -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 calculate average correlation coefficient of a correlation matrix ?
On 10/13/2009 10:13 AM, Amit Kumar wrote:
Hi! All,
I have large correlation matrix Cor. I wish to calculate average
correlation coefficient for this matrix.
Is there any function in R to do this?
Thanks in advance.
cormat - cor(iris[,1:4])
corlowtri - cormat[lower.tri(cormat)]
corlowtri
[1] -0.1175698 0.8717538 0.8179411 -0.4284401 -0.3661259 0.9628654
mean(corlowtri)
[1] 0.2900708
mean(abs(corlowtri))
[1] 0.594116
avgcor - function(x){mean(abs(x[lower.tri(x)]))}
avgcor(cor(iris[,1:4]))
[1] 0.594116
Amit
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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 column of frequency
On 9/29/2009 6:06 AM, abdul kudus wrote: Dear all, Given mypi mypi - c(0.1,0.2,0.2,0.1,0.3,0.4,0.4,0.4,0.4,0.2) I want to create myfreq as follows mypi myfreq 0.1 2 0.2 3 0.2 3 0.1 2 0.3 1 0.4 4 0.4 4 0.4 4 0.4 4 0.2 3 where myfreq is frequency of its corresponding observation. How to do that? mypi - c(0.1,0.2,0.2,0.1,0.3,0.4,0.4,0.4,0.4,0.2) ave(mypi, mypi, FUN=length) [1] 2 3 3 2 1 4 4 4 4 3 ?ave Thank you, Regards, A. Kudus Institute for Math Research Univ Putra Malaysia [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] xyplot lmline: error message.
On 9/28/2009 4:07 AM, Andrewjohnclose wrote:
Hi, I am trying to produce an xyplot with a regression line. The data should
be represented as log/log but when I fit the lmline I receive an error
message - the plot is fine without the log transformation, but the then the
plot is meaningless. I know it must be something simple, but I just can't
see it. Hope someone can help...
Thanks.
What do you expect to happen with the values of 0 in Pk? Have a look
at this:
log(rwpk$Pk)
[1] -1.108663 -1.698269 -2.551046 -2.796881
[5] -2.733368 -2.796881 -3.352407 -3.352407
[9] -4.074542 -3.244194 -4.074542 -4.342806
[13] -3.506558 -5.521461 -4.342806 -4.074542
[17] -4.342806 -3.506558 -4.342806 -4.342806
[21] -Inf -3.816713 -4.342806 -4.710531
[25] -4.342806 -Inf -Inf -4.074542
[29] -5.521461 -4.342806 -4.710531 -5.521461
[33] -5.521461 -5.521461 -5.521461 -5.521461
[37] -Inf -Inf -Inf -5.521461
[41] -Inf -Inf -5.521461 -Inf
[45] -Inf -Inf -5.521461 -Inf
[49] -Inf -Inf -Inf -Inf
[53] -Inf -5.521461 -Inf -Inf
[57] -Inf -Inf -Inf -Inf
[61] -Inf -Inf -Inf -5.521461
If you want to ignore the zeros, perhaps you could do this:
library(lattice)
rwpk - read.csv(file='http://www.nabble.com/file/p25641684/rwpk.csv')
xyplot(log(Pk) ~ log(k), data=subset(rwpk, Pk !=0), type=c('p','r'))
xyplot(log(Pk)~log(k),data=rwpk,cex=1,
panel=function(x,y){
panel.grid(h=-1, v=-1)
panel.xyplot(x,y,cex=1.0)
panel.lmline(x,y)
})
http://www.nabble.com/file/p25641684/rwpk.csv rwpk.csv
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Remove single entries
On 9/28/2009 12:03 PM, Raymond Danner wrote: Dear Community, I have a data set with two columns, bird number and mass. Individual birds were captured 1-13 times and weighed each time. I would like to remove those individuals that were captured only once, so that I can assess mass variability per bird. I¹ve tried many approaches with no success. Can anyone recommend a way to remove individuals that were captured only once? Thanks, Ray How about something like this? DF - data.frame(BIRD = rep(1:10, c(1,1,2,10,5,6,7,1,8,9)), MASS = rnorm(50,50,10)) DF$NOBS - with(DF, ave(MASS, BIRD, FUN=length)) subset(DF, NOBS 1) [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] cycles in a graphical model
On 9/16/2009 9:26 AM, shuva gupta wrote:
Hi,
Is there is any R package or existing codes in R to detect cycles in a
graphical model or DAG (Directed Acyclic Graph) ?
Thanks.
You might try this to search R packages:
library(sos)
findFn('Directed Acyclic Graph')
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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 merge the fitted values from a linear model?
On 9/1/2009 6:32 PM, kayj wrote: Hi All, I would like to run a linear model where the response is the duration of relief in days and the regressor is the drug dosage in mg. Then I would like compute the predicted values of the duration of relief from the model and merge it into the original data. I am not sure how the merge happens since if I have missing values in the data, R runs the resgression model but fitted values for some observations are not being calculated. Below is my R script Mydata-read.csv(file=”file1.csv”, header=T) Model-lm(y ~ x, data=Mydata) f-fitted(Model) Newdata-cbind(f , Mydata) Is Newdata merged correctly? Thanks for your help You might try something like this: DF - data.frame(Y = rnorm(20), X = sample(c(NA,0,1,2,3), size=20, replace=TRUE)) DF$f - fitted(lm(Y ~ X, data=DF, na.action=na.exclude)) DF Y X f 1 0.81371693 2 0.1116813 2 -0.36585221 0 -0.8160565 3 -1.07271855 0 -0.8160565 4 1.27182331 1 -0.3521876 5 -0.12492961 2 0.1116813 6 -1.84241736 0 -0.8160565 7 -0.28532869 1 -0.3521876 8 1.17361614 NA NA 9 0.88190221 3 0.5755502 10 0.92742858 1 -0.3521876 11 -1.18675102 0 -0.8160565 12 0.38076816 NA NA 13 -1.31518961 0 -0.8160565 14 -1.07973072 1 -0.3521876 15 0.00431749 3 0.5755502 16 0.49820163 3 0.5755502 17 -0.21377954 1 -0.3521876 18 -1.03107537 2 0.1116813 19 -1.23459162 0 -0.8160565 20 -0.05666561 0 -0.8160565 -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Lattice xyplot: modify line width of plot lines
On 8/24/2009 4:47 AM, ukoe...@med.uni-marburg.de wrote: # Hi all, # I want to increase the line width of the plotted lines # in a xy-lattice plot. My own attempts were all in vain. # Without the group option the line width is modified - # with the option it is funnily enough not. # Please have a look at my syntax. # # Many thanks in advance # Udo You need to change the superpose.line setting: xyplot(Choline ~ time, groups=BMI, data=data, type=l, scales=list(relation=free), auto.key=list(title=BMI Group, border=FALSE, lines=TRUE, points=FALSE), xlab=c(Point in Time), ylab=c(Concentration of Choline), par.settings = list(superpose.line = list(lwd=3))) library(lattice) data - data.frame(cbind(1:2,c(1,1,2,2), c(0.5,0.9,1.0,1.8))) names(data) - c(BMI,time,Choline) data$BMI - factor(data$BMI) levels(data$BMI) - c(=17.5,17.5) data$time - factor(data$time) levels(data$time) - c(Admission,Discharge) #Show names of settings names(trellis.par.get()) #Try to set the line width of the two plotted colored lines line-trellis.par.get(plot.line) line line$lwd=3 trellis.par.set(plot.line, line) line #Without group option: Line width is changed xyplot(Choline ~ time, data=data, type=l, scales=list(relation=free), auto.key=list(title=BMI Group, border=FALSE), xlab=c(Point in Time), ylab=c(Concentration of Choline)) #With group option: Line width is not changed xyplot(Choline ~ time, group=BMI, data=data, type=l, scales=list(relation=free), auto.key=list(title=BMI Group, border=FALSE), xlab=c(Point in Time), ylab=c(Concentration of Choline)) __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] vector replacement 1/0 to P/A
On 8/17/2009 10:22 AM, Lana Schaffer wrote:
Hi,
Can someone suggest an efficient way to substitute a vector/matrix
which contains 1's and 0's to P's and A's (resp.)?
Thanks,
Lana
Here is one approach:
mymat - matrix(rbinom(15, 1, .5), ncol=3)
mymat
[,1] [,2] [,3]
[1,]100
[2,]001
[3,]101
[4,]010
[5,]110
mymat[] - sapply(mymat, function(x){ifelse(x == 1, 'P', ifelse(x == 0,
'A', NA))})
mymat
[,1] [,2] [,3]
[1,] P A A
[2,] A A P
[3,] P A P
[4,] A P A
[5,] P P A
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] post hoc test after lme
On 8/14/2009 10:06 AM, Casimir wrote: Hi! I am quiet new with R and I have some problems to perform a posthoc test with an lme model. My model is the following: lme1-lme(eexp~meal+time, random=~1|id,na.action=na.omit) and then i try to get a post hoc test: summary(glht(lme1,linfct=mcp(meal=Tukey))) but I get a warning message: Erreur dans as.vector(x, mode) : argument 'mode' incorrect Thank you for your help Guillaume Use the data argument in the call to lme(). For example, this works: library(multcomp) library(nlme) Orthodont$AGECAT - as.factor(Orthodont$age) fm1 - lme(distance ~ AGECAT + Sex, data = Orthodont, random = ~ 1 | Subject) summary(glht(fm1, linfct=mcp(AGECAT = Tukey))) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: lme.formula(fixed = distance ~ AGECAT + Sex, data = Orthodont, random = ~1 | Subject) Linear Hypotheses: Estimate Std. Error z value Pr(|z|) 10 - 8 == 00.9815 0.3924 2.501 0.05984 . 12 - 8 == 02.4630 0.3924 6.277 0.001 *** 14 - 8 == 03.9074 0.3924 9.958 0.001 *** 12 - 10 == 0 1.4815 0.3924 3.776 0.001 *** 14 - 10 == 0 2.9259 0.3924 7.457 0.001 *** 14 - 12 == 0 1. 0.3924 3.681 0.00123 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method) But this does not work and returns the error you got: attach(Orthodont) fm2 - lme(distance ~ AGECAT + Sex, random = ~ 1 | Subject) summary(glht(fm2, linfct=mcp(AGECAT = Tukey))) Error in as.vector(x, mode) : invalid 'mode' argument -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] lme funcion in R
On 8/3/2009 1:15 PM, Hongwei Dong wrote: Hi, R users, I'm using the lme function in R to estimate a 2 level mixed effects model, in which the size of the subject groups are different. It turned out that It takes forever for R to converge. I also tried the same thing in SPSS and SPSS can give the results out within 20 minutes. Anyone can give me some advice on the lme function in R, especially why R does not converge? Thanks. Harry: You are much more likely to get helpful advice if you include the code you used to attempt to fit the model and a brief description of the data. For example, something along these lines but for your data and model: library(nlme) fm2 - lme(distance ~ age + Sex, data = Orthodont, random = ~ 1) str(Orthodont) Classes ‘nfnGroupedData’, ‘nfGroupedData’, ‘groupedData’ and 'data.frame': 108 obs. of 4 variables: $ distance: num 26 25 29 31 21.5 22.5 23 26.5 23 22.5 ... $ age : num 8 10 12 14 8 10 12 14 8 10 ... $ Subject : Ord.factor w/ 27 levels M16M05M02..: 15 15 15 15 3 3 3 3 7 7 ... $ Sex : Factor w/ 2 levels Male,Female: 1 1 1 1 1 1 1 1 1 1 ... - attr(*, outer)=Class 'formula' length 2 ~Sex .. ..- attr(*, .Environment)=environment: R_GlobalEnv - attr(*, formula)=Class 'formula' length 3 distance ~ age | Subject .. ..- attr(*, .Environment)=environment: R_GlobalEnv - attr(*, labels)=List of 2 ..$ x: chr Age ..$ y: chr Distance from pituitary to pterygomaxillary fissure - attr(*, units)=List of 2 ..$ x: chr (yr) ..$ y: chr (mm) - attr(*, FUN)=function (x) ..- attr(*, source)= chr function (x) max(x, na.rm = TRUE) - attr(*, order.groups)= logi TRUE hope this helps, Chuck Harry [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] edit.row.names taking row names from the edited dataframe
On 7/30/2009 2:46 PM, Ross Culloch wrote:
Hi all,
I am struggling to work out how to use the rownames from an edited dataframe
rather than the row names from the original dataframe. In my data set i'm
trying to extract several rows of data on particular individuals, i don't
doubt i'm using the long way round but what i have in the way of a script is
this:
##selecting the IDs from the dataframe individually
A1-mumpup[ID==A1,]
B1-mumpup[ID==B1,]
B2-mumpup[ID==B2,]
B3-mumpup[ID==B3,]
B4-mumpup[ID==B4,]
B6-mumpup[ID==B6,]
B7-mumpup[ID==B7,]
B8-mumpup[ID==B8,]
B9-mumpup[ID==B9,]
B13-mumpup[ID==B13,]
C1-mumpup[ID==C1,]
G1-mumpup[ID==G1,]
data-rbind(A1,B1,B2,B3,B4,B6,B7,B8,B9,B13,C1,G1)
It works fine to a certain extent, the only problem being that i get are all
the IDs in the original dataframe so if i use
summary(data$ID)
I get:
A1 B1 B10 B13 B2 B3 B4 B6 B7 B8 B9 C1 G1 G3 H2 H9 J1 J2 J3
K1
354 354 0 354 354 354 354 354 354 354 354 246 210 0 0 0 0 0 0
0
So it does take the data out of the dataframe but it keeps the IDs for some
reason.
I have looked at the edit.data.frame help and i understand why it is
happening, it is taking the rownames from the original dataframe and not the
edit and it seems i should use edit.row.names=T in my script, but i can't
get that to work.
Does anyone have any suggestions at all? Any help is much appreciated,
If I understand, you may want something like this:
myIDs - c('A1','B1','B2','B3','B4','B6','B7','B8','B9','B13','C1','G1')
subset(mumpup, ID %in% myIDs)
Best wishes,
Ross
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Welch Anova ?
On 7/24/2009 2:51 AM, Hardi wrote: Hi, I need to do factor analysis with non-constant variance. Is there a package that contains Welch ANOVA ? Thanks, ?oneway.test __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] metafor
On 7/24/2009 11:40 AM, Frank Pearson wrote: I had found the author's (Wolfgang Viechtbauer) earlier meta-analytic code in R, MiMa, useful. so I have been exploring metafor using an example dataset from MiMa. metafor provides a lot more. However, MiMa provided parameter estimates, standard errors, z values, etc. for individual moderators in the meta-analysis, but I don't see how to obtain these from metafor. Have you any help about how to get this? Frank Hi Frank: Which function in the metafor package are you using to fit the meta-regression model? The following works for me to reproduce one of the analyses in the MiMa Tutorial (http://www.wvbauer.com/downloads/mima_tutorial.pdf): library(metafor) f - STUDY N1 N2 YI VI MINUTES TRAINED MEANAGE 1 30 30 0.444 0.068 30 0 28 2 46 39 -0.495 0.049 10 0 42 3 15 15 0.195 0.134 20 1 31 4 10 10 0.546 0.207 40 1 39 5 12 12 0.840 0.181 20 1 17 6 10 10 0.105 0.200 30 1 51 7 26 24 0.472 0.082 15 1 26 8 18 14 -0.205 0.128 10 0 64 9 12 12 1.284 0.201 45 1 48 10 12 12 0.068 0.167 30 1 40 11 15 15 0.234 0.134 30 1 52 12 12 12 0.811 0.180 30 1 33 13 15 15 0.204 0.134 30 1 20 14 18 18 1.271 0.134 60 1 27 15 15 15 1.090 0.153 45 1 52 16 43 35 -0.059 0.052 10 1 61 madf - read.table(textConnection(f), sep= , header=TRUE) rma(YI, VI, mods=cbind(MINUTES, TRAINED, MEANAGE), data=madf) Mixed-Effects Model (k = 16; tau^2 estimator: REML) tau^2 (estimate of residual amount of heterogeneity): 0 (SE = 0.0472) tau (sqrt of the estimate of residual heterogeneity): 0 Test for Residual Heterogeneity: QE(df = 12) = 9.1238, p-val = 0.6923 Test of Moderators (coefficient(s) 2,3,4): QM(df = 3) = 28.8348, p-val = 0 Model Results: estimate se zvalpvalci.lb ci.ub intrcpt -0.2956 0.3488 -0.8473 0.3968 -0.9792 0.3881 MINUTES0.0250 0.0067 3.7149 0.0002 0.0118 0.0381 *** TRAINED0.3011 0.1953 1.5415 0.1232 -0.0817 0.6839 MEANAGE -0.0060 0.0063 -0.9518 0.3412 -0.0182 0.0063 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 sessionInfo() R version 2.9.1 (2009-06-26) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] metafor_0.5-0 hope this helps, Chuck [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Calculate weighted mean for each group
On 7/23/2009 5:18 PM, Alexis Maluendas wrote:
Hi R experts,
I need know how calculate a weighted mean by group in a data frame. I have
tried with aggragate() function:
data.frame(x=c(15,12,3,10,10),g=c(1,1,1,2,2,3,3),w=c(2,3,1,5,5,2,5)) - d
aggregate(d$x,by=list(d$g),weighted.mean,w=d$w)
Generating the following error:
Error en FUN(X[[1L]], ...) : 'x' and 'w' must have the same length
DF - data.frame(x=c(15,12, 3,10,10,14,12),
g=c( 1, 1, 1, 2, 2, 3, 3),
w=c( 2, 3, 1, 5, 5, 2, 5))
sapply(split(DF, DF$g), function(x){weighted.mean(x$x, x$w)})
123
11.5 10.0 12.57143
Thanks in advance
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Find multiple elements in a vector
On 7/22/2009 3:32 PM, Michael Knudsen wrote: Hi, Given a vector, say x=sample(0:9,10) x [1] 0 6 3 5 1 9 7 4 8 2 I can find the location of an element by which(x==2) [1] 10 but what if I want to find the location of more than one number? I could do c(which(x==2),which(x==3)) but isn't there something more streamlined? My first guess was y=c(2,3) which(x==y) integer(0) which doesn't work. I haven't found any clue in the R manual. How about this? which(x %in% c(2,3)) Thanks! -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 transform m/d/yyyy to yyyymmdd?
On 7/21/2009 1:16 PM, liujb wrote: Hello, I have a set of data that has a Date column looks like this: 12/9/2007 12/16/2007 1/1/2008 1/3/2008 1/12/2008 etc. I'd like the date to look something like the follow (so that I could sort by date easily). 20071209 20071216 20080101 20080103 20080112 How to do it? Thank you very much Julia dates - c(2/27/1992, 2/27/1992, 1/14/1992, 2/28/1992, 2/1/1992) as.character(as.Date(dates, %m/%d/%Y), %Y%m%d) [1] 19920227 19920227 19920114 19920228 19920201 ?as.Date -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] calculating median with a condition
On 7/20/2009 2:59 PM, Manisha Brahmachary wrote:
Hello,
I am trying to calculate the median of numbers across each row for the data
shown below , with the condition that if the number is negative, that it
should be ignored and the median should be taken of only the positive
numbers.
For eg: data is in Column A,B,C. Column D and E demonstrates what I want to
get as answer
A
B
C
Median
median value
-13.6688115
-32.50914055
-50.54011892
all negative, so ignore
NA
NA
-53.65656268
42.58599666
median C
42.58599666
33.30683089
18.93765489
-25.17024229
median A,B
26.12224289
The R script I have written is below( which doesnot do the job properly)
median.value- matrix(nrow=nrow(data),ncol=1)
for(k in 1:nrow(data)){
median.value[k]-median(data[which(data[k,]0)])}
Can someone suggest me the correct R script to do what I have explained
above.
X - as.data.frame(matrix(rnorm(100), ncol=10))
apply(X, 1, function(x){median(x[x 0])})
[1] 0.2297943 0.6476565 0.4699609 0.8744830
[5] 1.0242502 0.7800703 0.6648436 0.2930191
[9] 0.6001506 1.0767194
Thanks
Manisha
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] duplicate data points on a line graph
On 7/15/2009 9:56 PM, Carl Witthoft wrote:
If you want to take the second approach, it can be relatively easily
generalized by calculating the cex values based on the count of ordered
pairs in the original dataset.
Here's a data set:
xy
x y
[1,] 1 4
[2,] 1 5
[3,] 2 3
[4,] 3 3
[5,] 4 5
[6,] 5 2
[7,] 1 4
[8,] 2 3
Here's the same set fully sorted:
xy[order(x,y),]-xyord
x y
[1,] 1 4
[2,] 1 4
[3,] 1 5
[4,] 2 3
[5,] 2 3
[6,] 3 3
[7,] 4 5
[8,] 5 2
There's gotta be some very simple way to create a series of values for
cex but I'm missing it, other than a loop like
cexvec-rep(1,8)
for i in 2:8 {
if (xyord[i,1]==xyord[i-1,1] xyord[i,2]== xyord[i-1,2] ) {
cexvec[i]-cexvec[i-1]+1
}
}
How about using ave() like this:
x - sample(0:4, 60, replace=TRUE)
y - sample(0:4, 60, replace=TRUE)
xy - data.frame(x, y)
xy$freq - ave(xy$x, x, y, FUN=length)
with(xy, plot(x, y, cex=freq))
You get the idea, sort of :-)
Carl
On 7/15/2009 2:19 PM, NDC/jshipman wrote:
Hi,
I am new to R plot. I am trying to increase the data point
observation when duplicate data points exist
xy
110
110
23
45
9 8
in the about example 1, 10 would be displayed larger than the other
data points. Could someone give me some assistance with this problem
A couple of simple approaches:
x - c(1,1,2,4,9)
y - c(10,10,3,5,8)
plot(jitter(x), jitter(y))
plot(x, y, cex=c(2,2,1,1,1))
757-864-7114
LARC/J.L.Shipman/jshipman
Jeffery.L.Shipman at nasa.gov
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] duplicate data points on a line graph
On 7/15/2009 2:19 PM, NDC/jshipman wrote: Hi, I am new to R plot. I am trying to increase the data point observation when duplicate data points exist xy 110 110 23 45 9 8 in the about example 1, 10 would be displayed larger than the other data points. Could someone give me some assistance with this problem A couple of simple approaches: x - c(1,1,2,4,9) y - c(10,10,3,5,8) plot(jitter(x), jitter(y)) plot(x, y, cex=c(2,2,1,1,1)) 757-864-7114 LARC/J.L.Shipman/jshipman jeffery.l.ship...@nasa.gov __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] assessing data variation
On 7/10/2009 5:21 AM, e-letter wrote: I have data like so: time datum 3012 6024 9037 120 41 150 8 In addition to standard deviation, I want to measure the average of the differences in data for each time interval, i.e. average of 24-12, 37-24, 41-37, 8-41. Is there a statistical term for this task? Which package should I use please? I don't know a term for it, but you might try something like this: mydf - data.frame(time = c(30,60,90,120,150), datum = c(12,24,37,41,8)) diff(mydf$datum) [1] 12 13 4 -33 mean(abs(diff(mydf$datum))) [1] 15.5 ?diff rh...@conference.jabber.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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] variable driven summary of one column
On 6/25/2009 5:44 AM, Anne Skoeries wrote: Hello, how can I get a variable driven summary of one column of my data.frame? Usually I would do summary(data$columnname) to get a summary of column named columnname of my data.frame named data. In my case the columnname is not static but can be set dynamically. So I save the chosen columname in something like variable - columnname but how can I now get the summary of the specified column? summary(data$get(variable)) doesn't work. summary(paste(data$, variable, sep=) doesn't work either! and if I try summary(data[get(variable)] it gives me back a different result since the data isn't a factor anymore but a list. vname - Species summary(subset(iris, select=vname)) Species setosa:50 versicolor:50 virginica :50 vname - Sepal.Width summary(subset(iris, select=vname)) Sepal.Width Min. :2.000 1st Qu.:2.800 Median :3.000 Mean :3.057 3rd Qu.:3.300 Max. :4.400 Thanks for the help, Anne __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] distinguish regression lines in grouped, black and white lattice xyplot
On 6/24/2009 3:28 PM, Katharina May wrote:
Hi,
I've got the following problem which I cannot think of a solution right now:
if got a lattice xyplot in black and white and a grouping variable
with many (more than 8
values) and I plot it as regression lines (type=r), just like this
one (not reproducable but that's
I guess not the point here):
xyplot(log(AGWB) ~ log(BM_roots), data=sub_agwb_data, groups=species,
type=r, lty=c(1:6),panel=allo.panel.5)
The problem is that I've got 26 different values for the grouping
variable species and only 6 default values for the line type
lty (and according to the par {graphics} help page customizable to up
to 8 different line types).
Does anybody have any idea how these 26 different lines can be made
distinguishable from each other without the use
of colors?
If you need to distinguish individual regression lines, I would
consider 26 panels rather than attempting one panel with 26 regression
lines each of a different line type. Something like this:
xyplot(log(AGWB) ~ log(BM_roots) | species, data=sub_agwb_data,
type=r, panel=allo.panel.5)
Thanks,
Katharina
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] vector and NA
On 6/23/2009 5:41 AM, Alfredo Alessandrini wrote: Hi, I've a vector like this: -- inc[,5] [1]NANANANANANANA [8]NANANANANANANA [15]NANANANANANANA [22]NANANANANANANA [29]NANANANANANANA [36]NANANANANANANA [43]NANANANANANANA [50]NANANANANANANA [57]NANANANANANANA [64]NANANANANANANA [71]NANANANANANANA [78]NANANANA 13.095503 10.140119 7.989186 [85] 8.711888 7.201234 13.029250 14.430755 8.662832 8.810785 14.421302 [92] 7.614985 7.548091 9.843389 14.977402 20.875255 7.787543 2.005056 [99] 4.016916 3.601773 4.140390 7.241999 13.280794 18.038902 18.762169 [106] 4.280065 5.942021 6.292010 11.866446 19.450442 11.942362 6.224328 [113] 3.176050 5.456117 2.733487 3.992823 13.633171 19.514301 25.085256 [120] 5.640089 5.890486 12.421150 18.821420 22.478664 11.503805 7.051254 [127] 7.560921 12.000394 20.464875 16.147598 13.746290 9.416060 35.848221 [134] 36.739481 23.516759 7.317599 3.928247 10.371437 11.202935 12.574649 [141] 6.906980 9.191260 7.080267 2.810271 5.494705 10.617141 14.578020 [148] 10.981610 7.343975 2.179511 2.726651 10.794842 9.872493 19.842701 [155] 10.525064 16.134541 29.283385 18.352996 9.216318 6.253805 2.704267 [162] 4.274514 3.138237 12.296835 20.982433 13.001104 2.606328 3.333271 [169] 5.514425 2.179244 5.381514 6.848380 3.794428 5.114591 4.975830 [176] 3.809948 10.131608 14.145913 --- How can I extract a vector without the NA value? na.omit(inc[,5]) ?na.omit Regards, Alfredo __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] SAS-like method of recoding variables?
On 6/22/2009 2:27 PM, Mark Na wrote: Dear R-helpers, I am helping a SAS user run some analyses in R that she cannot do in SAS and she is complaining about R's peculiar (to her!) way of recoding variables. In particular, she is wondering if there is an R package that allows this kind of SAS recoding: IF TYPE='TRUCK' and count=12 THEN VEHICLES=TRUCK+((CAR+BIKE)/2.2); Thanks for any help or suggestions you might be able to provide! If the variables are in a data frame called mydf, she might do something like this: mydf$VEHICLE - with(mydf, ifelse(TYPE=='TRUCK' count==12, TRUCK+((CAR+BIKE)/2.2), NA)) or mydf - transform(mydf, VEHICLE = ifelse(TYPE=='TRUCK' count==12, TRUCK+((CAR+BIKE)/2.2), NA)) Mark Na __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] resampling the entire row
On 6/19/2009 10:59 AM, Seunghee Baek wrote: Hi, For bootstrapping method, I would like to resample the entire row instead of one column. What should I do? iris[sample(x=nrow(iris), replace=TRUE),] But I would look at the boot package or other packages related to bootstrapping. Thanks, Becky __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Hi
On 6/15/2009 5:42 PM, Oscar Bayona wrote:
Hi I have a simple question. I want to run a n times a simple linear
regession and save beta in a matrix but I´m not able.
Imagine:
Data.txt is a 10*5 file and want to run 4 different stimations always
regressing first column on the rest.
So I try this:
First I run Data on memory
This is my function
mrp - function(){
mr-matrix(0,4,1)
for(i in 1:4)
r(i)=lm(dat(,i+1)~dat(,1)
mr[i] - coefficients(r(i)))
}
I execute mrp usin source file choose option but nothing happens
Where I´m wrong?
It's hard to tell exactly what you want, but does this help?
mr - lm(as.matrix(cbind(dat[,2:ncol(dat)])) ~ dat[,1])
summary(mr)
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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 fix my nested conditional IF ELSE code?
On 6/14/2009 6:18 PM, Mark Na wrote:
Hi,
I've been struggling most of the morning with an IF ELSE problem, and I
wonder if someone might be able to sort me out.
Here's what I need to do (dummy example, my data are more complicated):
If type = A or B or C
and status = a then count = 1
and status = b then count = 2
and status = c then count = 3
Else if type = D or E or F
and status = a then count = 9
and status = b then count = 8
and status = c then count = 7
End
Seems simple when I write it like that, but the R code is escaping me.
mydf - data.frame(type = sample(LETTERS[1:6], 40, replace=TRUE),
status = sample(letters[1:3], 40, replace=TRUE))
mydf$count - with(mydf,
ifelse(type %in% c('A','B','C') status == 'a', 1,
ifelse(type %in% c('A','B','C') status == 'b', 2,
ifelse(type %in% c('A','B','C') status == 'c', 3,
ifelse(type %in% c('D','E','F') status == 'a', 9,
ifelse(type %in% c('D','E','F') status == 'b', 8,
ifelse(type %in% c('D','E','F') status == 'c', 7,
NA)))
mydf
type status count
1 C c 3
2 F b 8
3 B b 2
4 A a 1
5 C b 2
6 D c 7
7 C a 1
8 E c 7
9 F a 9
10E b 8
11F a 9
12D a 9
13A c 3
14B b 2
15D b 8
16C c 3
17E c 7
18A c 3
19B a 1
20E a 9
21E b 8
22E a 9
23C b 2
24A b 2
25F b 8
26D a 9
27B c 3
28E b 8
29E c 7
30F b 8
31B a 1
32F c 7
33F b 8
34E c 7
35C b 2
36C a 1
37F b 8
38D c 7
39F a 9
40D b 8
Thanks!
Mark Na
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Replacing 0s with NA
On 6/12/2009 4:55 AM, Christine Griffiths wrote:
Hello
I have a dataset in which I would like to replace 0s with NAs. There is
a lot of information on how to replace NAs with 0, but I have struggled
to find anything with regards to doing the reverse. Any recommendations
would be great.
X - as.data.frame(matrix(sample(0:9, 100, replace=TRUE), ncol=10))
X
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 7 5 8 9 6 6 4 8 5 8
2 6 2 9 6 8 5 9 1 1 3
3 4 1 5 8 9 5 3 2 1 4
4 4 8 7 7 4 1 1 4 9 8
5 9 2 5 8 4 8 4 8 6 0
6 3 4 2 8 2 0 6 4 8 5
7 3 5 0 2 7 7 9 9 3 1
8 7 3 3 4 8 3 9 2 7 1
9 4 7 9 1 5 4 8 2 1 9
10 7 7 6 1 0 9 0 5 7 0
X[] - lapply(X, function(x){replace(x, x == 0, NA)})
X
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 7 5 8 9 6 6 4 8 5 8
2 6 2 9 6 8 5 9 1 1 3
3 4 1 5 8 9 5 3 2 1 4
4 4 8 7 7 4 1 1 4 9 8
5 9 2 5 8 4 8 4 8 6 NA
6 3 4 2 8 2 NA 6 4 8 5
7 3 5 NA 2 7 7 9 9 3 1
8 7 3 3 4 8 3 9 2 7 1
9 4 7 9 1 5 4 8 2 1 9
10 7 7 6 1 NA 9 NA 5 7 NA
Cheers
Christine
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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 with if statements
On 6/9/2009 12:17 PM, Amit Patel wrote:
Hi
I am trying to create a column in a data frame which gives a sigificane score
from 0-7. It should read values from 7 different colums and add 1 to the
counter if the value is =0.05. I get an error message saying
Error in if (ALLRESULTS[i, 16] = 0.05) significance_count =
significance_count + :
missing value where TRUE/FALSE needed
The script is included below
it works if i convert the NA values to zero but this is not appropriate as it
includes the zero as significant.
ANY SUGGESTIONS
significance.count - rowSums(ALLRESULTS[,16:22] = .05, na.rm=TRUE)
?rowSums
#SCRIPT STARTS
for (i in 1:length(ALLRESULTS[,1])) {
significance_count = 0
if (ALLRESULTS[i,16] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
if (ALLRESULTS[i,17] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
if (ALLRESULTS[i,18] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
if (ALLRESULTS[i,19] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
if (ALLRESULTS[i,20] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
if (ALLRESULTS[i,21] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
if (ALLRESULTS[i,22] = 0.05 ) significance_count = significance_count +1
else significance_count = significance_count
ALLRESULTS[i,23] - significance_count}
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] if else
On 6/8/2009 1:48 PM, Cecilia Carmo wrote: Hi R-helpers! I have the following dataframe: firm-c(rep(1:3,4)) year-c(rep(2001:2003,4)) X1-rep(c(10,NA),6) X2-rep(c(5,NA,2),4) data-data.frame(firm, year,X1,X2) data So I want to obtain the same dataframe with a variable X3 that is: X1, if X2=NA X2, if X1=NA X1+X2 if X1 and X2 are not NA So my final data is X3-c(15,NA,12,5,10,2,15,NA,12,5,10,2) finaldata-data.frame(firm, year,X1,X2,X3) library(fortunes) fortune(dog) firm - c(rep(1:3, 4)) year - c(rep(2001:2003, 4)) X1 - rep(c(10, NA), 6) X2 - rep(c(5, NA, 2), 4) mydata - data.frame(firm, year, X1, X2) mydata$X3 - with(mydata, ifelse( is.na(X1) !is.na(X2), X2, ifelse(!is.na(X1) is.na(X2), X1, ifelse(!is.na(X1) !is.na(X2), X1 + X2, NA mydata firm year X1 X2 X3 1 1 2001 10 5 15 2 2 2002 NA NA NA 3 3 2003 10 2 12 4 1 2001 NA 5 5 5 2 2002 10 NA 10 6 3 2003 NA 2 2 7 1 2001 10 5 15 8 2 2002 NA NA NA 9 3 2003 10 2 12 101 2001 NA 5 5 112 2002 10 NA 10 123 2003 NA 2 2 I've tried this finaldata-ifelse(data$X1==NA,ifelse(data$X2==NA,NA,X2),ifelse(data$varvendas==NA,X1,X1+X2)) But I got just NA in X3. Anyone could help me with this? Thanks in advance, Cecília (Universidade de Aveiro - Portugal) __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] creating list with 200 identical objects
On 6/1/2009 6:08 AM, Rainer M Krug wrote:
Hi
I am doing an simulation, and I a large proportion of the simulation
time is taken up by memory allocations.
I am creating an object, and storing it in a list of those objects.
essentially:
x - list()
for (t in 1:500) {
x[1] - new(track)
}
I would like to initialize in one go, to avoid the continuous
reallocation of memory when a new track is added, and fill it wit
the object created by new(track).
How can I do this?
Does this help?
x - vector(list, 500)
for(i in 1:500){x[[i]] - runif(30)}
thanks
Rainer
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Sample size calculation proportions with EpiR: Discrepancy to other calculators
On 5/26/2009 2:53 AM, Karl Knoblick wrote: Hallo! I have done a sample size calculation for proportions with EpiR. The input is: treatment group rate p=0.65 control group rate p=0.50 significance level 0.95 power 0.80 two-sided ration group 1 and 2: 1.0 I have done this in the following way: library(epiR) epi.studysize(treat = 0.65, control = 0.5, n = NA, sigma = NA, power = 0.80, r = 1, conf.level = 0.95, sided.test = 2, method = proportions) Result: $n [1] 82 PASS 2002 and NQuery give both 170 subjects per group without continuity correction. With continuity correction 183 per group. Looking at http://statpages.org/proppowr.html I get 182 subjects per group (with continuity correction, I admit). What am I doing wrong? Can anybody explain this? epi.studysize(treat = .65, control = .50, n = NA, sigma = NA, power = 0.80, r = 1, conf.level = 0.95, sided.test = 2, method = cohort) gives the same sample size as PASS 2002 and NQuery (170 per group). Best wishes Karl __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Forcing a variableinto a model using stepAIC
On 5/22/2009 9:58 AM, Laura Bonnett wrote: Dear All, I am attempting to use forward and/or backward selection to determine the best model for the variables I have. Unfortunately, because I am dealing with patients and every patient is receiving treatment I need to force the variable for treatment into the model. Is there a way to do this using R? (Additionally, the model is stratified by randomisation period). I know that SAS can be used to do this but my SAS coding is poor and consequently it would be easier for me to use R, especially given the fractional polynomial transformations! Currently the model is as follows (without treatment). coxfita=coxph(Surv(rem.Remtime,rem.Rcens)~sind(nearma)+fsh(nearma)+fdr(nearma)+th1(nearma)+th2(nearma)+fp(cage)+fp(fint)+fp(tsb)+strata(rpa),data=nearma) Thank you for your help, Laura See the scope argument to stepAIC in the MASS package. You can specify a formula in the 'lower' component of scope which includes the treatment variable. That will force the treatment variable to remain in every model examined in the stepwise search. __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] SAS PROC SORT nodupkey
On 5/12/2009 8:29 AM, amor Gandhi wrote: Hi, I have the following data and I would like to delete douple names, it is almost similar to SAS PROC SORT nodupkey! Is there any function in R does this? x1 - rnorm(11,5,1) x2 - runif(11,0,1) nam -paste(A, c(1:4,4,5:9,9), sep=.) mydata - data.frame(x1,x2) crownames(mydata) - nam Many thanks in advance, Amor x1 - rnorm(11,5,1) x2 - runif(11,0,1) nam -paste(A, c(1:4,4,5:9,9), sep=.) mydata - data.frame(nam,x1,x2) mydata[!duplicated(mydata$nam),] nam x1 x2 1 A.1 5.418824 0.04867219 2 A.2 5.658403 0.18398337 3 A.3 4.458279 0.50975887 4 A.4 5.465920 0.16425144 6 A.5 3.769153 0.80380230 7 A.6 6.827979 0.13745441 8 A.7 5.353053 0.91418900 9 A.8 5.409866 0.41879708 10 A.9 5.041249 0.38226152 ?duplicated [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 for nice correlation output with significance symbols
On 5/10/2009 3:26 PM, Martin Batholdy wrote:
hi,
I am searching for a nice function which computes correlations out of a
data.frame and adds asterix-signs after each correlation when they are
significant...
...or a function which just show correlations greater than x in the output.
thanks!
This might be a starting point:
corstars - function(x){
require(Hmisc)
x - as.matrix(x)
R - rcorr(x)$r
p - rcorr(x)$P
mystars - ifelse(p .01, **|, ifelse(p .05, * |, |))
R - format(round(cbind(rep(-1.111, ncol(x)), R), 3))[,-1]
Rnew - matrix(paste(R, mystars, sep=), ncol=ncol(x))
diag(Rnew) - paste(diag(R), |, sep=)
rownames(Rnew) - colnames(x)
colnames(Rnew) - paste(colnames(x), |, sep=)
Rnew - as.data.frame(Rnew)
return(Rnew)
}
corstars(swiss[,1:4])
Fertility| Agriculture| Examination| Education|
Fertility 1.000 | 0.353* |-0.646**| -0.664**|
Agriculture 0.353* | 1.000 |-0.687**| -0.640**|
Examination -0.646**|-0.687**| 1.000 | 0.698**|
Education-0.664**|-0.640**| 0.698**| 1.000 |
At the very least, you could add a note that indicates what different
numbers of stars mean.
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] simple for loop question - how do you exit?
On 4/23/2009 2:29 PM, dre968 wrote: I have a loop and an if statement in the loop. once the if statement is true for 1 value in the loop i'd like to exit the loop. is there a command to do this? i know its going to be something like exit and i feel stupid asking this question Type ?if at the R prompt, which should load the help page for Control Flow. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Translate the elements of a dataframe
On 4/16/2009 2:58 PM, Juergen Rose wrote:
The second beginner question. I want to create a new dataframe, where
each element of the original dataframe is translated to 1 if it was +,
to 0 if it was - to -1 otherwise. I could do with:
Lines - abcd
+-+ +
+++ -
+1- '+ '
-++ +
+N- +
DF - read.table(textConnection(Lines), header = TRUE)
cnames - colnames(DF)
nrow -length(rownames(DF))
nc - length(cnames)
NDF - data.frame(matrix(c(rep(0,nc*nrow)),ncol=nc))
for (i in 1:length(cnames)) {
NDF[,i] - sapply(DF[,i],function(x) if (x==+) {1} else {if (x==-)
{0} else {-1}} )
}
colnames(NDF) - cnames
But this is shure one loop too much. Please give me the R way solution.
library(car)
DF[] - lapply(DF, function(x){recode(x, '+'=1; '-'=0; else=-1)})
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] excluding a column from a data frame
On 4/15/2009 1:38 AM, Erin Hodgess wrote: Dear R People: Suppose I have the following data frame: x1 x2 x3 1 -0.1582116 0.06635783 1.765448 2 -1.1407422 0.47235664 0.615931 3 0.8702362 2.32301341 2.653805 str(xx) 'data.frame': 3 obs. of 3 variables: $ x1: num -0.158 -1.141 0.87 $ x2: num 0.0664 0.4724 2.323 $ x3: num 1.765 0.616 2.654 I can exclude the second column nicely via: xx[,-2] x1 x3 1 -0.1582116 1.765448 2 -1.1407422 0.615931 3 0.8702362 2.653805 Now suppose I wanted to exclude the column called x2. If I try: xx[,-x2] Error in -x2 : invalid argument to unary operator things don't work. Is there a simple way to do this by name rather than number, please? Another way to do it is with subset(): subset(xx, select = -x2) Thanks, Erin -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] p-values of correlation matrix
On 4/13/2009 10:56 AM, thoeb wrote: Hello, I have a data frame containing several parameters. I want to investigate pair wise correlations between all of the parameters. For doing so I used the command cor(data.frame, method=”spearman”), the result is a matrix giving me the correlation coefficients of each pair, but not the p-values. Is it possible to get same matrix showing just the p-values instead of the correlation coefficients? As far as I know the command cor.test just works for two given vector, but not for a whole data frame. Thanks! library(Hmisc) ?rcorr - Tamara Hoebinger University of Vienna -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Efficiently create dummy
On 3/23/2009 5:44 AM, Rob Denniker wrote:
What's the neat way to create a dummy from a list?
The code below is not replicable, but hopefully self-explanatory...
d$treatment-rep(1,length(d))
notreat-c(AR, DE, MS, NY, TN, AK, LA, MD, NC, OK, UT,
VA)
#i would really like this to work:
d$treatment[d$st==any(notreat)]-0
#but instead i resort to this
for (i in 1:length(notreat)) {
temp.st - notreat[i]
d$treatment[d$st==temp.st]-0
i-i+1 }
d$treatment - as.numeric(!(d$st %in% notreat))
Thanks, 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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] lsmeans in R
On 3/11/2009 2:45 AM, suman Duvvuru wrote: I need help with calculating lsmeans (adjusted means) of different terms in a linear model including the main effect and the interaction effect terms. I use lm to run the linear models...I previously noted from literature that that effects package can be used to generate lsmeans. But I tried to use it but could not figure out which option to use to get means. If anyone can give an example of how to get lsmeans using lm object, that will very helpful. This R-help thread from March 2007 should help: http://finzi.psych.upenn.edu/R/Rhelp02/archive/95809.html Thanks, SUman [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Repeated values
On 3/11/2009 5:08 AM, Tammy Ma wrote: Hi,All. How to make a program to delete repeated value? a example c [1] 4 3 0 3 4 1 0 1 4 4 3 4 3 4 I want to get : 4 3 0 3 4 1 0 1 4 3 4 3 4 two 4 is being represented by one 4. x - c(4, 3, 0, 3, 4, 1, 0, 1, 4, 4, 3, 4, 3, 4) rle(x)$values [1] 4 3 0 3 4 1 0 1 4 3 4 3 4 ?rle Thanks. Tammy _ More than messages–check out the rest of the Windows Live™. http://www.microsoft.com/windows/windowslive/ [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Multilevel Modeling using glmmPQL
On 3/11/2009 11:29 AM, Howard Alper wrote:
Hi,
I'm trying to perform a power simulation for a simple multilevel
model, using the function glmmPQL in R version 2.8.1. I want to extract
the p-value for the fixed-effects portion of the regression, but I'm
having trouble doing that. I can extract the coefficients
(summary(fit)$coeff), and the covariance matrix (summary(fit)$varFix),
but I can't grab the p-value (or the t-statistic.) Could someone
explain how to do this?
Please send responses to hal...@health.nyc.gov.
library(MASS)
fm1 - glmmPQL(y ~ trt + I(week 2), random = ~ 1 | ID,
family = binomial, data = bacteria)
summary(fm1)$tTable[,5]
Look at str(summary(fm1)) .
Thanks in advance,
Howard
DISCLAIMER:\ Sample Disclaimer added in a VBScript.\ \...{{dropped:6}}
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Iterations of random sampling
On 3/11/2009 3:15 PM, René Pineda wrote: I have a univariate binary (1,0) 230,000 records, I need to make 20,000 iterations of random sampling of a fixed size. Where I put the result of the sum of selected records for each repetition X - rbinom(23, 1, .5) sample.sums - replicate(2, sum(sample(X, 10))) Thank's _ of your life [[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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] a few scatter plots for a specific correlation value
On 3/5/2009 6:46 AM, June Kim wrote:
Hello,
Is there a simple way to draw a few random sample scatter plots from a
given specific correlation coefficient(say, 0.18)?
scatter - function(n=100, r=.18){
require(MASS)
mymat - mvrnorm(n, mu=c(0,0), Sigma=matrix(c(1,r,r,1), ncol=2),
empirical=TRUE)
plot(mymat[,1], mymat[,2])
}
par(mfrow=c(2,2))
replicate(4, scatter())
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Doubt Linear Regression
On 3/5/2009 7:53 AM, Sueli Rodrigues wrote:
Hello. I have a file with 480 lines but each 6 lines corresponding just
one sample. How can can work out the linear regression to each 6 lines?
I use the model: model=lm(y~x)
mydf - data.frame(X = rnorm(480), Y = rnorm(480))
mydf$SAMPLE - rep(1:80, each=6)
by(mydf, mydf$SAMPLE, function(x){summary(lm(Y ~ X, data = x))})
OR
lapply(split(mydf, mydf$SAMPLE), function(x){summary(lm(Y ~ X, data = x))})
OR
library(nlme)
fm1 - lmList(Y ~ X | SAMPLE, mydf)
summary(fm1)
Sueli Rodrigues
Agronomy Eng. - UNESP
Master Degree - USP/ESALQ
PPG-Soils and Plants Nutrition
Phones(19)93442981
(19)33719762
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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 in writing my own function
On 2/23/2009 7:56 AM, tedzzx wrote: Dear all I am very intersted in writing my own function to deal with some complicated task, but I don't know how to start. I can't find detial material with examples teaching me how to write my own functions. Can anyone help me and recommend me some learning material? Chapter 10 of An Introduction to R would be a good place to start. http://cran.r-project.org/doc/manuals/R-intro.html Many Thanks Ted -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Doing pairwise comparisons using either Duncan, Tukey or LSD
On 2/19/2009 11:51 AM, Saeed Ahmadi wrote: Hi, I have a basic and simple question on how to code pairwise (multiple) mean compariosn between levels of a factor using one of the Duncan, Tukey or LSD. Here is one approach: library(multcomp) summary(glht(lm(Petal.Width ~ Species, data = iris), linfct = mcp(Species = Tukey))) Thanks in advance, Saeed -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Comparison of age categories using contrasts
On 2/16/2009 10:18 PM, Dylan Beaudette wrote:
On Mon, Feb 16, 2009 at 5:28 PM, Patrick Giraudoux
patrick.giraud...@univ-fcomte.fr wrote:
Greg Snow a écrit :
One approach is to create your own contrasts matrix:
mycmat - diag(8)
mycmat[ row(mycmat) == col(mycmat) + 1 ] - -1
mycmati - solve(mycmat)
contrasts(agefactor) - mycmati[,-1]
Now when you use agefactor, the intercept will be the first age group and
the slopes will be the differences between the pairs of groups (make sure
that the order of the levels of agefactor is correct).
The difference between this method and the contr.sdif function in MASS is
how the intercept will end up being interpreted (and the dimnames).
Hope this helps,
Actually, exactly what I needed including the reference to contr.sdif in
MASS I did not spot before (although I am a faithful reader of the
yellow book... but so many things still escape to me). Again thanks a lot.
Patrick
Glad you were able to solve your problem. Frank replied earlier with
the suggestion to use contrast.Design() to perform the tests after the
initial model has been fit. Perhaps I am a little denser than normal,
but I cannot figure out how to apply contrast.Design() to a similar
model with several levels of a factor. Example data:
# need these
library(lattice)
library(Design)
# replicate an important experimental dataset
set.seed(10101010)
x - rnorm(100)
y1 - x[1:25] * 2 + rnorm(25, mean=1)
y2 - x[26:50] * 2.6 + rnorm(25, mean=1.5)
y3 - x[51:75] * 2.9 + rnorm(25, mean=5)
y4 - x[76:100] * 3.5 + rnorm(25, mean=5.5)
d - data.frame(x=x, y=c(y1,y2,y3,y4), f=factor(rep(letters[1:4], each=25)))
# plot
xyplot(y ~ x, groups=f, data=d, auto.key=list(columns=4, title='Beard
Type', lines=TRUE, points=FALSE, cex=0.75), type=c('p','r'),
ylab='Number of Pirates', xlab='Distance from Land')
# standard comparison to base level of f
summary(lm(y ~ x * f, data=d))
# compare to level 4 of f
summary(lm(y ~ x * C(f, base=4), data=d))
# now try with contrast.Design():
dd - datadist(d)
options(datadist='dd')
l - ols(y ~ x * f, data=d)
# probably the wrong approach, as the results do not look too familiar:
contrast(l, list(f=levels(d$f)))
x f Contrast S.E. Lower Upper t Pr(|t|)
-0.3449623 a 0.3966682 0.1961059 0.007184856 0.7861515 2.02 0.0460
-0.3449623 b 0.5587395 0.1889173 0.183533422 0.9339456 2.96 0.0039
-0.3449623 c 4.1511677 0.1889170 3.775962254 4.5263730 21.97 0.
-0.3449623 d 4.3510407 0.120 3.975904726 4.7261766 23.04 0.
This is adjusting the output to a given value of 'x'... Am I trying to
use contrast.Design() for something that it was not intended for? Any
tips Frank?
Does this help?
# Compare with summary(lm(y ~ x * f, data=d))
contrast(l, a=list(f=levels(d$f)[2:4], x=0),
b=list(f=levels(d$f)[1], x=0))
f Contrast S.E. Lower Upper t Pr(|t|)
b 0.3673455 0.2724247 -0.1737135 0.9084046 1.35 0.1808
c 4.1310015 0.2714011 3.5919754 4.6700275 15.22 0.
d 4.4308653 0.2731223 3.8884207 4.9733098 16.22 0.
Error d.f.= 92
# Compare with summary(lm(y ~ x * C(f, base=4), data=d))
contrast(l, a=list(f=levels(d$f)[1:3], x=0),
b=list(f=levels(d$f)[4], x=0))
f Contrast S.E. Lower Upper t Pr(|t|)
a -4.4308653 0.2731223 -4.9733098 -3.8884207 -16.22 0.
b -4.0635198 0.2782704 -4.6161888 -3.5108507 -14.60 0.
c -0.2998638 0.2772684 -0.8505427 0.2508151 -1.08 0.2823
Error d.f.= 92
Cheers,
Dylan
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] wilcoxon test with bonferroni correction
On 2/1/2009 8:32 PM, Laura Lucia Prieto Godino wrote: Hi! I need to run a wilcoxon (Mann-whitly, in fact) test with bonferroni correction, as I am running 10 consecutive wilcoxon test not independent, and I know that bonferroni will partially correct for this problem, but I have no idea how to do it with R, I have been looking in the archive but couldn't understand how to do it. The format I am using at the moment is r4_o - [1] 1.05 2.60 1.57 3.07 1.20 1.00 2.11 1.10 0.10 r4_m - [1] 0 0 0 0 0 0 0 0 0 wilcoxon.test (r4_o, r3_m) Does any body know how to make the bonferroni correction when I compare them with the wilcoxon test? # Ten p-values X - seq(.001, .10, len=10) X [1] 0.001 0.012 0.023 0.034 0.045 0.056 0.067 0.078 0.089 0.100 # Same ten p-values adjusted by the Bonferroni method p.adjust(X, method=bonferroni) [1] 0.01 0.12 0.23 0.34 0.45 0.56 0.67 0.78 0.89 1.00 ?p.adjust Thank you very much. Lucia Lucia Prieto Godino PhD student. Department of Zoology, Downing street University of Cambridge. UK __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] bootstrapping in regression
On 1/29/2009 11:43 AM, Thomas Mang wrote: Hi, Please apologize if my questions sounds somewhat 'stupid' to the trained and experienced statisticians of you. Also I am not sure if I used all terms correctly, if not then corrections are welcome. I have asked myself the following question regarding bootstrapping in regression: Say for whatever reason one does not want to take the p-values for regression coefficients from the established test statistics distributions (t-distr for individual coefficients, F-values for whole-model-comparisons), but instead apply a more robust approach by bootstrapping. In the simple linear regression case, one possibility is to randomly rearrange the X/Y data pairs, estimate the model and take the beta1-coefficient. Do this many many times, and so derive the null distribution for beta1. Finally compare beta1 for the observed data against this null-distribution. What I now wonder is how the situation looks like in the multiple regression case. Assume there are two predictors, X1 and X2. Is it then possible to do the same, but just only rearranging the values of one predictor (the one of interest) at a time? Say I want again to test beta1. Is it then valid to many times randomly rearrange the X1 data (and keeping Y and X2 as observed), fit the model, take the beta1 coefficient, and finally compare the beta1 of the observed data against the distributions of these beta1s ? For X2, do the same, randomly rearrange X2 all the time while keeping Y and X1 as observed etc. Is this valid ? Second, if this is valid for the 'normal', fixed-effects only regression, is it also valid to derive null distributions for the regression coefficients of the fixed effects in a mixed model this way? Or does the quite different parameters estimation calculation forbid this approach (Forbid in the sense of bogus outcome) ? Thanks, Thomas Have a look at the following document by John Fox: http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-bootstrapping.pdf __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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 stepping through multiplie interactions
On 1/23/2009 12:44 PM, ppeetteerr wrote: I have a lm in R in the form model - lm( Z ~ A*B*C*D,data=mydata) I want to run the model and include all interactions expect the 4 way (A:B:C:D) is there an easy way of doing this? I then want to step down the model eliminating the non-significant terms I understand step() does this but how would I do it by hand? For the first part, try this: lm(Z ~ (A + B + C + D)^3, data = mydata) -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] creating a list of matrices or data frames
On 1/20/2009 11:34 AM, Simon Pickett wrote:
Hi all,
How would you create a list of data.frames within a loop, then bind all
the elements of the list using rbind?
take this example of matrices with differing numbers of rows
for(i in 1:3){
assign(paste(s,i, sep=),matrix(data = NA, nrow = i, ncol = 3, byrow
= FALSE, dimnames = NULL))
}
s1
s2
s3
I want to bind all the matrices at the end with do.call(rbind...)
rather than listing all the elements manually with rbind(s1,s2,s3...)
and so on.
thanks in advance.
df.list - vector(list, 3) # create list
for(i in 1:3){df.list[[i]] - matrix(data = NA,
nrow = i,
ncol = 3,
byrow = FALSE,
dimnames = NULL)}
do.call(rbind, df.list) # rbind list elements
Simon.
__
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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] lattice question: independent per-row or per-column scaling?
On 1/19/2009 8:51 AM, hadley wickham wrote: On Thu, Jan 8, 2009 at 11:25 AM, René J.V. Bertin rjvber...@gmail.com wrote: Hello - and happy newyear to all of you! I've got some data that I'm plotting with bwplot, a 3x2x3 design where the observable decreases with the principle independent factor, but at different rates. I'd like to get lattice to impose not a single set of axes ranges identical for all panels, but ranges that are identical for each panel row or each column. Effects will stand out much better like that. I've looked through the documentation of the latest lattice version, but I don't see a way to achieve this with a simple argument passed to bwplot. Can it be done otherwise and if so, how? The argument for xlim or ylim can be a list. Here is the key part of the help page for xyplot: xlim could also be a list, with as many components as the number of panels (recycled if necessary), with each component as described above. This is meaningful only when scales$x$relation is free or sliced, in which case these are treated as if they were the corresponding limit components returned by prepanel calculations. Here is a little example: library(lattice) mdf - data.frame(X1 = rep(LETTERS[1:3], each = 100*2*3), X2 = rep(c(J,K), 900), X3 = rep(LETTERS[24:26], 100*3*2), Y = c(runif(600, min=.01,max=.32), runif(600, min=.33,max=.65), runif(600, min=.66,max=.99))) bwplot(Y ~ X3 | X2*X1, data = mdf, layout=c(2,3,1), ylim=as.data.frame(matrix(c(.01,.32, .01,.32, .33,.65, .33,.65, .66,.99, .66,.99), nrow=2)), scales=list(y=list(relation=free))) It's not lattice, but you can do this with ggplot2 - see the examples for http://had.co.nz/ggplot2/facet_grid.html Regards, Hadley -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Extraction from an output
On 1/12/2009 6:42 AM, robert-mcfad...@o2.pl wrote: Hello, Would you tell my how to extract a result from a test - it's justified because I need to run this test many times. Here is an example from authors' test: library(coin) lungtumor - data.frame(dose = rep(c(0, 1, 2), c(40, 50, 48)), tumor = c(rep(c(0, 1), c(38, 2)), rep(c(0, 1), c(43, 7)), rep(c(0, 1), c(33, 15 ca.test-independence_test(tumor ~ dose, data = lungtumor, teststat = quad) ca.test Asymptotic General Independence Test data: tumor by dose chi-squared = 10.6381, df = 1, p-value = 0.001108 How to use ca.test and extract p-value and chi-squared. pvalue(ca.test) [1] 0.001107836 statistic(ca.test) [1] 10.63806 See the Value subsection on the help page for independence_test(). Robert __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] anova() or aov()?
On 1/12/2009 8:57 AM, j...@in.gr wrote: Dear all, I have a simple (I think) question that is troubling me lately: Is there any main difference between anova() command and aov() command when performing an ANOVA in Experimental design analyses? The main difference is that aov() *fits* a single model and anova() *summarizes* a single model or compares two or more models. Thank you for your time, Ismini __ 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. -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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] Applying 'lm' in each group
On 1/10/2009 12:55 PM, Bhargab Chattopadhyay wrote:
Hi,
I want to do regression in each group. I made the group the following way.
Now I want to run regression in each of the three groups..
Suppose ,
x-c(0.5578196,6.5411662,13.2728619,2.0217271,6.7216176,3.37220617,2.5773252,7.2600583,15.3731026,3.4140288,8.1335874,
15..6476637,4.3014084,9.1224379,18.5605355,4.7448394,11.9296663,18.5866354,12.3797674,18.7572812,2.70433816,2.88924220,
2.94688208,3.37154364,2.26311786,3.31002593)
y-c(18.63654,233.55387,152.61107,103.49646,234.55054,14.2767453,160.78742,150.35391,223.89225,168.55567,190.51031,
227.68339,152.42658,208.70115, 223.91982, 221.70702, 213.71135,168.0199,
222.69040,228.49353, 164.95750,243.18828,
229.94688,313.37154364,202.263786,139.31002593)
n-length(x)
m-3;
r-8;
c-r*(m-1);
g1-rep(1:(m-1),each=r);
g-c(g1,rep(m,n-c))
xg-split(x,g);
yg-split(y,g);
Now if I write lm(yg~xg) then it won't work. Please advice how to proceed.
mdf - data.frame(y,x,g)
lapply(split(mdf, g), function(X){lm(y ~ x, data = X)})
$`1`
Call:
lm(formula = y ~ x, data = X)
Coefficients:
(Intercept)x
76.1110.85
$`2`
Call:
lm(formula = y ~ x, data = X)
Coefficients:
(Intercept)x
166.6933.580
$`3`
Call:
lm(formula = y ~ x, data = X)
Coefficients:
(Intercept)x
218.412 -0.735
Thanks in advance.
Bhargab
[[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.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Merging two data frame with different nrows
On 12/21/2008 8:32 PM, Antonio Paredes wrote:
hello-
I'm trying to merge two data frames with different number of rows. I started
with a for loop, but it just taking to long . I wanted to ask if there is
any function in R to merge data frames with different number of rows.
As the Posting Guide asks, you can use RSiteSearch() or help.search()
before posting to search for functionality. For example:
RSiteSearch('merge', restrict='function')
The third hit points to merge(), which does what you want.
--
Chuck Cleland, Ph.D.
NDRI, Inc. (www.ndri.org)
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
__
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] Extract Data from a Webpage
Hi All: I would like to extract the provider name, address, and phone number from multiple webpages like this: http://oasasapps.oasas.state.ny.us/portal/pls/portal/oasasrep.providersearch.take_to_rpt?P1=3489P2=11490 Based on searching R-help archives, it seems like the XML package might have something useful for this task. I can load the XML package and supply the url as an argument to htmlTreeParse(), but I don't know how to go from there. thanks, Chuck Cleland sessionInfo() R version 2.8.0 Patched (2008-12-04 r47066) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] XML_1.98-1 -- Chuck Cleland, Ph.D. NDRI, Inc. (www.ndri.org) 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894 __ 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.
