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 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] 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] 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
> Vsa 0.0092181 0.0054564 1.68941 9.1140e-02 Sa <--> Sa
> Cisa -0.0063651 0.0051128 -1.24492 2.1316e-01 Sa <--> Ia
> Vip 0.0696837 0.0103795 6.71357 1.8991e-11 Ip <--> Ip
> Vsp 0.0284726 0.0089274 3.18936 1.4259e-03 Sp <--> Sp
> Cisp 0.0011771 0.0071251 0.16521 8.6878e-01 Sp <--> Ip
> Vd1 0.0480379 0.0063780 7.53177 4.9960e-14 ALC1 <--> ALC1
> Vd2 0.0762156 0.0044523 17.11821 0.e+00 ALC2 <--> ALC2
> Vd3 0.0762794 0.0097763 7.80249 5.9952e-15 ALC3 <--> ALC3
> Vd4 0.1057875 0.0108526 9.74770 0.e+00 PEER1 <--> PEER1
> Vd5 0.1712811 0.0087037 19.67904 0.e+00 PEER2 <--> PEER2
> Vd6 0.1289592 0.0177027 7.28471 3.2241e-13 PEER3 <--> PEER3
> Cd1 0.0109322 0.0061562 1.77578 7.5769e-02 PEER1 <--> ALC1
> Cd2 0.0339991 0.0046391 7.32874 2.3226e-13 PEER2 <--> ALC2
> Cd3 0.0374125 0.0101878 3.67229 2.4038e-04 PEER3 <--> ALC3
>
> Iterations = 139
>
> snip -
>
> Regards,
> John
>
>
> John Fox
> Senator William McMaster
> Professor of Social Statistics
> 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
>>
Re: [R] sem package and growth curves
gt; 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 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] 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 wrote:
>
>> From: Karsten Wolf
>> 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] 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] 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] 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] 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] 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] 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 "M16"<"M05"<"M02"<..: 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")= - attr(*, "formula")=Class 'formula' length 3 distance ~ age | Subject .. ..- attr(*, ".Environment")= - 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] 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] 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] 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] using "cochran.test()" as a "mcnemar.test()" ?
On 6/1/2009 6:01 PM, Tal Galili wrote:
> Hi Chuck,
> Thank you for the fast reply.
> I looked through that thread (that's how I came across the commands I
> wrote).
>
> 1) My question was if this command (cochran.test) could be used the same
> way as (mcnemar.test). It obviously can't be, but maybe there is another
> command I don't know about.
> 2) Why do you think I might not want this test ?
You posted a 5x5 matrix with values ranging from 0-3. That does not
fit with assumptions 1 and 3 described on the Wikipedia page you posted
(a large number of blocks and outcomes of the treatments coded as binary
responses).
hope this helps,
Chuck
> Thanks,
> Tal
>
>
>
>
> On Tue, Jun 2, 2009 at 12:56 AM, Chuck Cleland <mailto:cclel...@optonline.net>> wrote:
>
> On 6/1/2009 5:29 PM, Tal Galili wrote:
> > Hello all
> > I wish to perform a mcnemar.test() for a 5X5 matrix.
> > Wikipedia tells me (http://en.wikipedia.org/wiki/Cochran_test) I
> should turn
> > to cochran.test.
> > The only place I found it was in the "outliers" package, but the
> command
> > cochran.test() acts differently then mcnemar.test() , and doesn't
> take a
> > table as input.
> >
> > Any ideas on how to use it ?
> >
> > #Example code:
> > aa =
> >
> as.table(matrix(c(0,0,0,0,1,0,2,2,1,0,1,2,2,1,0,1,1,3,2,0,0,0,1,2,3),5,5))
> > mcnemar.test(aa)
> > # p-value = NA
> > install.packages("outliers")
> > require(outliers)
> > cochran.test(aa)
> > # Error in x$terms : $ operator is invalid for atomic vectors
> >
> > Thanks,
> > Tal
>
> I'm not sure you actually want Cochran's Q Test, but if you do see
> this thread:
>
> https://stat.ethz.ch/pipermail/r-help/2006-September/113139.html
>
> --
> Chuck Cleland, Ph.D.
> NDRI, Inc. (www.ndri.org <http://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
>
>
>
>
> --
> --
>
>
> My contact information:
> Tal Galili
> Phone number: 972-50-3373767
> FaceBook: Tal Galili
> My Blogs:
> http://www.r-statistics.com/
> http://www.talgalili.com
> http://www.biostatistics.co.il
--
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 "cochran.test()" as a "mcnemar.test()" ?
On 6/1/2009 5:29 PM, Tal Galili wrote:
> Hello all
> I wish to perform a mcnemar.test() for a 5X5 matrix.
> Wikipedia tells me (http://en.wikipedia.org/wiki/Cochran_test) I should turn
> to cochran.test.
> The only place I found it was in the "outliers" package, but the command
> cochran.test() acts differently then mcnemar.test() , and doesn't take a
> table as input.
>
> Any ideas on how to use it ?
>
> #Example code:
> aa =
> as.table(matrix(c(0,0,0,0,1,0,2,2,1,0,1,2,2,1,0,1,1,3,2,0,0,0,1,2,3),5,5))
> mcnemar.test(aa)
> # p-value = NA
> install.packages("outliers")
> require(outliers)
> cochran.test(aa)
> # Error in x$terms : $ operator is invalid for atomic vectors
>
> Thanks,
> Tal
I'm not sure you actually want Cochran's Q Test, but if you do see
this thread:
https://stat.ethz.ch/pipermail/r-help/2006-September/113139.html
--
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] 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] 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] 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] 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] 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] 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] 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
> 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 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] 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] 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] 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.
