subject:"\[R\] How to calculate\?"

[R] How to calculate relative risk from GAM model in mgcv package?

2023-12-06 Thread Eren Yalın

Hi R users,I am a beginner in the use of R. I need urgent help for my
thesis study.


I have daily air pollution parameters PM10, PM2.5 CO, NO2, SO2, and O3. I
also have daily hospital admission numbers. Taking into account the effect
of weekends and holidays, I would like to used generalised additive model
(GAM) to explore the relationship between daily patients admissions, and
air pollution parameters. I would like tu use mgcv package. How to get
overall relative risk and 95%CI for every pollutant?

I don't know if it's correct but here are the codes I used:

install.packages("mgcv")

library(mgcv)

data=read.csv2(file.choose(),header=TRUE)

data$date <- as.Date(data$date, format="%d.%m.%Y")


data$weekend <- factor(data$weekend)


data$holiday <- factor(data$holiday)


model <- gam(adm ~ s(PM10, k = 5) + s(PM2.5, k = 5) + s(CO, k = 5) + s(NO2,
k = 5) + s(SO2, k = 5) + s(O3, k = 5) + weekend + holiday, family =
quasipoisson(link = "log"), data = data, method = "REML")

pred <- predict.gam (model, type = "response")

relative_risk <- exp(pred$fit)

However, when I look at the results, it calculates RR for 365 days
separately.

How can I get a result like the table 4 in this article (
https://pubmed.ncbi.nlm.nih.gov/36161569/)? There is only one RR
calculation for each pollutant.

I would be very grateful if you could help me. Thank you.


Eren YALIN M.D., Research Assistant

 . University, Medical Faculty,  Department of Public Health

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the Mean by Multiple Groups in R

2023-10-24 Thread Gabor Grothendieck

A variation is to remove Well and then we can use dot to refer to the
remaining columns.

  aggregate(cbind(OD, ODnorm)  ~ . , subset(df, select = - Well), mean)


On Tue, Oct 24, 2023 at 8:32 AM Luigi Marongiu  wrote:
>
> Hello,
> I have a data frame with different groups (Time, Target, Conc) and
> each entry has a triplicate value of the measurements OD and ODnorm.
> How can I merge the triplicates into a single mean value?
> I tried the following:
> ```
> df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
> OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
> Target=rep("BACT", 9),
> Conc=c(1,1,1,2,2,2,3,3,3),
> ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
> stringsAsFactors = FALSE)
> aggregate(.~ODnorm, df, mean)
>
> > aggregate(.~ODnorm, df, mean)
>   ODnorm Time Well OD Target Conc
> 1  0   NA   NA NA NA   NA
> 2  6   NA   NA NA NA   NA
> 3  9   NA   NA NA NA   NA
> 4 45   NA   NA NA NA   NA
> 5 48   NA   NA NA NA   NA
> 6 82   NA   NA NA NA   NA
> 7138   NA   NA NA NA   NA
> 8158   NA   NA NA NA   NA
>
>  aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
>   ODnorm Time Target Conc
> 1  0   NA NA   NA
> 2  6   NA NA   NA
> 3  9   NA NA   NA
> 4 45   NA NA   NA
> 5 48   NA NA   NA
> 6 82   NA NA   NA
> 7138   NA NA   NA
> 8158   NA NA   NA
> ```
>
> Thank you.
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the Mean by Multiple Groups in R

2023-10-24 Thread Luigi Marongiu

Thank you

On Tue, Oct 24, 2023 at 3:01 PM peter dalgaard  wrote:
>
> Also,
>
> > aggregate(cbind(OD, ODnorm) ~ Time + Target + Conc, data = df, FUN = "mean")
>   Time Target Conc   ODODnorm
> 11   BACT1 765. 108.3
> 21   BACT2 745.  88.3
> 31   BACT3 675.  18.0
>
> (You might wish for "cbind(OD,ODnorm) ~ . - Well", but aggregate.formula is 
> not smart enough for that.)
>
> -pd
>
> > On 24 Oct 2023, at 14:40 , Sarah Goslee  wrote:
> >
> > Hi,
> >
> > I think you're misunderstanding which set of variables go on either
> > side of the formula.
> >
> > Is this what you're looking for?
> >
> >> aggregate(OD ~ Time + Target + Conc, data = df, FUN = "mean")
> >  Time Target Conc   OD
> > 11   BACT1 765.
> > 21   BACT2 745.
> > 31   BACT3 675.
> >> aggregate(ODnorm ~ Time + Target + Conc, data = df, FUN = "mean")
> >  Time Target ConcODnorm
> > 11   BACT1 108.3
> > 21   BACT2  88.3
> > 31   BACT3  18.0
> >
> > Or using a different form, that might be more straightforward to you:
> >
> >> aggregate(df[, c("OD", "ODnorm")], by = df[, c("Time", "Target", "Conc")], 
> >> data = df, FUN = "mean")
> >  Time Target Conc   ODODnorm
> > 11   BACT1 765. 108.3
> > 21   BACT2 745.  88.3
> > 31   BACT3 675.  18.0
> >
> > Sarah
> >
> > On Tue, Oct 24, 2023 at 8:31 AM Luigi Marongiu  
> > wrote:
> >>
> >> Hello,
> >> I have a data frame with different groups (Time, Target, Conc) and
> >> each entry has a triplicate value of the measurements OD and ODnorm.
> >> How can I merge the triplicates into a single mean value?
> >> I tried the following:
> >> ```
> >> df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
> >>OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
> >>Target=rep("BACT", 9),
> >>Conc=c(1,1,1,2,2,2,3,3,3),
> >>ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
> >>stringsAsFactors = FALSE)
> >> aggregate(.~ODnorm, df, mean)
> >>
> >>> aggregate(.~ODnorm, df, mean)
> >>  ODnorm Time Well OD Target Conc
> >> 1  0   NA   NA NA NA   NA
> >> 2  6   NA   NA NA NA   NA
> >> 3  9   NA   NA NA NA   NA
> >> 4 45   NA   NA NA NA   NA
> >> 5 48   NA   NA NA NA   NA
> >> 6 82   NA   NA NA NA   NA
> >> 7138   NA   NA NA NA   NA
> >> 8158   NA   NA NA NA   NA
> >>
> >> aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
> >>  ODnorm Time Target Conc
> >> 1  0   NA NA   NA
> >> 2  6   NA NA   NA
> >> 3  9   NA NA   NA
> >> 4 45   NA NA   NA
> >> 5 48   NA NA   NA
> >> 6 82   NA NA   NA
> >> 7138   NA NA   NA
> >> 8158   NA NA   NA
> >> ```
> >>
> >> Thank you.
> >>
> >> __
> >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide 
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >
> >
> >
> > --
> > Sarah Goslee (she/her)
> > http://www.numberwright.com
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd@cbs.dk  Priv: pda...@gmail.com
>


-- 
Best regards,
Luigi

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the Mean by Multiple Groups in R

2023-10-24 Thread peter dalgaard

Also,

> aggregate(cbind(OD, ODnorm) ~ Time + Target + Conc, data = df, FUN = "mean")
  Time Target Conc   ODODnorm
11   BACT1 765. 108.3
21   BACT2 745.  88.3
31   BACT3 675.  18.0

(You might wish for "cbind(OD,ODnorm) ~ . - Well", but aggregate.formula is not 
smart enough for that.)

-pd

> On 24 Oct 2023, at 14:40 , Sarah Goslee  wrote:
> 
> Hi,
> 
> I think you're misunderstanding which set of variables go on either
> side of the formula.
> 
> Is this what you're looking for?
> 
>> aggregate(OD ~ Time + Target + Conc, data = df, FUN = "mean")
>  Time Target Conc   OD
> 11   BACT1 765.
> 21   BACT2 745.
> 31   BACT3 675.
>> aggregate(ODnorm ~ Time + Target + Conc, data = df, FUN = "mean")
>  Time Target ConcODnorm
> 11   BACT1 108.3
> 21   BACT2  88.3
> 31   BACT3  18.0
> 
> Or using a different form, that might be more straightforward to you:
> 
>> aggregate(df[, c("OD", "ODnorm")], by = df[, c("Time", "Target", "Conc")], 
>> data = df, FUN = "mean")
>  Time Target Conc   ODODnorm
> 11   BACT1 765. 108.3
> 21   BACT2 745.  88.3
> 31   BACT3 675.  18.0
> 
> Sarah
> 
> On Tue, Oct 24, 2023 at 8:31 AM Luigi Marongiu  
> wrote:
>> 
>> Hello,
>> I have a data frame with different groups (Time, Target, Conc) and
>> each entry has a triplicate value of the measurements OD and ODnorm.
>> How can I merge the triplicates into a single mean value?
>> I tried the following:
>> ```
>> df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
>>OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
>>Target=rep("BACT", 9),
>>Conc=c(1,1,1,2,2,2,3,3,3),
>>ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
>>stringsAsFactors = FALSE)
>> aggregate(.~ODnorm, df, mean)
>> 
>>> aggregate(.~ODnorm, df, mean)
>>  ODnorm Time Well OD Target Conc
>> 1  0   NA   NA NA NA   NA
>> 2  6   NA   NA NA NA   NA
>> 3  9   NA   NA NA NA   NA
>> 4 45   NA   NA NA NA   NA
>> 5 48   NA   NA NA NA   NA
>> 6 82   NA   NA NA NA   NA
>> 7138   NA   NA NA NA   NA
>> 8158   NA   NA NA NA   NA
>> 
>> aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
>>  ODnorm Time Target Conc
>> 1  0   NA NA   NA
>> 2  6   NA NA   NA
>> 3  9   NA NA   NA
>> 4 45   NA NA   NA
>> 5 48   NA NA   NA
>> 6 82   NA NA   NA
>> 7138   NA NA   NA
>> 8158   NA NA   NA
>> ```
>> 
>> Thank you.
>> 
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
> 
> 
> -- 
> Sarah Goslee (she/her)
> http://www.numberwright.com
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the Mean by Multiple Groups in R

2023-10-24 Thread Luigi Marongiu

Thank you, the last is exactly what I was looking for.

On Tue, Oct 24, 2023 at 2:41 PM Sarah Goslee  wrote:
>
> Hi,
>
> I think you're misunderstanding which set of variables go on either
> side of the formula.
>
> Is this what you're looking for?
>
> > aggregate(OD ~ Time + Target + Conc, data = df, FUN = "mean")
>   Time Target Conc   OD
> 11   BACT1 765.
> 21   BACT2 745.
> 31   BACT3 675.
> > aggregate(ODnorm ~ Time + Target + Conc, data = df, FUN = "mean")
>   Time Target ConcODnorm
> 11   BACT1 108.3
> 21   BACT2  88.3
> 31   BACT3  18.0
>
> Or using a different form, that might be more straightforward to you:
>
> > aggregate(df[, c("OD", "ODnorm")], by = df[, c("Time", "Target", "Conc")], 
> > data = df, FUN = "mean")
>   Time Target Conc   ODODnorm
> 11   BACT1 765. 108.3
> 21   BACT2 745.  88.3
> 31   BACT3 675.  18.0
>
> Sarah
>
> On Tue, Oct 24, 2023 at 8:31 AM Luigi Marongiu  
> wrote:
> >
> > Hello,
> > I have a data frame with different groups (Time, Target, Conc) and
> > each entry has a triplicate value of the measurements OD and ODnorm.
> > How can I merge the triplicates into a single mean value?
> > I tried the following:
> > ```
> > df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
> > OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
> > Target=rep("BACT", 9),
> > Conc=c(1,1,1,2,2,2,3,3,3),
> > ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
> > stringsAsFactors = FALSE)
> > aggregate(.~ODnorm, df, mean)
> >
> > > aggregate(.~ODnorm, df, mean)
> >   ODnorm Time Well OD Target Conc
> > 1  0   NA   NA NA NA   NA
> > 2  6   NA   NA NA NA   NA
> > 3  9   NA   NA NA NA   NA
> > 4 45   NA   NA NA NA   NA
> > 5 48   NA   NA NA NA   NA
> > 6 82   NA   NA NA NA   NA
> > 7138   NA   NA NA NA   NA
> > 8158   NA   NA NA NA   NA
> >
> >  aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
> >   ODnorm Time Target Conc
> > 1  0   NA NA   NA
> > 2  6   NA NA   NA
> > 3  9   NA NA   NA
> > 4 45   NA NA   NA
> > 5 48   NA NA   NA
> > 6 82   NA NA   NA
> > 7138   NA NA   NA
> > 8158   NA NA   NA
> > ```
> >
> > Thank you.
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>
>
> --
> Sarah Goslee (she/her)
> http://www.numberwright.com



-- 
Best regards,
Luigi

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the Mean by Multiple Groups in R

2023-10-24 Thread Sarah Goslee

Hi,

I think you're misunderstanding which set of variables go on either
side of the formula.

Is this what you're looking for?

> aggregate(OD ~ Time + Target + Conc, data = df, FUN = "mean")
  Time Target Conc   OD
11   BACT1 765.
21   BACT2 745.
31   BACT3 675.
> aggregate(ODnorm ~ Time + Target + Conc, data = df, FUN = "mean")
  Time Target ConcODnorm
11   BACT1 108.3
21   BACT2  88.3
31   BACT3  18.0

Or using a different form, that might be more straightforward to you:

> aggregate(df[, c("OD", "ODnorm")], by = df[, c("Time", "Target", "Conc")], 
> data = df, FUN = "mean")
  Time Target Conc   ODODnorm
11   BACT1 765. 108.3
21   BACT2 745.  88.3
31   BACT3 675.  18.0

Sarah

On Tue, Oct 24, 2023 at 8:31 AM Luigi Marongiu  wrote:
>
> Hello,
> I have a data frame with different groups (Time, Target, Conc) and
> each entry has a triplicate value of the measurements OD and ODnorm.
> How can I merge the triplicates into a single mean value?
> I tried the following:
> ```
> df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
> OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
> Target=rep("BACT", 9),
> Conc=c(1,1,1,2,2,2,3,3,3),
> ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
> stringsAsFactors = FALSE)
> aggregate(.~ODnorm, df, mean)
>
> > aggregate(.~ODnorm, df, mean)
>   ODnorm Time Well OD Target Conc
> 1  0   NA   NA NA NA   NA
> 2  6   NA   NA NA NA   NA
> 3  9   NA   NA NA NA   NA
> 4 45   NA   NA NA NA   NA
> 5 48   NA   NA NA NA   NA
> 6 82   NA   NA NA NA   NA
> 7138   NA   NA NA NA   NA
> 8158   NA   NA NA NA   NA
>
>  aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
>   ODnorm Time Target Conc
> 1  0   NA NA   NA
> 2  6   NA NA   NA
> 3  9   NA NA   NA
> 4 45   NA NA   NA
> 5 48   NA NA   NA
> 6 82   NA NA   NA
> 7138   NA NA   NA
> 8158   NA NA   NA
> ```
>
> Thank you.
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



-- 
Sarah Goslee (she/her)
http://www.numberwright.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to Calculate the Mean by Multiple Groups in R

2023-10-24 Thread Luigi Marongiu

Hello,
I have a data frame with different groups (Time, Target, Conc) and
each entry has a triplicate value of the measurements OD and ODnorm.
How can I merge the triplicates into a single mean value?
I tried the following:
```
df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
Target=rep("BACT", 9),
Conc=c(1,1,1,2,2,2,3,3,3),
ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
stringsAsFactors = FALSE)
aggregate(.~ODnorm, df, mean)

> aggregate(.~ODnorm, df, mean)
  ODnorm Time Well OD Target Conc
1  0   NA   NA NA NA   NA
2  6   NA   NA NA NA   NA
3  9   NA   NA NA NA   NA
4 45   NA   NA NA NA   NA
5 48   NA   NA NA NA   NA
6 82   NA   NA NA NA   NA
7138   NA   NA NA NA   NA
8158   NA   NA NA NA   NA

 aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
  ODnorm Time Target Conc
1  0   NA NA   NA
2  6   NA NA   NA
3  9   NA NA   NA
4 45   NA NA   NA
5 48   NA NA   NA
6 82   NA NA   NA
7138   NA NA   NA
8158   NA NA   NA
```

Thank you.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the derivatives at each data point?

2023-01-31 Thread konstantinos christodoulou

Thank you Petr!

On Tue, Jan 31, 2023 at 11:58 AM PIKAL Petr  wrote:

> Hi Konstantinos
>
> Not exactly derivative but
> > diff(df[,2])
> [1] -0.01 -0.01 -0.01 -0.01  0.00  0.01 -0.02 -0.03 -0.02
>
> May be enaough for you.
>
> Cheers
> Petr
>
> >
> > -Original Message-
> > From: R-help  On Behalf Of konstantinos
> > christodoulou
> > Sent: Tuesday, January 31, 2023 10:16 AM
> > To: r-help mailing list 
> > Subject: [R] How to calculate the derivatives at each data point?
> >
> > Hi everyone,
> >
> > I have a vector with atmospheric measurements (x-axis) that is
> > obtained/calculated at different altitudes (y-axis). The altitude is
> uniformly
> > distributed every 7 meters.
> > For example my dataframe is:
> > df <- dataframe(
> > *altitude* = c(1005, 1012, 1019, 1026, 1033, 1040, 1047, 1054, 1061,
> 1068),
> > *atm_values* = c(1.41, 1.40, 1.39, 1.38, 1.37, 1.37, 1.38, 1.36, 1.33,
> 1.31)
> >  )
> >
> > How can I find the derivatives of the atmospheric measurements at each
> > altitude?
> >
> > I look forward to hearing from you!
> >
> > Thanks,
> > Kostas
> >
> >   [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the derivatives at each data point?

2023-01-31 Thread konstantinos christodoulou

Hi Ivan!

Thank you for your valuable insights! I look forward to learning more about
numerical differentiation and about this subject.

The pracma package and the fornberg() function is impressive. I got some
really good approximations on my derivatives.

Thank you!
Kostas

On Tue, Jan 31, 2023 at 12:18 PM Ivan Krylov  wrote:

> В Tue, 31 Jan 2023 11:16:21 +0200
> konstantinos christodoulou 
> пишет:
>
> > How can I find the derivatives of the atmospheric measurements at each
> > altitude?
>
> Welcome to the world of finite difference methods! If you can find a
> good textbook on them, it may be a good idea to skim it.
>
> pracma::fornberg() will give you a numerically stable approximation
> (otherwise the Vandermonde matrix required to obtain the Taylor series
> coefficients may get hard to solve) to the derivative values you're
> interested in, but do note that they are only approximations. In
> particular, there's less information for the values at the ends of
> the altitude range than for the points in the middle.
>
> >   [[alternative HTML version deleted]]
>
> P.S. Please compose your messages to R-help in plain text:
> https://www.r-project.org/posting-guide.html
> https://stat.ethz.ch/pipermail/r-help/2023-January/476845.html
>
> --
> Best regards,
> Ivan
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the derivatives at each data point?

2023-01-31 Thread Ivan Krylov

В Tue, 31 Jan 2023 11:16:21 +0200
konstantinos christodoulou 
пишет:

> How can I find the derivatives of the atmospheric measurements at each
> altitude?

Welcome to the world of finite difference methods! If you can find a
good textbook on them, it may be a good idea to skim it.

pracma::fornberg() will give you a numerically stable approximation
(otherwise the Vandermonde matrix required to obtain the Taylor series
coefficients may get hard to solve) to the derivative values you're
interested in, but do note that they are only approximations. In
particular, there's less information for the values at the ends of
the altitude range than for the points in the middle.

>   [[alternative HTML version deleted]]

P.S. Please compose your messages to R-help in plain text:
https://www.r-project.org/posting-guide.html
https://stat.ethz.ch/pipermail/r-help/2023-January/476845.html

-- 
Best regards,
Ivan

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the derivatives at each data point?

2023-01-31 Thread PIKAL Petr

Hi Konstantinos

Not exactly derivative but
> diff(df[,2])
[1] -0.01 -0.01 -0.01 -0.01  0.00  0.01 -0.02 -0.03 -0.02

May be enaough for you.

Cheers
Petr

>
> -Original Message-
> From: R-help  On Behalf Of konstantinos
> christodoulou
> Sent: Tuesday, January 31, 2023 10:16 AM
> To: r-help mailing list 
> Subject: [R] How to calculate the derivatives at each data point?
> 
> Hi everyone,
> 
> I have a vector with atmospheric measurements (x-axis) that is
> obtained/calculated at different altitudes (y-axis). The altitude is
uniformly
> distributed every 7 meters.
> For example my dataframe is:
> df <- dataframe(
> *altitude* = c(1005, 1012, 1019, 1026, 1033, 1040, 1047, 1054, 1061,
1068),
> *atm_values* = c(1.41, 1.40, 1.39, 1.38, 1.37, 1.37, 1.38, 1.36, 1.33,
1.31)
>  )
> 
> How can I find the derivatives of the atmospheric measurements at each
> altitude?
> 
> I look forward to hearing from you!
> 
> Thanks,
> Kostas
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate the derivatives at each data point?

2023-01-31 Thread konstantinos christodoulou

Hi everyone,

I have a vector with atmospheric measurements (x-axis) that is
obtained/calculated at different altitudes (y-axis). The altitude is
uniformly distributed every 7 meters.
For example my dataframe is:
df <- dataframe(
*altitude* = c(1005, 1012, 1019, 1026, 1033, 1040, 1047, 1054, 1061, 1068),
*atm_values* = c(1.41, 1.40, 1.39, 1.38, 1.37, 1.37, 1.38, 1.36, 1.33, 1.31)
 )

How can I find the derivatives of the atmospheric measurements at each
altitude?

I look forward to hearing from you!

Thanks,
Kostas

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate area and volume of an image with R

2021-02-28 Thread Martin Møller Skarbiniks Pedersen

On Sun, Feb 28, 2021, 17:49 Paul Bernal  wrote:

> Hello everyone,
>
> Is there a way to calculate volume and area of an image with R?

How can an image have a volume or an area?
I think you need to be more specific.

Regards
Martin

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate area and volume of an image with R

2021-02-28 Thread Paul Bernal

Hello everyone,

Is there a way to calculate volume and area of an image with R?

Any guidance will be greatly appreciated.

Best regards,

Paul

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 odd ratios with R?

2020-07-07 Thread Minato Nakazawa



A review (sorry in Japanese) on the calculation of odds ratios 
with confidence intervals using several packages in R is given
by Prof. Okumura, Mie Univ.

https://oku.edu.mie-u.ac.jp/~okumura/stat/2by2.html

x <- matrix(c(6, 3, 4, 5), 2)
# using vcd
library(vcd)
res <- oddsratio(x)
exp(confint(res))
summary(res)
# Note: summary(oddsratio(x, log=FALSE)) gives inappropriate p-value.
# using fmsb
library(fmsb)
oddsratio(x, p.calc.by.independence=FALSE)

Best,
Minato Nakazawa

On Mon, 6 Jul 2020 15:01:34 +0200
Luigi Marongiu  wrote:

> Hello,
> Is it possible to calculate with a single function the odd ratios?
> Now I can use this implement:
> ```
> or <- (De/He)/(Dn/Hn) # Disease exposed, Healthy non-exposed
> logo <- log(or)
> x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
> lower_ci = exp(logo - 1.96*x)
> upper_ci = exp(logo + 1.96*x)
> cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), 
> ")",
> spe = "")
> ```
> for instance,
> ```
> De <-6
> Dn <-3
> He <-4
> Hn <-5
> or <- (De/He)/(Dn/Hn)
> logo <- log(or)
> x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
> lower_ci = exp(logo - 1.96*x)
> upper_ci = exp(logo + 1.96*x)
> cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), 
> ")",
> spe = "")
> > OR: 2.5 ( 0.37 - 16.889 )
> ```
> Is there a simple function from some package that can also add a
> p-value to this test? Or how can I calculate the p-value on my own?
> -- 
> Best regards,
> Luigi
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 


-- 
Minato Nakazawa 
Professor, Division of Global Health, Department of Public Health,
Kobe University Graduate School of Health Sciences
[web] http://minato.sip21c.org/
[phone] +81-78-796-4551
[mobile e-mail] minatonakaz...@gmail.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 odd ratios with R?

2020-07-06 Thread Stefan Evert

fisher.test() computes exact confidence intervals for the odds ratio.

> On 6 Jul 2020, at 15:01, Luigi Marongiu  wrote:
> 
> Is there a simple function from some package that can also add a
> p-value to this test? Or how can I calculate the p-value on my own?

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 odd ratios with R?

2020-07-06 Thread Michael Dewey


Dear Luigi

You could try the epitools package which gives a large number of ways of 
doing this. I would have thought that using Wald intervals for the log 
odds ration was not optimal with small frequencies.


Michael

On 06/07/2020 14:01, Luigi Marongiu wrote:

Hello,
Is it possible to calculate with a single function the odd ratios?
Now I can use this implement:
```
or <- (De/He)/(Dn/Hn) # Disease exposed, Healthy non-exposed
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
 spe = "")
```
for instance,
```
De <-6
Dn <-3
He <-4
Hn <-5
or <- (De/He)/(Dn/Hn)
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
 spe = "")

OR: 2.5 ( 0.37 - 16.889 )

```
Is there a simple function from some package that can also add a
p-value to this test? Or how can I calculate the p-value on my own?



--
Michael
http://www.dewey.myzen.co.uk/home.html

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 odd ratios with R?

2020-07-06 Thread Sorkin, John

Luigi,
Odds ratios can be produced using a logistic regression, which can be performed 
using the glm function. The following has a detailed description of how 
logistic regression can be performed using R:

https://stats.idre.ucla.edu/r/dae/logit-regression/

John


John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric 
Medicine
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)



From: R-help  on behalf of Luigi Marongiu 

Sent: Monday, July 6, 2020 9:01 AM
To: r-help 
Subject: [R] how to calculate odd ratios with R?

Hello,
Is it possible to calculate with a single function the odd ratios?
Now I can use this implement:
```
or <- (De/He)/(Dn/Hn) # Disease exposed, Healthy non-exposed
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
spe = "")
```
for instance,
```
De <-6
Dn <-3
He <-4
Hn <-5
or <- (De/He)/(Dn/Hn)
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
spe = "")
> OR: 2.5 ( 0.37 - 16.889 )
```
Is there a simple function from some package that can also add a
p-value to this test? Or how can I calculate the p-value on my own?
--
Best regards,
Luigi

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://nam03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-helpdata=02%7C01%7C%7C2cd90411f0e34701d7f308d821acceb7%7C717009a620de461a88940312a395cac9%7C0%7C0%7C637296373389265031sdata=GVQUX4DafNEu29kEupPbDrNblkQyas3LquN%2FVahCFPw%3Dreserved=0
PLEASE do read the posting guide 
https://nam03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.r-project.org%2Fposting-guide.htmldata=02%7C01%7C%7C2cd90411f0e34701d7f308d821acceb7%7C717009a620de461a88940312a395cac9%7C0%7C0%7C637296373389265031sdata=9hLF%2F0BrW3Tn%2BleeHoXaoRXY0NNxeXVQZAJQEvLBn7E%3Dreserved=0
and provide commented, minimal, self-contained, reproducible code.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate odd ratios with R?

2020-07-06 Thread Luigi Marongiu

Hello,
Is it possible to calculate with a single function the odd ratios?
Now I can use this implement:
```
or <- (De/He)/(Dn/Hn) # Disease exposed, Healthy non-exposed
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
spe = "")
```
for instance,
```
De <-6
Dn <-3
He <-4
Hn <-5
or <- (De/He)/(Dn/Hn)
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
spe = "")
> OR: 2.5 ( 0.37 - 16.889 )
```
Is there a simple function from some package that can also add a
p-value to this test? Or how can I calculate the p-value on my own?
-- 
Best regards,
Luigi

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate p value and correlation coefficient for Spearman’s correlation of differential expression data with 40000 permutations?

2019-10-31 Thread Ana Marija

Hello,

I have 3 groups,let's call them g1, g2, g3. Each of them is a result
of analysis in between groups of conditions, and g1 looks like this

   geneSymbol  logFC t  P.Value
adj.P.Val Beta
EXykpF1BRREdXnv9Xk  MKI67 -0.3115880 -5.521186 5.772137e-07
0.008986062 4.3106665
0Tm7hdRJxd9zoevPlA CCL3L3  0.1708020  4.162115 9.109798e-05
0.508784638 0.6630544
u_M5UdFdhg3lZ.qe64 UBE2G1 -0.1528149 -4.031466 1.430822e-04
0.508784638 0.3354065
lkkLCXcnzL9NXFXTl4 SEL1L3 -0.2138729 -3.977482 1.720517e-04
0.508784638 0.2015945
0Uu3XrB6Bd14qoNeuc  ZFP36  0.1667330  3.944917 1.921715e-04
0.508784638 0.1213335
3h7Sgq2i3sAUkxL_n8  ITGB5  0.3419488  3.938960 1.960886e-04
0.508784638 0.1066896

g2 and g2 look the same and each has  15568 entries (genes)

How to calculate p value and correlation coefficient for Spearman’s
correlation for this data for 4 permutations?

I joined all 3 groups, g1, g2, g3, and extracted only Beta (B)

I got this data frame (d), with matching 15568 entries:

 B.x   B.y B
EXykpF1BRREdXnv9Xk -4.970533 -4.752771 -5.404054
0Tm7hdRJxd9zoevPlA -4.862168 -5.147294 -3.909654
u_M5UdFdhg3lZ.qe64 -5.368846 -5.396183 -5.405330
lkkLCXcnzL9NXFXTl4 -4.367704 -4.847795 -5.148524
0Uu3XrB6Bd14qoNeuc -5.286592 -4.949305 -5.278798
3h7Sgq2i3sAUkxL_n8 -4.579528 -2.403240 -4.710600

To calculate Spearman’s I could use in R:

> cor(d,use="pairwise.complete.obs",method="spearman")
B.x  B.yB
B.x 1.0  0.234171932  0.002474729
B.y 0.234171932  1.0 -0.005469126
B   0.002474729 -0.005469126  1.0

Can someone please tell me what would be the method to use to get
correlation coefficient and p value taken in account number of
permutations? And am I am correct to use Beta in order to do
correlation in between these 3 groups?

Thanks!

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 multiple meta p values

2019-10-31 Thread Ana Marija

Can you please get back to me about this, I need this meta p values
for manuscript I have to submit next week

On Wed, Oct 30, 2019 at 5:35 PM Ana Marija  wrote:
>
> I also tried to do it this way:
>
> d$META <- sapply(seq_len(nrow(d)), function(rn) {
>   unlist(sumz(as.matrix(d[,.(LCL,Retina)])[rn,], weights =
> as.vector(d[,.(wl,wr)])[rn,],
>   na.action=na.fail)["p"])
> })
>
> but again I am getting error:
> Error in sumz(as.matrix(d[, .(LCL, Retina)])[rn, ], weights = as.vector(d[,  :
>   Must have at least two valid p values
>
> for this reference these are details about my data frame:
> > head(d)
> chrpos gene_id LCL  Retina
>wl   wr
> 1: chr1 775930 ENSG0237094 0.3559520 9.72251e-05 31.62278 21.2838
> 2: chr1 815963 ENSG0237094 0.2648080 3.85837e-06 31.62278 21.2838
> 3: chr1 816376 ENSG0237094 0.3313120 3.85824e-06 31.62278 21.2838
> 4: chr1 817186 ENSG0237094 0.0912854 3.75134e-06 31.62278 21.2838
> 5: chr1 817341 ENSG0237094 0.1020520 3.75134e-06 31.62278 21.2838
> 6: chr1 817514 ENSG0237094 0.0831412 3.82866e-06 31.62278 21.2838
> > sapply(d,class)
> chr pos gene_id LCL  Retina  wl
> "character" "character" "character"   "numeric"   "numeric"   "numeric"
>  wr
>   "numeric"
> > sum(is.na(d$LCL))
> [1] 0
> > sum(is.na(d$Retina))
> [1] 0
> > sum(is.na(d$wl))
> [1] 0
> > sum(is.na(d$wr))
> [1] 0
> > dim(d)
> [1] 1668837   7
>
> On Wed, Oct 30, 2019 at 4:52 PM Ana Marija  
> wrote:
> >
> > Hi Michael,
> >
> > this still doesn't work, by data frame has a few less columns now, but
> > the principle is still the same:
> >
> > > head(d)
> > chrpos gene_id LCL
> > Retina   wl   wr
> > 1: chr1 775930 ENSG0237094 0.3559520 9.72251e-05 31.62278 21.2838
> > 2: chr1 815963 ENSG0237094 0.2648080 3.85837e-06 31.62278 21.2838
> > 3: chr1 816376 ENSG0237094 0.3313120 3.85824e-06 31.62278 21.2838
> > 4: chr1 817186 ENSG0237094 0.0912854 3.75134e-06 31.62278 21.2838
> > 5: chr1 817341 ENSG0237094 0.1020520 3.75134e-06 31.62278 21.2838
> > 6: chr1 817514 ENSG0237094 0.0831412 3.82866e-06 31.62278 21.2838
> >
> > so solution for the first row should be:
> > > sumz(c(0.3559520,9.72251e-05), weights = c(31.62278,21.2838), na.action = 
> > > na.fail)
> > sumz =  2.386896 p =  0.008495647
> >
> > when I run what you proposed in the last email:
> >
> > helper <- function(x) {
> >   p <- sumz(as.numeric(x[4:5]), weights = as.numeric(x[6:7]))$p
> >   p
> > }
> >
> > d$META <- apply(d, MARGIN = 1, helper)
> >
> > I am getting:
> >
> > Error in sumz(as.numeric(x[4:5]), weights = as.numeric(x[6:7])) :
> >   Must have at least two valid p values
> >
> > Please advise,
> > Ana
> >
> > On Wed, Oct 30, 2019 at 5:02 AM Michael Dewey  
> > wrote:
> > >
> > > Dear Ana
> > >
> > > Yes, when apply coerces q to a matrix it does so as a character matrix
> > > because of the values in the first column. So you need to wrap the
> > > references to x in helper in as.numeric() tat is to day like
> > > as.numeric(x[2:4]) and similarly for the other one. Sorry about that, I
> > > should have thought of it before.
> > >
> > > When I next update metap I will try to get it to degrade more gracefully
> > > when it finds an error.
> > >
> > > Michael
> > >
> > > On 28/10/2019 19:06, Ana Marija wrote:
> > > > Hi Michael,
> > > >
> > > > I tried what you proposed with my data frame q:
> > > >
> > > >> head(q)
> > > > IDP G  E
> > > >   wb  wg   we
> > > > 1:  rs1029830 0.0979931 0.0054060 0.39160 580.6436 40.6325 35.39774
> > > > 2:  rs1029832 0.1501820 0.0028140 0.39320 580.6436 40.6325 35.39774
> > > > 3: rs11078374 0.1701250 0.0009805 0.49730 580.6436 40.6325 35.39774
> > > > 4:  rs1124961 0.1710150 0.7252000 0.05737 580.6436 40.6325 35.39774
> > > > 5:  rs1135237 0.1493650 0.6851000 0.06354 580.6436 40.6325 35.39774
> > > > 6: rs11867934 0.0757972 0.0006140 0.00327 580.6436 40.6325 35.39774
> > > >
> > > > so the solution of the first row would be this:
> > > >> sumz(c(0.0979931,0.0054060,0.39160), weights = 
> > > >> c(580.6436,40.6325,35.39774), na.action = na.fail)
> > > > sumz =  1.481833 p =  0.06919239
> > > >
> > > > I tried applying the function you wrote:
> > > > helper <- function(x) {
> > > >p <- sumz(x[2:4], weights = x[5:7])$p
> > > >p
> > > > }
> > > >
> > > > With:
> > > >
> > > > q$META <- apply(q, MARGIN = 1, helper)
> > > >
> > > > # I want to make a new column in q named META with results
> > > > but I got this error:
> > > >   Error in sumz(x[2:4], weights = x[5:7]) :
> > > >Must have at least two valid p values
> > > >
> > > > Please advise,
> > > > Ana
> > > >
> > > > On Sun, Oct 27, 2019 at 9:49 AM Michael Dewey  
> > > > wrote:
> > > >>
> > > >> Dear Ana
> > > >>
> > > >> There must be several ways of doing this but see

Re: [R] how to calculate multiple meta p values

2019-10-30 Thread Ana Marija

I also tried to do it this way:

d$META <- sapply(seq_len(nrow(d)), function(rn) {
  unlist(sumz(as.matrix(d[,.(LCL,Retina)])[rn,], weights =
as.vector(d[,.(wl,wr)])[rn,],
  na.action=na.fail)["p"])
})

but again I am getting error:
Error in sumz(as.matrix(d[, .(LCL, Retina)])[rn, ], weights = as.vector(d[,  :
  Must have at least two valid p values

for this reference these are details about my data frame:
> head(d)
chrpos gene_id LCL  Retina
   wl   wr
1: chr1 775930 ENSG0237094 0.3559520 9.72251e-05 31.62278 21.2838
2: chr1 815963 ENSG0237094 0.2648080 3.85837e-06 31.62278 21.2838
3: chr1 816376 ENSG0237094 0.3313120 3.85824e-06 31.62278 21.2838
4: chr1 817186 ENSG0237094 0.0912854 3.75134e-06 31.62278 21.2838
5: chr1 817341 ENSG0237094 0.1020520 3.75134e-06 31.62278 21.2838
6: chr1 817514 ENSG0237094 0.0831412 3.82866e-06 31.62278 21.2838
> sapply(d,class)
chr pos gene_id LCL  Retina  wl
"character" "character" "character"   "numeric"   "numeric"   "numeric"
 wr
  "numeric"
> sum(is.na(d$LCL))
[1] 0
> sum(is.na(d$Retina))
[1] 0
> sum(is.na(d$wl))
[1] 0
> sum(is.na(d$wr))
[1] 0
> dim(d)
[1] 1668837   7

On Wed, Oct 30, 2019 at 4:52 PM Ana Marija  wrote:
>
> Hi Michael,
>
> this still doesn't work, by data frame has a few less columns now, but
> the principle is still the same:
>
> > head(d)
> chrpos gene_id LCL
> Retina   wl   wr
> 1: chr1 775930 ENSG0237094 0.3559520 9.72251e-05 31.62278 21.2838
> 2: chr1 815963 ENSG0237094 0.2648080 3.85837e-06 31.62278 21.2838
> 3: chr1 816376 ENSG0237094 0.3313120 3.85824e-06 31.62278 21.2838
> 4: chr1 817186 ENSG0237094 0.0912854 3.75134e-06 31.62278 21.2838
> 5: chr1 817341 ENSG0237094 0.1020520 3.75134e-06 31.62278 21.2838
> 6: chr1 817514 ENSG0237094 0.0831412 3.82866e-06 31.62278 21.2838
>
> so solution for the first row should be:
> > sumz(c(0.3559520,9.72251e-05), weights = c(31.62278,21.2838), na.action = 
> > na.fail)
> sumz =  2.386896 p =  0.008495647
>
> when I run what you proposed in the last email:
>
> helper <- function(x) {
>   p <- sumz(as.numeric(x[4:5]), weights = as.numeric(x[6:7]))$p
>   p
> }
>
> d$META <- apply(d, MARGIN = 1, helper)
>
> I am getting:
>
> Error in sumz(as.numeric(x[4:5]), weights = as.numeric(x[6:7])) :
>   Must have at least two valid p values
>
> Please advise,
> Ana
>
> On Wed, Oct 30, 2019 at 5:02 AM Michael Dewey  wrote:
> >
> > Dear Ana
> >
> > Yes, when apply coerces q to a matrix it does so as a character matrix
> > because of the values in the first column. So you need to wrap the
> > references to x in helper in as.numeric() tat is to day like
> > as.numeric(x[2:4]) and similarly for the other one. Sorry about that, I
> > should have thought of it before.
> >
> > When I next update metap I will try to get it to degrade more gracefully
> > when it finds an error.
> >
> > Michael
> >
> > On 28/10/2019 19:06, Ana Marija wrote:
> > > Hi Michael,
> > >
> > > I tried what you proposed with my data frame q:
> > >
> > >> head(q)
> > > IDP G  E
> > >   wb  wg   we
> > > 1:  rs1029830 0.0979931 0.0054060 0.39160 580.6436 40.6325 35.39774
> > > 2:  rs1029832 0.1501820 0.0028140 0.39320 580.6436 40.6325 35.39774
> > > 3: rs11078374 0.1701250 0.0009805 0.49730 580.6436 40.6325 35.39774
> > > 4:  rs1124961 0.1710150 0.7252000 0.05737 580.6436 40.6325 35.39774
> > > 5:  rs1135237 0.1493650 0.6851000 0.06354 580.6436 40.6325 35.39774
> > > 6: rs11867934 0.0757972 0.0006140 0.00327 580.6436 40.6325 35.39774
> > >
> > > so the solution of the first row would be this:
> > >> sumz(c(0.0979931,0.0054060,0.39160), weights = 
> > >> c(580.6436,40.6325,35.39774), na.action = na.fail)
> > > sumz =  1.481833 p =  0.06919239
> > >
> > > I tried applying the function you wrote:
> > > helper <- function(x) {
> > >p <- sumz(x[2:4], weights = x[5:7])$p
> > >p
> > > }
> > >
> > > With:
> > >
> > > q$META <- apply(q, MARGIN = 1, helper)
> > >
> > > # I want to make a new column in q named META with results
> > > but I got this error:
> > >   Error in sumz(x[2:4], weights = x[5:7]) :
> > >Must have at least two valid p values
> > >
> > > Please advise,
> > > Ana
> > >
> > > On Sun, Oct 27, 2019 at 9:49 AM Michael Dewey  
> > > wrote:
> > >>
> > >> Dear Ana
> > >>
> > >> There must be several ways of doing this but see below for an idea with
> > >> comments in-line.
> > >>
> > >> On 26/10/2019 00:31, Ana Marija wrote:
> > >>> Hello,
> > >>>
> > >>> I would like to use this package metap
> > >>> to calculate multiple o values
> > >>>
> > >>> I have my data frame with 3 p values
> >  head(tt)
> > >>> RSG   E  B
> > >>> 1: rs2089177   0.9986   0.7153   0.604716
> > >>> 2: rs4360974   0.9738   0.7838   0.430228
> >

Re: [R] how to calculate multiple meta p values

2019-10-30 Thread Ana Marija

Hi Michael,

this still doesn't work, by data frame has a few less columns now, but
the principle is still the same:

> head(d)
chrpos gene_id LCL
Retina   wl   wr
1: chr1 775930 ENSG0237094 0.3559520 9.72251e-05 31.62278 21.2838
2: chr1 815963 ENSG0237094 0.2648080 3.85837e-06 31.62278 21.2838
3: chr1 816376 ENSG0237094 0.3313120 3.85824e-06 31.62278 21.2838
4: chr1 817186 ENSG0237094 0.0912854 3.75134e-06 31.62278 21.2838
5: chr1 817341 ENSG0237094 0.1020520 3.75134e-06 31.62278 21.2838
6: chr1 817514 ENSG0237094 0.0831412 3.82866e-06 31.62278 21.2838

so solution for the first row should be:
> sumz(c(0.3559520,9.72251e-05), weights = c(31.62278,21.2838), na.action = 
> na.fail)
sumz =  2.386896 p =  0.008495647

when I run what you proposed in the last email:

helper <- function(x) {
  p <- sumz(as.numeric(x[4:5]), weights = as.numeric(x[6:7]))$p
  p
}

d$META <- apply(d, MARGIN = 1, helper)

I am getting:

Error in sumz(as.numeric(x[4:5]), weights = as.numeric(x[6:7])) :
  Must have at least two valid p values

Please advise,
Ana

On Wed, Oct 30, 2019 at 5:02 AM Michael Dewey  wrote:
>
> Dear Ana
>
> Yes, when apply coerces q to a matrix it does so as a character matrix
> because of the values in the first column. So you need to wrap the
> references to x in helper in as.numeric() tat is to day like
> as.numeric(x[2:4]) and similarly for the other one. Sorry about that, I
> should have thought of it before.
>
> When I next update metap I will try to get it to degrade more gracefully
> when it finds an error.
>
> Michael
>
> On 28/10/2019 19:06, Ana Marija wrote:
> > Hi Michael,
> >
> > I tried what you proposed with my data frame q:
> >
> >> head(q)
> > IDP G  E
> >   wb  wg   we
> > 1:  rs1029830 0.0979931 0.0054060 0.39160 580.6436 40.6325 35.39774
> > 2:  rs1029832 0.1501820 0.0028140 0.39320 580.6436 40.6325 35.39774
> > 3: rs11078374 0.1701250 0.0009805 0.49730 580.6436 40.6325 35.39774
> > 4:  rs1124961 0.1710150 0.7252000 0.05737 580.6436 40.6325 35.39774
> > 5:  rs1135237 0.1493650 0.6851000 0.06354 580.6436 40.6325 35.39774
> > 6: rs11867934 0.0757972 0.0006140 0.00327 580.6436 40.6325 35.39774
> >
> > so the solution of the first row would be this:
> >> sumz(c(0.0979931,0.0054060,0.39160), weights = 
> >> c(580.6436,40.6325,35.39774), na.action = na.fail)
> > sumz =  1.481833 p =  0.06919239
> >
> > I tried applying the function you wrote:
> > helper <- function(x) {
> >p <- sumz(x[2:4], weights = x[5:7])$p
> >p
> > }
> >
> > With:
> >
> > q$META <- apply(q, MARGIN = 1, helper)
> >
> > # I want to make a new column in q named META with results
> > but I got this error:
> >   Error in sumz(x[2:4], weights = x[5:7]) :
> >Must have at least two valid p values
> >
> > Please advise,
> > Ana
> >
> > On Sun, Oct 27, 2019 at 9:49 AM Michael Dewey  
> > wrote:
> >>
> >> Dear Ana
> >>
> >> There must be several ways of doing this but see below for an idea with
> >> comments in-line.
> >>
> >> On 26/10/2019 00:31, Ana Marija wrote:
> >>> Hello,
> >>>
> >>> I would like to use this package metap
> >>> to calculate multiple o values
> >>>
> >>> I have my data frame with 3 p values
>  head(tt)
> >>> RSG   E  B
> >>> 1: rs2089177   0.9986   0.7153   0.604716
> >>> 2: rs4360974   0.9738   0.7838   0.430228
> >>> 3: rs6502526   0.9744   0.7839   0.429160
> >>> 4: rs8069906   0.7184   0.4918   0.521452
> >>> 5: rs9905280   0.7205   0.4861   0.465758
> >>> 6: rs4313843   0.9804   0.8522   0.474313
> >>>
> >>> and data frame with corresponding weights for each of the p values
> >>> from the tt data frame
> >>>
>  head(df)
> >>>  wg   we wbRS
> >>> 1 40.6325 35.39774 580.6436 rs2089177
> >>> 2 40.6325 35.39774 580.6436 rs4360974
> >>> 3 40.6325 35.39774 580.6436 rs6502526
> >>> 4 40.6325 35.39774 580.6436 rs8069906
> >>> 5 40.6325 35.39774 580.6436 rs9905280
> >>> 6 40.6325 35.39774 580.6436 rs4313843
> >>>
> >>> RS column is the same in df and tt
> >>>
> >>
> >> So you can create a new data-frame with merge()
> >>
> >> newdata <- merge(tt, df)
> >>
> >> which will use RS as the key to merge them on.
> >>
> >> The write a function of one argument, a seven element vector, which
> >> picks out the p-values and the weights and feeds them to sumz().
> >> Something like
> >>
> >> helper <- function(x) {
> >>p <- sumz(x[2:4], weights = x[5:7])$p
> >>p
> >> }
> >> Note you need to check that 2:4 and 5:7 are actually where they are in
> >> the row of newdat.
> >>
> >> Then use apply() to apply that to the rows of newdat.
> >>
> >> I have not tested any of this but the general idea should be OK even if
> >> the details are wrong.
> >>
> >> Michael
> >>
> >>
> >>> How to use this sunz() function to create a new data frame which would
> >>> look the same as tt only it would have

Re: [R] how to calculate multiple meta p values

2019-10-30 Thread Michael Dewey


Dear Ana

Yes, when apply coerces q to a matrix it does so as a character matrix 
because of the values in the first column. So you need to wrap the 
references to x in helper in as.numeric() tat is to day like 
as.numeric(x[2:4]) and similarly for the other one. Sorry about that, I 
should have thought of it before.


When I next update metap I will try to get it to degrade more gracefully 
when it finds an error.


Michael

On 28/10/2019 19:06, Ana Marija wrote:

Hi Michael,

I tried what you proposed with my data frame q:


head(q)

IDP G  E
  wb  wg   we
1:  rs1029830 0.0979931 0.0054060 0.39160 580.6436 40.6325 35.39774
2:  rs1029832 0.1501820 0.0028140 0.39320 580.6436 40.6325 35.39774
3: rs11078374 0.1701250 0.0009805 0.49730 580.6436 40.6325 35.39774
4:  rs1124961 0.1710150 0.7252000 0.05737 580.6436 40.6325 35.39774
5:  rs1135237 0.1493650 0.6851000 0.06354 580.6436 40.6325 35.39774
6: rs11867934 0.0757972 0.0006140 0.00327 580.6436 40.6325 35.39774

so the solution of the first row would be this:

sumz(c(0.0979931,0.0054060,0.39160), weights = c(580.6436,40.6325,35.39774), 
na.action = na.fail)

sumz =  1.481833 p =  0.06919239

I tried applying the function you wrote:
helper <- function(x) {
   p <- sumz(x[2:4], weights = x[5:7])$p
   p
}

With:

q$META <- apply(q, MARGIN = 1, helper)

# I want to make a new column in q named META with results
but I got this error:
  Error in sumz(x[2:4], weights = x[5:7]) :
   Must have at least two valid p values

Please advise,
Ana

On Sun, Oct 27, 2019 at 9:49 AM Michael Dewey  wrote:


Dear Ana

There must be several ways of doing this but see below for an idea with
comments in-line.

On 26/10/2019 00:31, Ana Marija wrote:

Hello,

I would like to use this package metap
to calculate multiple o values

I have my data frame with 3 p values

head(tt)

RSG   E  B
1: rs2089177   0.9986   0.7153   0.604716
2: rs4360974   0.9738   0.7838   0.430228
3: rs6502526   0.9744   0.7839   0.429160
4: rs8069906   0.7184   0.4918   0.521452
5: rs9905280   0.7205   0.4861   0.465758
6: rs4313843   0.9804   0.8522   0.474313

and data frame with corresponding weights for each of the p values
from the tt data frame


head(df)

 wg   we wbRS
1 40.6325 35.39774 580.6436 rs2089177
2 40.6325 35.39774 580.6436 rs4360974
3 40.6325 35.39774 580.6436 rs6502526
4 40.6325 35.39774 580.6436 rs8069906
5 40.6325 35.39774 580.6436 rs9905280
6 40.6325 35.39774 580.6436 rs4313843

RS column is the same in df and tt



So you can create a new data-frame with merge()

newdata <- merge(tt, df)

which will use RS as the key to merge them on.

The write a function of one argument, a seven element vector, which
picks out the p-values and the weights and feeds them to sumz().
Something like

helper <- function(x) {
   p <- sumz(x[2:4], weights = x[5:7])$p
   p
}
Note you need to check that 2:4 and 5:7 are actually where they are in
the row of newdat.

Then use apply() to apply that to the rows of newdat.

I have not tested any of this but the general idea should be OK even if
the details are wrong.

Michael



How to use this sunz() function to create a new data frame which would
look the same as tt only it would have additional column, say named
"META" which has calculated meta p values for each row

This i s example of how much would be p value in the first row:


sumz(c(0.9986,0.7153,0.604716), weights = c(40.6325,35.39774,580.6436), 
na.action = na.fail)

p =  0.6940048

Thanks
Ana

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



--
Michael
http://www.dewey.myzen.co.uk/home.html




--
Michael
http://www.dewey.myzen.co.uk/home.html

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 multiple meta p values

2019-10-28 Thread Ana Marija

Hi Michael,

I tried what you proposed with my data frame q:

> head(q)
   IDP G  E
 wb  wg   we
1:  rs1029830 0.0979931 0.0054060 0.39160 580.6436 40.6325 35.39774
2:  rs1029832 0.1501820 0.0028140 0.39320 580.6436 40.6325 35.39774
3: rs11078374 0.1701250 0.0009805 0.49730 580.6436 40.6325 35.39774
4:  rs1124961 0.1710150 0.7252000 0.05737 580.6436 40.6325 35.39774
5:  rs1135237 0.1493650 0.6851000 0.06354 580.6436 40.6325 35.39774
6: rs11867934 0.0757972 0.0006140 0.00327 580.6436 40.6325 35.39774

so the solution of the first row would be this:
> sumz(c(0.0979931,0.0054060,0.39160), weights = c(580.6436,40.6325,35.39774), 
> na.action = na.fail)
sumz =  1.481833 p =  0.06919239

I tried applying the function you wrote:
helper <- function(x) {
  p <- sumz(x[2:4], weights = x[5:7])$p
  p
}

With:

q$META <- apply(q, MARGIN = 1, helper)

# I want to make a new column in q named META with results
but I got this error:
 Error in sumz(x[2:4], weights = x[5:7]) :
  Must have at least two valid p values

Please advise,
Ana

On Sun, Oct 27, 2019 at 9:49 AM Michael Dewey  wrote:
>
> Dear Ana
>
> There must be several ways of doing this but see below for an idea with
> comments in-line.
>
> On 26/10/2019 00:31, Ana Marija wrote:
> > Hello,
> >
> > I would like to use this package metap
> > to calculate multiple o values
> >
> > I have my data frame with 3 p values
> >> head(tt)
> >RSG   E  B
> > 1: rs2089177   0.9986   0.7153   0.604716
> > 2: rs4360974   0.9738   0.7838   0.430228
> > 3: rs6502526   0.9744   0.7839   0.429160
> > 4: rs8069906   0.7184   0.4918   0.521452
> > 5: rs9905280   0.7205   0.4861   0.465758
> > 6: rs4313843   0.9804   0.8522   0.474313
> >
> > and data frame with corresponding weights for each of the p values
> > from the tt data frame
> >
> >> head(df)
> > wg   we wbRS
> > 1 40.6325 35.39774 580.6436 rs2089177
> > 2 40.6325 35.39774 580.6436 rs4360974
> > 3 40.6325 35.39774 580.6436 rs6502526
> > 4 40.6325 35.39774 580.6436 rs8069906
> > 5 40.6325 35.39774 580.6436 rs9905280
> > 6 40.6325 35.39774 580.6436 rs4313843
> >
> > RS column is the same in df and tt
> >
>
> So you can create a new data-frame with merge()
>
> newdata <- merge(tt, df)
>
> which will use RS as the key to merge them on.
>
> The write a function of one argument, a seven element vector, which
> picks out the p-values and the weights and feeds them to sumz().
> Something like
>
> helper <- function(x) {
>   p <- sumz(x[2:4], weights = x[5:7])$p
>   p
> }
> Note you need to check that 2:4 and 5:7 are actually where they are in
> the row of newdat.
>
> Then use apply() to apply that to the rows of newdat.
>
> I have not tested any of this but the general idea should be OK even if
> the details are wrong.
>
> Michael
>
>
> > How to use this sunz() function to create a new data frame which would
> > look the same as tt only it would have additional column, say named
> > "META" which has calculated meta p values for each row
> >
> > This i s example of how much would be p value in the first row:
> >
> >> sumz(c(0.9986,0.7153,0.604716), weights = c(40.6325,35.39774,580.6436), 
> >> na.action = na.fail)
> > p =  0.6940048
> >
> > Thanks
> > Ana
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
> --
> Michael
> http://www.dewey.myzen.co.uk/home.html

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 multiple meta p values

2019-10-27 Thread Michael Dewey


Dear Ana

There must be several ways of doing this but see below for an idea with 
comments in-line.


On 26/10/2019 00:31, Ana Marija wrote:

Hello,

I would like to use this package metap
to calculate multiple o values

I have my data frame with 3 p values

head(tt)

   RSG   E  B
1: rs2089177   0.9986   0.7153   0.604716
2: rs4360974   0.9738   0.7838   0.430228
3: rs6502526   0.9744   0.7839   0.429160
4: rs8069906   0.7184   0.4918   0.521452
5: rs9905280   0.7205   0.4861   0.465758
6: rs4313843   0.9804   0.8522   0.474313

and data frame with corresponding weights for each of the p values
from the tt data frame


head(df)

wg   we wbRS
1 40.6325 35.39774 580.6436 rs2089177
2 40.6325 35.39774 580.6436 rs4360974
3 40.6325 35.39774 580.6436 rs6502526
4 40.6325 35.39774 580.6436 rs8069906
5 40.6325 35.39774 580.6436 rs9905280
6 40.6325 35.39774 580.6436 rs4313843

RS column is the same in df and tt



So you can create a new data-frame with merge()

newdata <- merge(tt, df)

which will use RS as the key to merge them on.

The write a function of one argument, a seven element vector, which 
picks out the p-values and the weights and feeds them to sumz(). 
Something like


helper <- function(x) {
 p <- sumz(x[2:4], weights = x[5:7])$p
 p
}
Note you need to check that 2:4 and 5:7 are actually where they are in 
the row of newdat.


Then use apply() to apply that to the rows of newdat.

I have not tested any of this but the general idea should be OK even if 
the details are wrong.


Michael



How to use this sunz() function to create a new data frame which would
look the same as tt only it would have additional column, say named
"META" which has calculated meta p values for each row

This i s example of how much would be p value in the first row:


sumz(c(0.9986,0.7153,0.604716), weights = c(40.6325,35.39774,580.6436), 
na.action = na.fail)

p =  0.6940048

Thanks
Ana

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



--
Michael
http://www.dewey.myzen.co.uk/home.html

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 multiple meta p values

2019-10-25 Thread Ana Marija

this is the function I was referring to:
https://www.rdocumentation.org/packages/metap/versions/1.1/topics/sumz

On Fri, Oct 25, 2019 at 6:31 PM Ana Marija  wrote:
>
> Hello,
>
> I would like to use this package metap
> to calculate multiple o values
>
> I have my data frame with 3 p values
> > head(tt)
>   RSG   E  B
> 1: rs2089177   0.9986   0.7153   0.604716
> 2: rs4360974   0.9738   0.7838   0.430228
> 3: rs6502526   0.9744   0.7839   0.429160
> 4: rs8069906   0.7184   0.4918   0.521452
> 5: rs9905280   0.7205   0.4861   0.465758
> 6: rs4313843   0.9804   0.8522   0.474313
>
> and data frame with corresponding weights for each of the p values
> from the tt data frame
>
> > head(df)
>wg   we wbRS
> 1 40.6325 35.39774 580.6436 rs2089177
> 2 40.6325 35.39774 580.6436 rs4360974
> 3 40.6325 35.39774 580.6436 rs6502526
> 4 40.6325 35.39774 580.6436 rs8069906
> 5 40.6325 35.39774 580.6436 rs9905280
> 6 40.6325 35.39774 580.6436 rs4313843
>
> RS column is the same in df and tt
>
> How to use this sunz() function to create a new data frame which would
> look the same as tt only it would have additional column, say named
> "META" which has calculated meta p values for each row
>
> This i s example of how much would be p value in the first row:
>
> > sumz(c(0.9986,0.7153,0.604716), weights = c(40.6325,35.39774,580.6436), 
> > na.action = na.fail)
> p =  0.6940048
>
> Thanks
> Ana

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate multiple meta p values

2019-10-25 Thread Ana Marija

Hello,

I would like to use this package metap
to calculate multiple o values

I have my data frame with 3 p values
> head(tt)
  RSG   E  B
1: rs2089177   0.9986   0.7153   0.604716
2: rs4360974   0.9738   0.7838   0.430228
3: rs6502526   0.9744   0.7839   0.429160
4: rs8069906   0.7184   0.4918   0.521452
5: rs9905280   0.7205   0.4861   0.465758
6: rs4313843   0.9804   0.8522   0.474313

and data frame with corresponding weights for each of the p values
from the tt data frame

> head(df)
   wg   we wbRS
1 40.6325 35.39774 580.6436 rs2089177
2 40.6325 35.39774 580.6436 rs4360974
3 40.6325 35.39774 580.6436 rs6502526
4 40.6325 35.39774 580.6436 rs8069906
5 40.6325 35.39774 580.6436 rs9905280
6 40.6325 35.39774 580.6436 rs4313843

RS column is the same in df and tt

How to use this sunz() function to create a new data frame which would
look the same as tt only it would have additional column, say named
"META" which has calculated meta p values for each row

This i s example of how much would be p value in the first row:

> sumz(c(0.9986,0.7153,0.604716), weights = c(40.6325,35.39774,580.6436), 
> na.action = na.fail)
p =  0.6940048

Thanks
Ana

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 True Positive Rate in R?

2019-09-24 Thread Jim Lemon

Hi Ana,
Your minimum value in pvalR is very small and may be causing trouble.
As the qvalue package seems to be in Bioconductor, perhaps posting to
that help list would get an answer.

Jim

On Wed, Sep 25, 2019 at 1:48 AM Ana Marija  wrote:
>
> Hello,
>
> I tried using qvalue function:
>
> library(qvalue)
> qval_obj=qvalue(pvalR)
> pi1=1-qval_obj$pi0
>
> but after running:
>
> qval_obj=qvalue(pvalR)
> Error in smooth.spline(lambda, pi0, df = smooth.df) :
>   missing or infinite values in inputs are not allowed
>
> or
>
>  qval_obj=qvalue(pvalR,lambda=0.5)
> Error in pi0est(p, ...) :
>   ERROR: The estimated pi0 <= 0. Check that you have valid
> p-values or use a different range of lambda.
>
> I checked:
>
> max(pvalR)
> [1] 0.000352731
> min(pvalR)
> [1] 1.84872e-127
> sum(is.na(pvalR))
>[1] 0
>sum(is.infinite(pvalR))
>[1] 0
>
> Please advise,
>
> Thanks
> Ana
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 True Positive Rate in R?

2019-09-24 Thread Eric Berger

https://cran.r-project.org/web/packages/qvalue/index.html
states
 Package ‘qvalue’ was removed from the CRAN repository.
Just an FYI

On Tue, Sep 24, 2019 at 6:48 PM Ana Marija 
wrote:

> Hello,
>
> I tried using qvalue function:
>
> library(qvalue)
> qval_obj=qvalue(pvalR)
> pi1=1-qval_obj$pi0
>
> but after running:
>
> qval_obj=qvalue(pvalR)
> Error in smooth.spline(lambda, pi0, df = smooth.df) :
>   missing or infinite values in inputs are not allowed
>
> or
>
>  qval_obj=qvalue(pvalR,lambda=0.5)
> Error in pi0est(p, ...) :
>   ERROR: The estimated pi0 <= 0. Check that you have valid
> p-values or use a different range of lambda.
>
> I checked:
>
> max(pvalR)
> [1] 0.000352731
> min(pvalR)
> [1] 1.84872e-127
> sum(is.na(pvalR))
>[1] 0
>sum(is.infinite(pvalR))
>[1] 0
>
> Please advise,
>
> Thanks
> Ana
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate True Positive Rate in R?

2019-09-24 Thread Ana Marija

Hello,

I tried using qvalue function:

library(qvalue)
qval_obj=qvalue(pvalR)
pi1=1-qval_obj$pi0

but after running:

qval_obj=qvalue(pvalR)
Error in smooth.spline(lambda, pi0, df = smooth.df) :
  missing or infinite values in inputs are not allowed

or

 qval_obj=qvalue(pvalR,lambda=0.5)
Error in pi0est(p, ...) :
  ERROR: The estimated pi0 <= 0. Check that you have valid
p-values or use a different range of lambda.

I checked:

max(pvalR)
[1] 0.000352731
min(pvalR)
[1] 1.84872e-127
sum(is.na(pvalR))
   [1] 0
   sum(is.infinite(pvalR))
   [1] 0

Please advise,

Thanks
Ana

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate the monthly average from daily data

2019-09-12 Thread Subhamitra Patra

Dear R-users,

I have daily data from 03-01-1994 to 03-08-2017. In my datafile, The first
column is date and the second and third columns are the returns and volume
data. I want to estimate the average returns, and volumes for each month,
then, I want to export the results into excel.

So please help me in this regard.

Please find the attached datasheet.

Thank you.

-- 
*Best Regards,*
*Subhamitra Patra*
*Phd. Research Scholar*
*Department of Humanities and Social Sciences*
*Indian Institute of Technology, Kharagpur*
*INDIA*

[image: Mailtrack]

Sender
notified by
Mailtrack

09/12/19,
03:40:34 PM
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 Rolling mean for a List object

2017-05-14 Thread Gabor Grothendieck

Try this code which does not use rollapply:

w <- 3
Mean <- function(L) Reduce("+", L) / length(L)
lapply(w:length(Data), function(i) Mean(Data[seq(to = i, length = w)]))

On Sun, May 14, 2017 at 6:44 PM, Christofer Bogaso
 wrote:
> Hi again,
>
> I am looking to find a way on how to calculate Rolling average for the
> elements of a list. For example consider below object 'Data'. This is
> a list, where each elements are a Matrix. Basically, I am trying to
> get Rolling average of those Matrices with rolling window as 5.
>
> Data = structure(list(`2017-03-01` = structure(c(1.24915216491479e-06,
> -2.0209685810767e-06, -6.64165527006046e-06, -2.0209685810767e-06,
> 3.26966891657893e-06, 1.07453495291747e-05, -6.64165527006046e-06,
> 1.07453495291747e-05, 3.53132196103035e-05), .Dim = c(3L, 3L)),
> `2017-03-02` = structure(c(0.00863066441403338, -7.25585852047094e-05,
> -0.000950715788640005, -7.25585852047094e-05, 6.10004981580403e-07,
> 7.99273256915577e-06, -0.000950715788640005, 7.99273256915577e-06,
> 0.000104726642980084), .Dim = c(3L, 3L)), `2017-03-03` =
> structure(c(0.000785677680557358,
> 0.000283148300122928, 0.000170319078518317, 0.000283148300122928,
> 0.000102043066573597, 6.13808419844048e-05, 0.000170319078518317,
> 6.13808419844048e-05, 3.6921741860797e-05), .Dim = c(3L,
> 3L)), `2017-03-06` = structure(c(0.000100715163251975, 
> 1.80035062425799e-06,
> -5.05489732985851e-07, 1.80035062425799e-06, 3.21824665284709e-08,
> -9.03596565752718e-09, -5.05489732985851e-07, -9.03596565752718e-09,
> 2.53705461922188e-09), .Dim = c(3L, 3L)), `2017-03-07` =
> structure(c(0.000640065014281149,
> -0.000110994847091752, -0.000231235438845606, -0.000110994847091752,
> 1.92478198402357e-05, 4.00989612058198e-05, -0.000231235438845606,
> 4.00989612058198e-05, 8.35381203238728e-05), .Dim = c(3L,
> 3L)), `2017-03-08` = structure(c(7.72648041923266e-06,
> -2.11571338014623e-05,
> 7.82052544997182e-06, -2.11571338014623e-05, 5.79337921544145e-05,
> -2.14146538093767e-05, 7.82052544997182e-06, -2.14146538093767e-05,
> 7.91571517626794e-06), .Dim = c(3L, 3L)), `2017-03-09` =
> structure(c(4.43321118550061e-05,
> 1.90242249279913e-05, 5.68672547310199e-05, 1.90242249279913e-05,
> 8.16385953582618e-06, 2.44034267661023e-05, 5.68672547310199e-05,
> 2.44034267661023e-05, 7.29467766214148e-05), .Dim = c(3L,
> 3L)), `2017-03-10` = structure(c(0.000100081081692311, 
> 1.39245218598852e-05,
> 2.0935583168872e-05, 1.39245218598852e-05, 1.93735225227204e-06,
> 2.91281809264057e-06, 2.0935583168872e-05, 2.91281809264057e-06,
> 4.3794355057858e-06), .Dim = c(3L, 3L)), `2017-03-14` =
> structure(c(7.82185299651879e-06,
> -3.05963602958646e-05, -4.65590052688468e-05, -3.05963602958646e-05,
> 0.00011968228804236, 0.000182122586662866, -4.65590052688468e-05,
> 0.000182122586662866, 0.000277139058045361), .Dim = c(3L,
> 3L)), `2017-03-15` = structure(c(4.02156693772954e-05, 
> -2.2362610665311e-05,
> -2.08706726432905e-05, -2.2362610665311e-05, 1.24351120722764e-05,
> 1.16054944222453e-05, -2.08706726432905e-05, 1.16054944222453e-05,
> 1.08312253240602e-05), .Dim = c(3L, 3L)), `2017-03-16` =
> structure(c(2.64254966198469e-05,
> 5.78730550194069e-06, 5.0445603894268e-05, 5.78730550194069e-06,
> 1.26744656702641e-06, 1.10478196556107e-05, 5.0445603894268e-05,
> 1.10478196556107e-05, 9.62993804379875e-05), .Dim = c(3L,
> 3L)), `2017-03-17` = structure(c(0.000138433807049962, 
> 8.72005344938308e-05,
> 0.00014374477881467, 8.72005344938308e-05, 5.49282966209652e-05,
> 9.05459570205481e-05, 0.00014374477881467, 9.05459570205481e-05,
> 0.000149259504428865), .Dim = c(3L, 3L)), `2017-03-20` =
> structure(c(3.92058275846982e-05,
> 1.24332187386233e-05, -1.24235553811814e-05, 1.24332187386233e-05,
> 3.94290690251335e-06, -3.93984239286701e-06, -1.24235553811814e-05,
> -3.93984239286701e-06, 3.93678026502162e-06), .Dim = c(3L,
> 3L)), `2017-03-21` = structure(c(0.000407544227952838,
> -6.22427018306449e-05,
> 1.90596071859105e-05, -6.22427018306449e-05, 9.50609446890975e-06,
> -2.9109023406881e-06, 1.90596071859105e-05, -2.9109023406881e-06,
> 8.91360007491622e-07), .Dim = c(3L, 3L)), `2017-03-22` =
> structure(c(0.000220297355944482,
> 0.000282600064158173, 8.26030839524992e-05, 0.000282600064158173,
> 0.000362522718077154, 0.00010596421697645, 8.26030839524992e-05,
> 0.00010596421697645, 3.09729976068491e-05), .Dim = c(3L,
> 3L)), `2017-03-23` = structure(c(1.19559010537042e-05, 
> 3.56054556562106e-05,
> 5.51130473489473e-06, 3.56054556562106e-05, 0.000106035376739222,
> 1.64130261253175e-05, 5.51130473489473e-06, 1.64130261253175e-05,
> 2.54054292892148e-06), .Dim = c(3L, 3L)), `2017-03-24` =
> structure(c(0.000573948692221572,
> -7.36566239512158e-05, 5.40736580500709e-05,

[R] How to calculate Rolling mean for a List object

2017-05-14 Thread Christofer Bogaso

Hi again,

I am looking to find a way on how to calculate Rolling average for the
elements of a list. For example consider below object 'Data'. This is
a list, where each elements are a Matrix. Basically, I am trying to
get Rolling average of those Matrices with rolling window as 5.

Data = structure(list(`2017-03-01` = structure(c(1.24915216491479e-06,
-2.0209685810767e-06, -6.64165527006046e-06, -2.0209685810767e-06,
3.26966891657893e-06, 1.07453495291747e-05, -6.64165527006046e-06,
1.07453495291747e-05, 3.53132196103035e-05), .Dim = c(3L, 3L)),
`2017-03-02` = structure(c(0.00863066441403338, -7.25585852047094e-05,
-0.000950715788640005, -7.25585852047094e-05, 6.10004981580403e-07,
7.99273256915577e-06, -0.000950715788640005, 7.99273256915577e-06,
0.000104726642980084), .Dim = c(3L, 3L)), `2017-03-03` =
structure(c(0.000785677680557358,
0.000283148300122928, 0.000170319078518317, 0.000283148300122928,
0.000102043066573597, 6.13808419844048e-05, 0.000170319078518317,
6.13808419844048e-05, 3.6921741860797e-05), .Dim = c(3L,
3L)), `2017-03-06` = structure(c(0.000100715163251975, 1.80035062425799e-06,
-5.05489732985851e-07, 1.80035062425799e-06, 3.21824665284709e-08,
-9.03596565752718e-09, -5.05489732985851e-07, -9.03596565752718e-09,
2.53705461922188e-09), .Dim = c(3L, 3L)), `2017-03-07` =
structure(c(0.000640065014281149,
-0.000110994847091752, -0.000231235438845606, -0.000110994847091752,
1.92478198402357e-05, 4.00989612058198e-05, -0.000231235438845606,
4.00989612058198e-05, 8.35381203238728e-05), .Dim = c(3L,
3L)), `2017-03-08` = structure(c(7.72648041923266e-06,
-2.11571338014623e-05,
7.82052544997182e-06, -2.11571338014623e-05, 5.79337921544145e-05,
-2.14146538093767e-05, 7.82052544997182e-06, -2.14146538093767e-05,
7.91571517626794e-06), .Dim = c(3L, 3L)), `2017-03-09` =
structure(c(4.43321118550061e-05,
1.90242249279913e-05, 5.68672547310199e-05, 1.90242249279913e-05,
8.16385953582618e-06, 2.44034267661023e-05, 5.68672547310199e-05,
2.44034267661023e-05, 7.29467766214148e-05), .Dim = c(3L,
3L)), `2017-03-10` = structure(c(0.000100081081692311, 1.39245218598852e-05,
2.0935583168872e-05, 1.39245218598852e-05, 1.93735225227204e-06,
2.91281809264057e-06, 2.0935583168872e-05, 2.91281809264057e-06,
4.3794355057858e-06), .Dim = c(3L, 3L)), `2017-03-14` =
structure(c(7.82185299651879e-06,
-3.05963602958646e-05, -4.65590052688468e-05, -3.05963602958646e-05,
0.00011968228804236, 0.000182122586662866, -4.65590052688468e-05,
0.000182122586662866, 0.000277139058045361), .Dim = c(3L,
3L)), `2017-03-15` = structure(c(4.02156693772954e-05, -2.2362610665311e-05,
-2.08706726432905e-05, -2.2362610665311e-05, 1.24351120722764e-05,
1.16054944222453e-05, -2.08706726432905e-05, 1.16054944222453e-05,
1.08312253240602e-05), .Dim = c(3L, 3L)), `2017-03-16` =
structure(c(2.64254966198469e-05,
5.78730550194069e-06, 5.0445603894268e-05, 5.78730550194069e-06,
1.26744656702641e-06, 1.10478196556107e-05, 5.0445603894268e-05,
1.10478196556107e-05, 9.62993804379875e-05), .Dim = c(3L,
3L)), `2017-03-17` = structure(c(0.000138433807049962, 8.72005344938308e-05,
0.00014374477881467, 8.72005344938308e-05, 5.49282966209652e-05,
9.05459570205481e-05, 0.00014374477881467, 9.05459570205481e-05,
0.000149259504428865), .Dim = c(3L, 3L)), `2017-03-20` =
structure(c(3.92058275846982e-05,
1.24332187386233e-05, -1.24235553811814e-05, 1.24332187386233e-05,
3.94290690251335e-06, -3.93984239286701e-06, -1.24235553811814e-05,
-3.93984239286701e-06, 3.93678026502162e-06), .Dim = c(3L,
3L)), `2017-03-21` = structure(c(0.000407544227952838,
-6.22427018306449e-05,
1.90596071859105e-05, -6.22427018306449e-05, 9.50609446890975e-06,
-2.9109023406881e-06, 1.90596071859105e-05, -2.9109023406881e-06,
8.91360007491622e-07), .Dim = c(3L, 3L)), `2017-03-22` =
structure(c(0.000220297355944482,
0.000282600064158173, 8.26030839524992e-05, 0.000282600064158173,
0.000362522718077154, 0.00010596421697645, 8.26030839524992e-05,
0.00010596421697645, 3.09729976068491e-05), .Dim = c(3L,
3L)), `2017-03-23` = structure(c(1.19559010537042e-05, 3.56054556562106e-05,
5.51130473489473e-06, 3.56054556562106e-05, 0.000106035376739222,
1.64130261253175e-05, 5.51130473489473e-06, 1.64130261253175e-05,
2.54054292892148e-06), .Dim = c(3L, 3L)), `2017-03-24` =
structure(c(0.000573948692221572,
-7.36566239512158e-05, 5.40736580500709e-05, -7.36566239512158e-05,
9.45258404700116e-06, -6.93944101735685e-06, 5.40736580500709e-05,
-6.93944101735685e-06, 5.0944632064554e-06), .Dim = c(3L,
3L)), `2017-03-27` = structure(c(6.50931905856128e-06, -6.3937553506226e-07,
3.58314387213273e-06, -6.3937553506226e-07, 6.28024331206322e-08,
-3.51953024554351e-07, 3.58314387213273e-06, -3.51953024554351e-07,
1.97239064376729e-06), .Dim = c(3L, 3L)), `2017-03-28` =

Re: [R] How to calculate row means while ignore NAs

2016-10-28 Thread lily li

My apologize, it has been solved. Just include w inside of select, such as:
select = c(w, x, y)

On Fri, Oct 28, 2016 at 12:06 PM, lily li  wrote:

> Hi R users,
>
> I have the dataframe as below:
>
> w=c(5,6,7,8)
> x=c(1,2,3,4)
> y=c(1,2,3)
> length(y)=4
> z=data.frame(w,x,y)
>
> z$mean1 <- rowMeans(subset(z, select = c(x, y)), na.rm = T)
> z$mean2 <- rowMeans(subset(z, select = c(x, y)), na.rm=F)
>
>   w x  y mean1 mean2
> 1 5 1  1 1 1
> 2 6 2  2 2 2
> 3 7 3  3 3 3
> 4 8 4 NA 4NA
>
> How to calculate mean value for the three columns, while default removing
> NAs.
> For example, for the fourth row, the mean value should be (8+4)/2 = 6
>
> Thanks for your help.
>

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate row means while ignore NAs

2016-10-28 Thread lily li

Hi R users,

I have the dataframe as below:

w=c(5,6,7,8)
x=c(1,2,3,4)
y=c(1,2,3)
length(y)=4
z=data.frame(w,x,y)

z$mean1 <- rowMeans(subset(z, select = c(x, y)), na.rm = T)
z$mean2 <- rowMeans(subset(z, select = c(x, y)), na.rm=F)

  w x  y mean1 mean2
1 5 1  1 1 1
2 6 2  2 2 2
3 7 3  3 3 3
4 8 4 NA 4NA

How to calculate mean value for the three columns, while default removing
NAs.
For example, for the fourth row, the mean value should be (8+4)/2 = 6

Thanks for your help.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate spatial correlation between two data sets

2016-03-19 Thread Yogesh Tiwari via R-help

Dear R Users,

I am using R on Windows.

1) How to calculate spatial correlation between point observation at one
location and model simulated data over some particular region. For example,
observation is at only one location 18N, 72 E and model is at 0.5x0.5 grid.
So how to calculate correlation between every grid of model (over the
region 50:100E, 10S:40N) and this point observation at 18N, 72 E. Model
data is in netcdf format and observation is in text format. I would like to
save output fine in netcdf format. We have extracted model data as same day
as observation.

2) Its similar to above but, how to calculate spatial correlation between
 two model simulated data (i.e. grid to grid correlation over the
region50:100E, 10S:40N). One model has 0.25x0.25 and other is has 0.5x0.5
horizontal resolution.  I would like to save output fine in netcdf format

Great Thanks,

Best regards,
Yogesh


-- 
Yogesh K. Tiwari (Dr.rer.nat),
Scientist,
Centre for Climate Change Research,
Indian Institute of Tropical Meteorology,
Homi Bhabha Road,
Pashan,
Pune

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate the prediction interval without knowing the functional form

2016-01-03 Thread Wensui Liu

If I have predictions derived empirically without knowing the functional
form, is there a way to calculate the prediction interval?

Thanks


-- 
WenSui Liu
https://statcompute.wordpress.com/

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the prediction interval without knowing the functional form

2016-01-03 Thread Bert Gunter

Standard answer: bootstrap. However, "derived empirically" is too
vague to know whether the standard answer applies.

As this appears to be primarily a statistics, not an R issue, I
suggest that you post on a statistics list like
stats.stackexchange.com instead, perhaps with some more details on
what "derived empirically" means.

Cheers,
Bert

Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )

On Sun, Jan 3, 2016 at 2:57 PM, Wensui Liu  wrote:
> If I have predictions derived empirically without knowing the functional
> form, is there a way to calculate the prediction interval?
>
> Thanks
>
>
> --
> WenSui Liu
> https://statcompute.wordpress.com/
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the prediction interval without knowing the functional form

2016-01-03 Thread Duncan Murdoch


On 03/01/2016 5:57 PM, Wensui Liu wrote:

If I have predictions derived empirically without knowing the functional
form, is there a way to calculate the prediction interval?


As Bert said, this isn't really an R question, but I'd say the answer is 
"probably not".  The prediction interval depends on the uncertainty of 
individual observations, and that isn't usually reflected in 
predictions.  For example, if y_i is N(mu, sigma^2), and we fit the 
model to N observations, all of our predictions will be ybar, and there 
will be no indication of sigma.


Duncan Murdoch

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 variance on multiple numbers at once?

2015-11-12 Thread Luigi Marongiu

Essentially the problem is related to this problem. I am measuring the
fluorescence (flu) generated over time (clock) from 5 samples (samp)
which are grouped by two targets (targ).
I can calculate the variance of the fluorescence for each sample at a
given time (in this example, for sample 1); but if I consider a target
at the time there would be multiple readings at once to consider.
would the variance formula still be OK to use.
In this example:
>>>
samp <- c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4,
5, 5, 5, 5, 5)
clock <- rep(1:5,5)
targ <- c(rep("A", 5), rep("B", 5), rep("A", 5), rep("A", 5), rep("B", 5))
flu <- c(-0.012, -0.01, -0.008, -0.002, -0.001, -0.007, -0.008,
-0.009, -0.009, -0.012, -0.002, -0.003, -0.003, 0.001, 0.002, -0.006,
-0.001, 0.001, 0.002, 0.002, -0.002, -0.003, -0.003, -0.002, -0.001)
my.data <- data.frame(well, clock, targ, flu, stringsAsFactors = FALSE)

# variance with individual numbers
sub <- subset(my.data, samp == 1)
plot(sub$flu ~ sub$clock)
abline(lm(sub$flu ~ sub$clock))
for (i in 2:nrow(sub)) {
  X <- subset(sub, clock <= i)
  V <- var(X$flu)
  cat("variance at clock ", i, " = ", V, "\n", sep="")
}
# variance with multiple numbers
sub <- subset(my.data, targ == 'A')
plot(sub$flu ~ sub$clock)
abline(lm(sub$flu ~ sub$clock))
for (i in 2:max(sub$clock)) {
  X <- subset(sub, clock <= i)
  V <- var(X$flu)
  cat("variance at clock ", i, " = ", V, "\n", sep="")
}

the results for the individual numbers are:
variance at clock 2 = 2e-06
variance at clock 3 = 4e-06
variance at clock 4 = 1.87e-05
variance at clock 5 = 2.38e-05

while the results for the multiple numbers are:
variance at clock 2 = 2.026667e-05
variance at clock 3 = 1.91e-05
variance at clock 4 = 2.026515e-05
variance at clock 5 = 1.995238e-05

shall I accept these latter values?
Thanks

On Thu, Nov 12, 2015 at 10:59 AM, John Kane <jrkrid...@inbox.com> wrote:
> I still don't understand what you are doing. I think we need to see your 
> actual data and the code your are using.  If S.Ellison's post has not shown 
> up it is below my signature.
>
> Have a look at 
> http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example
>  and/or http://adv-r.had.co.nz/Reproducibility.html for suggestions on how to 
> pose a question here. In particular data should be supplied in dput() format 
> as it gives us a copy of exactly how the data is formatted on your machine.
>
>
> John Kane
> Kingston ON Canada
>
> ## ===
>> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Luigi
>> Marongiu
>> if I have a sample set of the following numbers x1=0.09, x2=0.94, x3=0.48,
>> x4=0.74, x5=0.04 I can calculate the variance easily.
> Not without concatenating them into a vector, you can't. You need them in a 
> vector, as in
> var( c(x1, x2, x3, x4, x5) )
>
>> But if each x is actually a subset of multiple values, what would be the 
>> formula
>> to calculate the variance? and it is possible to implement such mathematical
>> function in R?
> This is what R wants anyway, so the function you are looking for is var()
>
>> For instance if I have the following: x1=(0.77, 0.22, 0.44), x2=(0.26, 0.89, 
>> 0.58),
>> x3=(0.20, 0.25, 0.91), x4=(0.06, 0.13, 0.26) and x5=(0.65, 0.16, 0.72) how 
>> can i
>> calculate the variance for each x?
> var(x1)
> var(x2)
> 
>
> or, if you want to be a bit more slick about it and do it in one line
>
> lapply(list( x1, x2, x3, ...), var  )
>
> (or sapply() if you want a vector result)
>
> ## ===
>
>
>> -Original Message-
>> From: marongiu.lu...@gmail.com
>> Sent: Wed, 11 Nov 2015 23:20:26 +
>> To: jrkrid...@inbox.com
>> Subject: Re: [R] How to calculate variance on multiple numbers at once?
>>
>> Thank you for the reply. For clarification, let's say that I can
>> calculate the sum of variances of the individual x numbers in five
>> consecutive steps (although of course there are better
>> implementations) where each step the sum is incremented by (x -
>> mean)^2.
>> In the case I am handling, at each step i have to consider 3 values at
>> once and I don't know how to relate them neither with the mathematical
>> formula nor with the R implementation.
>> Besides I haven't seen Ellison's answer...
>> Best regards
>> L
>>
>> On Wed, Nov 11, 2015 at 1:04 PM, John Kane <jrkrid...@inbox.com> wrote:
>>> I really don't understand what you are looking for but if S. Ellison's
>>

[R] How to calculate variance on multiple numbers at once?

2015-11-11 Thread Luigi Marongiu

Dear all,

if I have a sample set of the following numbers x1=0.09, x2=0.94,
x3=0.48, x4=0.74, x5=0.04 I can calculate the variance easily.
But if each x is actually a subset of multiple values, what would be
the formula to calculate the variance? and it is possible to implement
such mathematical function in R?

For instance if I have the following: x1=(0.77, 0.22, 0.44), x2=(0.26,
0.89, 0.58), x3=(0.20, 0.25, 0.91), x4=(0.06, 0.13, 0.26) and
x5=(0.65, 0.16, 0.72) how can i calculate the variance for each x?

Thank you

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 variance on multiple numbers at once?

2015-11-11 Thread S Ellison



> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Luigi
> Marongiu
> if I have a sample set of the following numbers x1=0.09, x2=0.94, x3=0.48,
> x4=0.74, x5=0.04 I can calculate the variance easily.
Not without concatenating them into a vector, you can't. You need them in a 
vector, as in
var( c(x1, x2, x3, x4, x5) )

> But if each x is actually a subset of multiple values, what would be the 
> formula
> to calculate the variance? and it is possible to implement such mathematical
> function in R?
This is what R wants anyway, so the function you are looking for is var()

> For instance if I have the following: x1=(0.77, 0.22, 0.44), x2=(0.26, 0.89, 
> 0.58),
> x3=(0.20, 0.25, 0.91), x4=(0.06, 0.13, 0.26) and x5=(0.65, 0.16, 0.72) how 
> can i
> calculate the variance for each x?
var(x1)
var(x2)


or, if you want to be a bit more slick about it and do it in one line

lapply(list( x1, x2, x3, ...), var  ) 

(or sapply() if you want a vector result)






***
This email and any attachments are confidential. Any use...{{dropped:8}}

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the shared area of 2 plots?

2015-10-31 Thread Viechtbauer Wolfgang (STAT)

Something like this?

http://stats.stackexchange.com/questions/12209/percentage-of-overlapping-regions-of-two-normal-distributions

Best,
Wolfgang

--
Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com

From: R-help [r-help-boun...@r-project.org] On Behalf Of Jackson Rodrigues 
[jacksonmrodrig...@gmail.com]
Sent: Saturday, October 31, 2015 5:18 PM
To: r-help@r-project.org
Subject: [R] How to calculate the shared area of 2 plots?

Dear all,

I have 2 data sets (A and B) with normal distributions. When I plot them
together, they overlay/share one part of plot.
I want to extract from each A and B that shared area.
Actually I want to sum A+B and then to compare the result with the
individual shared area of A and B to know how much (in %) each one
represents. (Am I clear?)

The folowing codes do not represent exactly my data, however the final
result (plot) is very similar to what I get in the end with my data.
I can provide my original data if necessary.

A<-seq(-4,4,.01)
B<-seq(-4,4,.01)

densities<-dnorm(A, 0,1)
densities2<-dnorm(B, 0,1)

plot(A,densities, type="l", col="red")
lines(B+2,densities2, type="l")

I am very grateful for any support, help, advise  and etc ... :)


best wishes

Jackson
--

Jackson M. Rodrigues
Department of Palynology and Climate Dynamics
Albrecht-von-Haller-Institute for Plant Sciences
Georg-August-University Göttingen
Untere Karspuele 2
37073 Göttingen/Germany
Tel.:   0049 (0) 176 8186 4994

*"In order to succeed, we must first believe that we can."*
*Nikos Kazantzakis*

[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate the shared area of 2 plots?

2015-10-31 Thread Jackson Rodrigues

Dear all,

I have 2 data sets (A and B) with normal distributions. When I plot them
together, they overlay/share one part of plot.
I want to extract from each A and B that shared area.
Actually I want to sum A+B and then to compare the result with the
individual shared area of A and B to know how much (in %) each one
represents. (Am I clear?)

The folowing codes do not represent exactly my data, however the final
result (plot) is very similar to what I get in the end with my data.
I can provide my original data if necessary.

A<-seq(-4,4,.01)
B<-seq(-4,4,.01)

densities<-dnorm(A, 0,1)
densities2<-dnorm(B, 0,1)

plot(A,densities, type="l", col="red")
lines(B+2,densities2, type="l")

I am very grateful for any support, help, advise  and etc ... :)


best wishes

Jackson
-- 

Jackson M. Rodrigues
Department of Palynology and Climate Dynamics
Albrecht-von-Haller-Institute for Plant Sciences
Georg-August-University Göttingen
Untere Karspuele 2
37073 Göttingen/Germany
Tel.:   0049 (0) 176 8186 4994

*"In order to succeed, we must first believe that we can."*
*Nikos Kazantzakis*

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 standard error of estimate (S) for my non-linear regression model?

2015-09-26 Thread peter dalgaard

This is one area in which terminology in (computational) statistics has gone a 
bit crazy. The thing some call "standard error of estimate" is actually the 
residual standard deviation in the regression model, not to be confused with 
the standard errors that are associated with parameter estimates. In 
summary(nls(...)) (and summary(lm()) for that matter), you'll find it as 
"residual standard error",  and even that is a bit of a misnomer.

-pd 

> On 26 Sep 2015, at 07:08 , Michael Eisenring  wrote:
> 
> Hi all,
> 
> I am looking for something that indicates the goodness of fit for my non
> linear regression model (since R2 is not very reliable).
> 
> I read that the standard error of estimate (also known as standard error of
> the regression) is a good alternative.
> 
> 
> 
> The standard error of estimate is described on this page (including the
> formula) http://onlinestatbook.com/2/regression/accuracy.html
>  atbook.com%2F2%2Fregression%2Faccuracy.html> 
> 
> Unfortunately however, I have no clue how to programm it in R. Does anyone
> know and could help me?
> 
> Thank you very much.
> 
> 
> 
> I added an example of my model and a dput() of my data
> 
> #CODE
> 
> dta<-read.csv("Regression_exp2.csv",header=T, sep = ",")
> attach(dta)  # tells R to do the following analyses on this dataset
> head(dta)
> 
> 
> 
> # loading packages: analysis of mixed effect models
> library(nls2)#model
> 
> #Aim: fit equation to data: y~yo+a*(1-b^x) : Two parameter exp. single rise
> to the maximum
> # y =Gossypol (from my data set) x= Damage_cm (from my data set)
> #The other 3 parameters are unknown: yo=Intercept, a= assymptote ans b=slope
> 
> plot(Gossypol~Damage_cm, dta)
> # Looking at the plot, 0 is a plausible estimate for y0:
> # a+y0 is the asymptote, so estimate about 4000;
> # b is between 0 and 1, so estimate .5
> dta.nls <- nls(Gossypol~y0+a*(1-b^Damage_cm), dta,
>   start=list(y0=0, a=4000, b=.5))
> 
> xval <- seq(0, 10, 0.1)
> lines(xval, predict(dta.nls, data.frame(Damage_cm=xval)))
> profile(dta.nls, alpha= .05)
> 
> 
> summary(dta.nls)
> 
> 
> 
> 
> 
> 
> 
> #INPUT
> 
> structure(list(Gossypol = c(948.2418407, 1180.171957, 3589.187889, 
> 450.7205451, 349.0864019, 592.3403778, 723.885643, 2005.919344, 
> 720.9785449, 1247.806111, 1079.846532, 1500.863038, 4198.569251, 
> 3618.448997, 4140.242559, 1036.331811, 1013.807628, 2547.326207, 
> 2508.417927, 2874.651764, 1120.955, 1782.864308, 1517.045807, 
> 2287.228752, 4171.427741, 3130.376482, 1504.491931, 6132.876396, 
> 3350.203452, 5113.942098, 1989.576826, 3470.09352, 4576.787021, 
> 4854.985845, 1414.161257, 2608.716056, 910.8879471, 2228.522959, 
> 2952.931863, 5909.068158, 1247.806111, 6982.035521, 2867.610671, 
> 5629.979049, 6039.995102, 3747.076592, 3743.331903, 4274.324792, 
> 3378.151945, 3736.144027, 5654.858696, 5972.926124, 3723.629772, 
> 3322.115942, 3575.043632, 2818.419785), Treatment = structure(c(5L, 
> 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 
> 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 3L, 
> 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 4L, 2L, 4L, 4L, 2L, 4L, 2L, 2L, 
> 4L, 4L, 4L, 4L, 4L, 4L, 2L), .Label = c("1c_2d", "1c_7d", "3c_2d", 
> "9c_2d", "C"), class = "factor"), Damage_cm = c(0, 0, 0, 0, 0, 
> 0, 0, 0, 0, 0, 0, 0.142, 0.4035, 0.4435, 0.491, 0.4955, 0.578, 
> 0.5895, 0.6925, 0.6965, 0.756, 0.8295, 1.0475, 1.313, 1.516, 
> 1.573, 1.62, 1.8115, 1.8185, 1.8595, 1.989, 2.129, 2.171, 2.3035, 
> 2.411, 2.559, 2.966, 2.974, 3.211, 3.2665, 3.474, 3.51, 3.547, 
> 4.023, 4.409, 4.516, 4.7245, 4.809, 4.9835, 5.568, 5.681, 5.683, 
> 7.272, 8.043, 9.437, 9.7455), Damage_groups = c(0.278, 1.616, 
> 2.501, 3.401, 4.577, 5.644, 7.272, 8.043, 9.591, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), Gossypol_Averaged =
> c(1783.211, 
> 3244.129, 2866.307, 3991.809, 4468.809, 5121.309, 3723.629772, 
> 3322.115942, 3196.731, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA), Groups = c(42006L, 42038L, 42067L, 42099L, 
> 42130L, 42162L, 42193L, 42225L, 42257L, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
> NA, NA, NA, NA, NA, NA, NA, NA, NA)), .Names = c("Gossypol", 
> "Treatment", "Damage_cm", "Damage_groups", "Gossypol_Averaged", 
> "Groups"), class = "data.frame", row.names = c(NA, -56L))
> 
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>

Re: [R] How to calculate standard error of estimate (S) for my non-linear regression model?

2015-09-26 Thread Bert Gunter

Michael:

You appear to be laboring under the illusion that a single numeric
summary (**any summary**)is a useful measure of model adequacy. It is
not; for details about why not, consult any applied statistics text
(e.g. on regression) and/or post on a statistics site, like
stats.stackexchange.com. Better yet, consult a local statistician.

Incidentally, this is even more the case for NON-linear models. Again,
consult appropriate statistical resources. Even googling on "R^2
inadequate for nonlinear models" brought up some interesting
resources, among them:

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2892436/

Cheers,
Bert

Bert Gunter

"Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom."
   -- Clifford Stoll


On Sat, Sep 26, 2015 at 8:56 AM, peter dalgaard <pda...@gmail.com> wrote:
>
>> On 26 Sep 2015, at 16:46 , Michael Eisenring <michael.eisenr...@gmx.ch> 
>> wrote:
>>
>> Dear Peter,
>> Thank you for your answer.
>> If I look at my summary I see there a Residual standard error: 1394 on 53
>> degrees of freedom.
>> This number is very high (the fit of the curve is pretty bad I know but
>> still...). Are you sure the residual standard error given in the summary is
>> the same as the one described on this page:
>> http://onlinestatbook.com/2/regression/accuracy.html
>
> Sure I'm sure (& I did check!)... But notice that unlike R^2, Residual SE is 
> not dimensionless. Switch from millimeters to meters in your response measure 
> and the Residual SE instantly turns into 1.394. It's basically saying that 
> your model claims to be able to predict Y within +/-2800 units. How good or 
> bad that is, is your judgment.
>
>> I am basically just looking for a value that describes the goodness of fit
>> for my non-linear regression model.
>>
>>
>> This is probably a pretty obvious question, but I am not a statistician and
>> as you said the terminology is sometimes pretty confusing.
>> Thanks mike
>>
>> -Ursprüngliche Nachricht-
>> Von: peter dalgaard [mailto:pda...@gmail.com]
>> Gesendet: Samstag, 26. September 2015 01:43
>> An: Michael Eisenring <michael.eisenr...@gmx.ch>
>> Cc: r-help@r-project.org
>> Betreff: Re: [R] How to calculate standard error of estimate (S) for my
>> non-linear regression model?
>>
>> This is one area in which terminology in (computational) statistics has gone
>> a bit crazy. The thing some call "standard error of estimate" is actually
>> the residual standard deviation in the regression model, not to be confused
>> with the standard errors that are associated with parameter estimates. In
>> summary(nls(...)) (and summary(lm()) for that matter), you'll find it as
>> "residual standard error",  and even that is a bit of a misnomer.
>>
>> -pd
>>
>>> On 26 Sep 2015, at 07:08 , Michael Eisenring <michael.eisenr...@gmx.ch>
>> wrote:
>>>
>>> Hi all,
>>>
>>> I am looking for something that indicates the goodness of fit for my
>>> non linear regression model (since R2 is not very reliable).
>>>
>>> I read that the standard error of estimate (also known as standard
>>> error of the regression) is a good alternative.
>>>
>>>
>>>
>>> The standard error of estimate is described on this page (including
>>> the
>>> formula) http://onlinestatbook.com/2/regression/accuracy.html
>>> <https://3c.gmx.net/mail/client/dereferrer?redirectUrl=http%3A%2F%2Fon
>>> linest atbook.com%2F2%2Fregression%2Faccuracy.html>
>>>
>>> Unfortunately however, I have no clue how to programm it in R. Does
>>> anyone know and could help me?
>>>
>>> Thank you very much.
>>>
>>>
>>>
>>> I added an example of my model and a dput() of my data
>>>
>>> #CODE
>>>
>>> dta<-read.csv("Regression_exp2.csv",header=T, sep = ",")
>>> attach(dta)  # tells R to do the following analyses on this dataset
>>> head(dta)
>>>
>>>
>>>
>>> # loading packages: analysis of mixed effect models
>>> library(nls2)#model
>>>
>>> #Aim: fit equation to data: y~yo+a*(1-b^x) : Two parameter exp. single
>>> rise to the maximum # y =Gossypol (from my data set) x= Damage_cm
>>> (from my data set) #The other 3 parameters are unknown: yo=Intercept,
>>> a= assymptote ans b=slope
>>>
>>> plot(Gossypol~Damage_cm, dta)
>>> # Looking at the plot, 0 is a plausible estim

Re: [R] How to calculate standard error of estimate (S) for my non-linear regression model?

2015-09-26 Thread Michael Eisenring

Dear Peter,
Thank you for your answer.
If I look at my summary I see there a Residual standard error: 1394 on 53
degrees of freedom.
This number is very high (the fit of the curve is pretty bad I know but
still...). Are you sure the residual standard error given in the summary is
the same as the one described on this page:
http://onlinestatbook.com/2/regression/accuracy.html
I am basically just looking for a value that describes the goodness of fit
for my non-linear regression model.


This is probably a pretty obvious question, but I am not a statistician and
as you said the terminology is sometimes pretty confusing.
Thanks mike

-Ursprüngliche Nachricht-
Von: peter dalgaard [mailto:pda...@gmail.com] 
Gesendet: Samstag, 26. September 2015 01:43
An: Michael Eisenring <michael.eisenr...@gmx.ch>
Cc: r-help@r-project.org
Betreff: Re: [R] How to calculate standard error of estimate (S) for my
non-linear regression model?

This is one area in which terminology in (computational) statistics has gone
a bit crazy. The thing some call "standard error of estimate" is actually
the residual standard deviation in the regression model, not to be confused
with the standard errors that are associated with parameter estimates. In
summary(nls(...)) (and summary(lm()) for that matter), you'll find it as
"residual standard error",  and even that is a bit of a misnomer.

-pd 

> On 26 Sep 2015, at 07:08 , Michael Eisenring <michael.eisenr...@gmx.ch>
wrote:
> 
> Hi all,
> 
> I am looking for something that indicates the goodness of fit for my 
> non linear regression model (since R2 is not very reliable).
> 
> I read that the standard error of estimate (also known as standard 
> error of the regression) is a good alternative.
> 
> 
> 
> The standard error of estimate is described on this page (including 
> the
> formula) http://onlinestatbook.com/2/regression/accuracy.html
> <https://3c.gmx.net/mail/client/dereferrer?redirectUrl=http%3A%2F%2Fon
> linest atbook.com%2F2%2Fregression%2Faccuracy.html>
> 
> Unfortunately however, I have no clue how to programm it in R. Does 
> anyone know and could help me?
> 
> Thank you very much.
> 
> 
> 
> I added an example of my model and a dput() of my data
> 
> #CODE
> 
> dta<-read.csv("Regression_exp2.csv",header=T, sep = ",")
> attach(dta)  # tells R to do the following analyses on this dataset
> head(dta)
> 
> 
> 
> # loading packages: analysis of mixed effect models 
> library(nls2)#model
> 
> #Aim: fit equation to data: y~yo+a*(1-b^x) : Two parameter exp. single 
> rise to the maximum # y =Gossypol (from my data set) x= Damage_cm 
> (from my data set) #The other 3 parameters are unknown: yo=Intercept, 
> a= assymptote ans b=slope
> 
> plot(Gossypol~Damage_cm, dta)
> # Looking at the plot, 0 is a plausible estimate for y0:
> # a+y0 is the asymptote, so estimate about 4000; # b is between 0 and 
> 1, so estimate .5 dta.nls <- nls(Gossypol~y0+a*(1-b^Damage_cm), dta,
>   start=list(y0=0, a=4000, b=.5))
> 
> xval <- seq(0, 10, 0.1)
> lines(xval, predict(dta.nls, data.frame(Damage_cm=xval))) 
> profile(dta.nls, alpha= .05)
> 
> 
> summary(dta.nls)
> 
> 
> 
> 
> 
> 
> 
> #INPUT
> 
> structure(list(Gossypol = c(948.2418407, 1180.171957, 3589.187889, 
> 450.7205451, 349.0864019, 592.3403778, 723.885643, 2005.919344, 
> 720.9785449, 1247.806111, 1079.846532, 1500.863038, 4198.569251, 
> 3618.448997, 4140.242559, 1036.331811, 1013.807628, 2547.326207, 
> 2508.417927, 2874.651764, 1120.955, 1782.864308, 1517.045807, 
> 2287.228752, 4171.427741, 3130.376482, 1504.491931, 6132.876396, 
> 3350.203452, 5113.942098, 1989.576826, 3470.09352, 4576.787021, 
> 4854.985845, 1414.161257, 2608.716056, 910.8879471, 2228.522959, 
> 2952.931863, 5909.068158, 1247.806111, 6982.035521, 2867.610671, 
> 5629.979049, 6039.995102, 3747.076592, 3743.331903, 4274.324792, 
> 3378.151945, 3736.144027, 5654.858696, 5972.926124, 3723.629772, 
> 3322.115942, 3575.043632, 2818.419785), Treatment = structure(c(5L, 
> 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
> 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
> 3L, 3L, 3L, 2L, 2L, 2L, 4L, 2L, 4L, 4L, 2L, 4L, 2L, 2L, 4L, 4L, 4L, 
> 4L, 4L, 4L, 2L), .Label = c("1c_2d", "1c_7d", "3c_2d", "9c_2d", "C"), 
> class = "factor"), Damage_cm = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
> 0.142, 0.4035, 0.4435, 0.491, 0.4955, 0.578, 0.5895, 0.6925, 0.6965, 
> 0.756, 0.8295, 1.0475, 1.313, 1.516, 1.573, 1.62, 1.8115, 1.8185, 
> 1.8595, 1.989, 2.129, 2.171, 2.3035, 2.411, 2.559, 2.966, 2.974, 
> 3.211, 3.2665, 3.474, 3.51, 3.547, 4.023, 4.409, 4.516, 4

Re: [R] How to calculate standard error of estimate (S) for my non-linear regression model?

2015-09-26 Thread peter dalgaard


> On 26 Sep 2015, at 16:46 , Michael Eisenring <michael.eisenr...@gmx.ch> wrote:
> 
> Dear Peter,
> Thank you for your answer.
> If I look at my summary I see there a Residual standard error: 1394 on 53
> degrees of freedom.
> This number is very high (the fit of the curve is pretty bad I know but
> still...). Are you sure the residual standard error given in the summary is
> the same as the one described on this page:
> http://onlinestatbook.com/2/regression/accuracy.html

Sure I'm sure (& I did check!)... But notice that unlike R^2, Residual SE is 
not dimensionless. Switch from millimeters to meters in your response measure 
and the Residual SE instantly turns into 1.394. It's basically saying that your 
model claims to be able to predict Y within +/-2800 units. How good or bad that 
is, is your judgment.

> I am basically just looking for a value that describes the goodness of fit
> for my non-linear regression model.
> 
> 
> This is probably a pretty obvious question, but I am not a statistician and
> as you said the terminology is sometimes pretty confusing.
> Thanks mike
> 
> -Ursprüngliche Nachricht-
> Von: peter dalgaard [mailto:pda...@gmail.com] 
> Gesendet: Samstag, 26. September 2015 01:43
> An: Michael Eisenring <michael.eisenr...@gmx.ch>
> Cc: r-help@r-project.org
> Betreff: Re: [R] How to calculate standard error of estimate (S) for my
> non-linear regression model?
> 
> This is one area in which terminology in (computational) statistics has gone
> a bit crazy. The thing some call "standard error of estimate" is actually
> the residual standard deviation in the regression model, not to be confused
> with the standard errors that are associated with parameter estimates. In
> summary(nls(...)) (and summary(lm()) for that matter), you'll find it as
> "residual standard error",  and even that is a bit of a misnomer.
> 
> -pd 
> 
>> On 26 Sep 2015, at 07:08 , Michael Eisenring <michael.eisenr...@gmx.ch>
> wrote:
>> 
>> Hi all,
>> 
>> I am looking for something that indicates the goodness of fit for my 
>> non linear regression model (since R2 is not very reliable).
>> 
>> I read that the standard error of estimate (also known as standard 
>> error of the regression) is a good alternative.
>> 
>> 
>> 
>> The standard error of estimate is described on this page (including 
>> the
>> formula) http://onlinestatbook.com/2/regression/accuracy.html
>> <https://3c.gmx.net/mail/client/dereferrer?redirectUrl=http%3A%2F%2Fon
>> linest atbook.com%2F2%2Fregression%2Faccuracy.html>
>> 
>> Unfortunately however, I have no clue how to programm it in R. Does 
>> anyone know and could help me?
>> 
>> Thank you very much.
>> 
>> 
>> 
>> I added an example of my model and a dput() of my data
>> 
>> #CODE
>> 
>> dta<-read.csv("Regression_exp2.csv",header=T, sep = ",")
>> attach(dta)  # tells R to do the following analyses on this dataset
>> head(dta)
>> 
>> 
>> 
>> # loading packages: analysis of mixed effect models 
>> library(nls2)#model
>> 
>> #Aim: fit equation to data: y~yo+a*(1-b^x) : Two parameter exp. single 
>> rise to the maximum # y =Gossypol (from my data set) x= Damage_cm 
>> (from my data set) #The other 3 parameters are unknown: yo=Intercept, 
>> a= assymptote ans b=slope
>> 
>> plot(Gossypol~Damage_cm, dta)
>> # Looking at the plot, 0 is a plausible estimate for y0:
>> # a+y0 is the asymptote, so estimate about 4000; # b is between 0 and 
>> 1, so estimate .5 dta.nls <- nls(Gossypol~y0+a*(1-b^Damage_cm), dta,
>>  start=list(y0=0, a=4000, b=.5))
>> 
>> xval <- seq(0, 10, 0.1)
>> lines(xval, predict(dta.nls, data.frame(Damage_cm=xval))) 
>> profile(dta.nls, alpha= .05)
>> 
>> 
>> summary(dta.nls)
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> #INPUT
>> 
>> structure(list(Gossypol = c(948.2418407, 1180.171957, 3589.187889, 
>> 450.7205451, 349.0864019, 592.3403778, 723.885643, 2005.919344, 
>> 720.9785449, 1247.806111, 1079.846532, 1500.863038, 4198.569251, 
>> 3618.448997, 4140.242559, 1036.331811, 1013.807628, 2547.326207, 
>> 2508.417927, 2874.651764, 1120.955, 1782.864308, 1517.045807, 
>> 2287.228752, 4171.427741, 3130.376482, 1504.491931, 6132.876396, 
>> 3350.203452, 5113.942098, 1989.576826, 3470.09352, 4576.787021, 
>> 4854.985845, 1414.161257, 2608.716056, 910.8879471, 2228.522959, 
>> 2952.931863, 5909.068158, 1247.806111, 6982.035521, 2867.610671, 
>> 5629.979

Re: [R] How to calculate standard error of estimate (S) for my non-linear regression model?

2015-09-26 Thread Jeff Newmiller

You may not be a statistician, but you should at least learn about the 
calculations you are making. You cannot expect to convince others that 
your calculations are right just because "Peter on the internet said they 
were right".


To give you a gentle push in this direction, I have reproduced the 
calculations on that reference web page using R, so you can get a head 
start on understanding how to perform them on your data. (Hint: about all 
you should need to do is give your dta.nls to predict and use your 
column names instead of X and Y.) Once you have convinced yourself that
the pre-defined functions are doing what you expect, then you can omit the 
do-it-yourself option with confidence in the results.


Please note that no use of attach is included here... that usually ends in 
unhappiness at some point, so prefer to use the data argument instead.


###
dta <- data.frame( X=c(1,2,3,4,5), Y=c(1,2,1.3,3.75,2.25) )
nrow( dta )
# linear regression
y.lm <- lm( Y~X, data=dta )
# compute predictions from original data
## Be aware that you can also "predict" using a nonlinear regression fit
## Also be aware that you can compute estimates using different data if you
## specify the "newdata" argument... see the help for predict.lm
## ?predict.lm
dta$Yprime <- predict( y.lm )
# ?with
dta$YlessYprime <- with( dta, Y - Yprime )
dta$YlessYprime2 <- with( dta, YlessYprime^2 )

# confirm sums
# ?sum
sum( dta$X )
sum( dta$Y )
sum( dta$Yprime )
sum( dta$YlessYprime )
sum( dta$YlessYprime2 )

# standard error of the estimate for population data
sigma.est <- sqrt( sum( dta$YlessYprime2 ) / nrow( dta ) )
sigma.est
# sd function assumes sample standard deviation, can correct the result if 
# you want

# ?sd
# ?all.equal
all.equal( sigma.est, sd( dta$YlessYprime ) * sqrt( ( nrow( dta ) - 1 ) / 
nrow( dta ) ) )


# alternate formulation
SSY <- sum( ( dta$Y - mean( dta$Y ) )^2 )
rho <- with( dta, cor( Y, X ) )
all.equal( sigma.est, sqrt( (1-rho^2)*SSY/nrow(dta) ) )

# when working with a sample...
s.est <- sqrt( sum( dta$YlessYprime2 ) / ( nrow( dta ) - 2 ) )
s.est



dta <- data.frame( X=c(1,2,3,4,5), Y=c(1,2,1.3,3.75,2.25) )
nrow( dta )

[1] 5

# linear regression
y.lm <- lm( Y~X, data=dta )
# compute predictions from original data
## Be aware that you can also "predict" using a nonlinear regression fit
## Also be aware that you can compute estimates using different data if you
## specify the "newdata" argument... see the help for predict.lm
## ?predict.lm
dta$Yprime <- predict( y.lm )
# ?with
dta$YlessYprime <- with( dta, Y - Yprime )
dta$YlessYprime2 <- with( dta, YlessYprime^2 )

# confirm sums
# ?sum
sum( dta$X )

[1] 15

sum( dta$Y )

[1] 10.3

sum( dta$Yprime )

[1] 10.3

sum( dta$YlessYprime )

[1] 2.220446e-16

sum( dta$YlessYprime2 )

[1] 2.79075


# standard error of the estimate for population data
sigma.est <- sqrt( sum( dta$YlessYprime2 ) / nrow( dta ) )
sigma.est

[1] 0.7470944
# sd function assumes sample standard deviation, can correct the result 
# if you want

# ?sd
# ?all.equal
all.equal( sigma.est, sd( dta$YlessYprime ) * sqrt( ( nrow( dta ) - 1 ) 

/ nrow( dta ) ) )
[1] TRUE


# alternate formulation
SSY <- sum( ( dta$Y - mean( dta$Y ) )^2 )
rho <- with( dta, cor( Y, X ) )
all.equal( sigma.est, sqrt( (1-rho^2)*SSY/nrow(dta) ) )

[1] TRUE


# when working with a sample...
s.est <- sqrt( sum( dta$YlessYprime2 ) / ( nrow( dta ) - 2 ) )
s.est

[1] 0.9644947





On Sat, 26 Sep 2015, Michael Eisenring wrote:


Dear Peter,
Thank you for your answer.
If I look at my summary I see there a Residual standard error: 1394 on 53
degrees of freedom.
This number is very high (the fit of the curve is pretty bad I know but
still...). Are you sure the residual standard error given in the summary is
the same as the one described on this page:
http://onlinestatbook.com/2/regression/accuracy.html
I am basically just looking for a value that describes the goodness of fit
for my non-linear regression model.


This is probably a pretty obvious question, but I am not a statistician and
as you said the terminology is sometimes pretty confusing.
Thanks mike

-Urspr?ngliche Nachricht-
Von: peter dalgaard [mailto:pda...@gmail.com]
Gesendet: Samstag, 26. September 2015 01:43
An: Michael Eisenring <michael.eisenr...@gmx.ch>
Cc: r-help@r-project.org
Betreff: Re: [R] How to calculate standard error of estimate (S) for my
non-linear regression model?

This is one area in which terminology in (computational) statistics has gone
a bit crazy. The thing some call "standard error of estimate" is actually
the residual standard deviation in the regression model, not to be confused
with the standard errors that are associated with parameter estimates. In
summary(nls(...)) (and summary(lm()) for that matter), you'll find it as
"residual standard error",  and even that is a bit of a misnomer.

-pd


On

[R] How to calculate standard error of estimate (S) for my non-linear regression model?

2015-09-25 Thread Michael Eisenring

Hi all,

I am looking for something that indicates the goodness of fit for my non
linear regression model (since R2 is not very reliable).

I read that the standard error of estimate (also known as standard error of
the regression) is a good alternative.

 

The standard error of estimate is described on this page (including the
formula) http://onlinestatbook.com/2/regression/accuracy.html
 

Unfortunately however, I have no clue how to programm it in R. Does anyone
know and could help me?

Thank you very much.

 

I added an example of my model and a dput() of my data

#CODE

dta<-read.csv("Regression_exp2.csv",header=T, sep = ",")
attach(dta)  # tells R to do the following analyses on this dataset
head(dta)

 

# loading packages: analysis of mixed effect models
library(nls2)#model

#Aim: fit equation to data: y~yo+a*(1-b^x) : Two parameter exp. single rise
to the maximum
# y =Gossypol (from my data set) x= Damage_cm (from my data set)
#The other 3 parameters are unknown: yo=Intercept, a= assymptote ans b=slope

plot(Gossypol~Damage_cm, dta)
# Looking at the plot, 0 is a plausible estimate for y0:
# a+y0 is the asymptote, so estimate about 4000;
# b is between 0 and 1, so estimate .5
dta.nls <- nls(Gossypol~y0+a*(1-b^Damage_cm), dta,
   start=list(y0=0, a=4000, b=.5))

xval <- seq(0, 10, 0.1)
lines(xval, predict(dta.nls, data.frame(Damage_cm=xval)))
profile(dta.nls, alpha= .05)


summary(dta.nls)

 

 

 

#INPUT

structure(list(Gossypol = c(948.2418407, 1180.171957, 3589.187889, 
450.7205451, 349.0864019, 592.3403778, 723.885643, 2005.919344, 
720.9785449, 1247.806111, 1079.846532, 1500.863038, 4198.569251, 
3618.448997, 4140.242559, 1036.331811, 1013.807628, 2547.326207, 
2508.417927, 2874.651764, 1120.955, 1782.864308, 1517.045807, 
2287.228752, 4171.427741, 3130.376482, 1504.491931, 6132.876396, 
3350.203452, 5113.942098, 1989.576826, 3470.09352, 4576.787021, 
4854.985845, 1414.161257, 2608.716056, 910.8879471, 2228.522959, 
2952.931863, 5909.068158, 1247.806111, 6982.035521, 2867.610671, 
5629.979049, 6039.995102, 3747.076592, 3743.331903, 4274.324792, 
3378.151945, 3736.144027, 5654.858696, 5972.926124, 3723.629772, 
3322.115942, 3575.043632, 2818.419785), Treatment = structure(c(5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 4L, 2L, 4L, 4L, 2L, 4L, 2L, 2L, 
4L, 4L, 4L, 4L, 4L, 4L, 2L), .Label = c("1c_2d", "1c_7d", "3c_2d", 
"9c_2d", "C"), class = "factor"), Damage_cm = c(0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0.142, 0.4035, 0.4435, 0.491, 0.4955, 0.578, 
0.5895, 0.6925, 0.6965, 0.756, 0.8295, 1.0475, 1.313, 1.516, 
1.573, 1.62, 1.8115, 1.8185, 1.8595, 1.989, 2.129, 2.171, 2.3035, 
2.411, 2.559, 2.966, 2.974, 3.211, 3.2665, 3.474, 3.51, 3.547, 
4.023, 4.409, 4.516, 4.7245, 4.809, 4.9835, 5.568, 5.681, 5.683, 
7.272, 8.043, 9.437, 9.7455), Damage_groups = c(0.278, 1.616, 
2.501, 3.401, 4.577, 5.644, 7.272, 8.043, 9.591, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), Gossypol_Averaged =
c(1783.211, 
3244.129, 2866.307, 3991.809, 4468.809, 5121.309, 3723.629772, 
3322.115942, 3196.731, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA), Groups = c(42006L, 42038L, 42067L, 42099L, 
42130L, 42162L, 42193L, 42225L, 42257L, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA)), .Names = c("Gossypol", 
"Treatment", "Damage_cm", "Damage_groups", "Gossypol_Averaged", 
"Groups"), class = "data.frame", row.names = c(NA, -56L))

 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the average direct effect, average total effect and average indirect effect for spatial regression models with spatial lag of dependent variable

2015-07-24 Thread Linus Holtermann

Hello,

the command impacts in the spdep-packeage should help you.

Mit freundlichen Grüßen

Linus Holtermann
Hamburgisches WeltWirtschaftsInstitut gemeinnützige GmbH (HWWI)
Heimhuder Straße 71
20148 Hamburg
Tel +49-(0)40-340576-336
Fax+49-(0)40-340576-776
Internet: www.hwwi.org
Email: holterm...@hwwi.org
Amtsgericht Hamburg HRB 94303
Geschäftsführer: PD Dr. Christian Growitsch | Prof. Dr. Henning Vöpel
Prokura: Dipl. Kauffrau Alexis Malchin
Umsatzsteuer-ID: DE 241849425

-Ursprüngliche Nachricht-
Von: R-help [mailto:r-help-boun...@r-project.org] Im Auftrag von wenyueyang
Gesendet: Donnerstag, 23. Juli 2015 17:35
An: r-help@r-project.org
Betreff: [R] How to calculate the average direct effect, average total effect
and average indirect effect for spatial regression models with spatial lag of
dependent variable

Hi,

I am using four spatial regression models (SAR, SEM, SAC and SDM) to evaluate
the spillover effect of some factors. LeSage and Pace (2009) pointed out that
when the spatial lag of the dependent variable is included in the model,
parameter estimates lose their conventional interpretation as marginal effects,
because the spatial lag gives rise to a series of feedback loops and spillover
effects across regions. Therefore, I need to calculate the three different
marginal effects: average direct effect, average total effect and average
indirect (spillover) effect, and this is what I don't know how to perform in R.
Can you teach me about how to calculate the average direct effect, average
total effect and average indirect (spillover) effect for the spatial regression
models in R and tell me the related R code?
I would like to express my great appreciation to you!

Thank you very much and best regards.

Yours sincerely,
Wenyue Yang

--
View this message in context:
http://r.789695.n4.nabble.com/How-to-calculate-the-average-direct-effect-average-total-effect-and-average-indirect-effect-for-spate-tp4710253.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 the average direct effect, average total effect and average indirect effect for spatial regression models with spatial lag of dependent variable

2015-07-24 Thread wenyueyang

Dear Linus Holtermann,

thank a lot for telling me the key to calculate the effects. I will try it.
thank you for your help!

best regards,

Wenyue Yang





--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-calculate-the-average-direct-effect-average-total-effect-and-average-indirect-effect-for-spate-tp4710253p4710315.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate the average direct effect, average total effect and average indirect effect for spatial regression models with spatial lag of dependent variable

2015-07-23 Thread wenyueyang

Hi,

I am using four spatial regression models (SAR, SEM, SAC and SDM) to
evaluate the spillover effect of some factors. LeSage and Pace (2009)
pointed out that when the spatial lag of the dependent variable is included
in the model, parameter estimates lose their conventional interpretation as
marginal effects, because the spatial lag gives rise to a series of feedback
loops and spillover effects across regions. Therefore, I need to calculate
the three different marginal effects: average direct effect, average total
effect and average indirect (spillover) effect, and this is what I don't
know how to perform in R.
Can you teach me about how to calculate the average direct effect, average
total effect and average indirect (spillover) effect for the spatial
regression models in R and tell me the related R code? 
I would like to express my great appreciation to you!

Thank you very much and best regards.

Yours sincerely,
Wenyue Yang



--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-calculate-the-average-direct-effect-average-total-effect-and-average-indirect-effect-for-spate-tp4710253.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 vif of each term of model in R?

2015-04-27 Thread PIKAL Petr

, 2015 3:11 PM
 To: Methekar, Pushpa (GE Transportation, Non-GE); r-help@r-
project.org
 Subject: RE: How to calculate vif of each term of model in R?

 Hi

 I did not see any answer so I try.

 Your question lacks some info:

 Which vif - car or HH?

 answers and comments in line

  -Original Message-
  From: R-help [mailto:r-help-boun...@r-project.org] On Behalf
  Of Methekar, Pushpa (GE Transportation, Non-GE)
  Sent: Wednesday, April 08, 2015 10:24 AM
  To: r-help@r-project.org
  Subject: [R] How to calculate vif of each term of model in R?

  I am beginner in R doing modelling in R, I loaded excel sheet
  in
  R, i have chosen x elements and y elements then fitted model
  for
linear
 and
  second order regression. Now I have both models. I am bit
   confused
 how
  to calculate vif for each term in model like

  e.g model1-lm(y1~x1+x2+.x9) when I am using rms package
  then
 it's
  giving me like

  vif(model1)

 x1 x2 x3 x4 x5
  x6
  x7

   6.679692   1.520271   1.667125   3.618439   4.931810
  2.073879
  13.870630

  x8 x9

 220.969628 214.034135

  now i want to compare each term with std vif as vif=10 and
  which
 will
  satisfy this condition i want to delete that term and update
model1.
 i
  have done something like this

  fun = function(model1) {

   for(i in 1:length(model1)){

v=vif(model1)

   ss=any(v[i]=10)

 here you select only one item from vif, Why do you use any?

  if(ss==1){update(model1,.~.,-v[i])}

  else{print(no update)}

 Why do you change i here?

   i-i+1

  }

  return(model1)

}

 if you want to get rid of all terms bigger than some threshold
 in
once
 you can use

 sel - which(vif(model1)10)

 and select values for update possibly by

 update(model1,.~. - names(vif(model1))[sel])

 or if you want to get rid one by one you can use

 vmax - which.max(vif(model1))
 and check if max vif value is bigger than 10.

 vif(model1)[vmax]=10

 If it is just update with

 - names(vif(model1))[vmax])

 if it is not do not update.

 All of this untested.

 Cheers
 Petr

  fun(model1)

  but giving error as

  Error in if (ss == 1) { : missing value where TRUE/FALSE
  needed.

  please tell me how do i solve this problem.

Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny 
pouze jeho adresátům.
Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
svého systému.
Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
zpožděním přenosu e-mailu.

V případě, že je tento e-mail součástí obchodního jednání:
- vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a 
to z jakéhokoliv důvodu i bez uvedení důvodu.
- a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce 
s dodatkem či odchylkou.
- trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
dosažením shody na všech jejích náležitostech.
- odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně 
pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně 
osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi 
či osobě jím zastoupené známá.

This e-mail and any documents attached to it may be confidential and are 
intended only for its intended recipients.
If you received this e-mail by mistake, please immediately inform its sender. 
Delete the contents of this e-mail with all attachments and its copies from 
your system.
If you are not the intended recipient of this e-mail, you are not authorized to 
use, disseminate, copy or disclose this e-mail in any manner.
The sender of this e-mail shall not be liable for any possible damage caused by 
modifications of the e-mail or by delay with transfer of the email.

In case that this e-mail forms part of business dealings:
- the sender reserves the right to end negotiations about entering into a 
contract in any time, for any

Re: [R] How to calculate vif of each term of model in R?

2015-04-23 Thread PIKAL Petr

: Methekar, Pushpa (GE Transportation, Non-GE)
[mailto:pushpa.methe...@ge.com]
Sent: Friday, April 17, 2015 11:43 AM
To: PIKAL Petr
Subject: RE: How to calculate vif of each term of model in R?

Car package

-Original Message-
From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
Sent: Friday, April 17, 2015 3:11 PM
To: Methekar, Pushpa (GE Transportation, Non-GE); r-help@r-
   project.org
Subject: RE: How to calculate vif of each term of model in R?

Hi

I did not see any answer so I try.

Your question lacks some info:

Which vif - car or HH?

answers and comments in line

 -Original Message-
 From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
 Methekar, Pushpa (GE Transportation, Non-GE)
 Sent: Wednesday, April 08, 2015 10:24 AM
 To: r-help@r-project.org
 Subject: [R] How to calculate vif of each term of model in R?

 I am beginner in R doing modelling in R, I loaded excel sheet
 in
 R, i have chosen x elements and y elements then fitted model
 for
   linear
and
 second order regression. Now I have both models. I am bit
  confused
how
 to calculate vif for each term in model like

 e.g model1-lm(y1~x1+x2+.x9) when I am using rms package
 then
it's
 giving me like

 vif(model1)

x1 x2 x3 x4 x5
 x6
 x7

  6.679692   1.520271   1.667125   3.618439   4.931810
 2.073879
 13.870630

 x8 x9

220.969628 214.034135

 now i want to compare each term with std vif as vif=10 and
 which
will
 satisfy this condition i want to delete that term and update
   model1.
i
 have done something like this

 fun = function(model1) {

  for(i in 1:length(model1)){

   v=vif(model1)

  ss=any(v[i]=10)

here you select only one item from vif, Why do you use any?

 if(ss==1){update(model1,.~.,-v[i])}

 else{print(no update)}

Why do you change i here?

  i-i+1

 }

 return(model1)

   }

if you want to get rid of all terms bigger than some threshold in
   once
you can use

sel - which(vif(model1)10)

and select values for update possibly by

update(model1,.~. - names(vif(model1))[sel])

or if you want to get rid one by one you can use

vmax - which.max(vif(model1))
and check if max vif value is bigger than 10.

vif(model1)[vmax]=10

If it is just update with

- names(vif(model1))[vmax])

if it is not do not update.

All of this untested.

Cheers
Petr

 fun(model1)

 but giving error as

 Error in if (ss == 1) { : missing value where TRUE/FALSE
 needed.

 please tell me how do i solve this problem.

   [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list -- To UNSUBSCRIBE and more,
 see https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-
  project.org/posting-
 guide.html and provide commented, minimal, self-contained,
 reproducible code.

 Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou
 určeny pouze jeho adresátům.
 Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě
 neprodleně jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho
 kopie vymažte ze svého systému.
 Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento
 email jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
 Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou
 modifikacemi či zpožděním přenosu e-mailu.

 V případě, že je tento e-mail součástí obchodního jednání:
 - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření
 smlouvy, a to z jakéhokoliv důvodu i bez uvedení důvodu.
 - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně
 přijmout; Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky
 ze strany příjemce s dodatkem či odchylkou.
 - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve
 výslovným dosažením shody na všech jejích náležitostech.
 - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za
 společnost žádné smlouvy s výjimkou případů, kdy k tomu byl písemně
 zmocněn nebo písemně pověřen a takové pověření nebo plná moc byly
 adresátovi tohoto emailu případně osobě, kterou adresát zastupuje,
 předloženy nebo jejich existence je adresátovi či osobě jím zastoupené
 známá.

 This e-mail and any documents attached to it may be confidential and
 are intended only for its intended recipients.
 If you

Re: [R] How to calculate vif of each term of model in R?

2015-04-17 Thread PIKAL Petr

Hi

I did not see any answer so I try.

Your question lacks some info:

Which vif - car or HH?

answers and comments in line

 -Original Message-
 From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
 Methekar, Pushpa (GE Transportation, Non-GE)
 Sent: Wednesday, April 08, 2015 10:24 AM
 To: r-help@r-project.org
 Subject: [R] How to calculate vif of each term of model in R?


 I am beginner in R doing modelling in R, I loaded excel sheet in R, i
 have chosen x elements and y elements then fitted model for linear and
 second order regression. Now I have both models. I am bit confused how
 to calculate vif for each term in model like

 e.g model1-lm(y1~x1+x2+.x9) when I am using rms package then it's
 giving me like

 vif(model1)

x1 x2 x3 x4 x5 x6
 x7

  6.679692   1.520271   1.667125   3.618439   4.931810   2.073879
 13.870630

 x8 x9

220.969628 214.034135

 now i want to compare each term with std vif as vif=10 and which will
 satisfy this condition i want to delete that term and update model1. i
 have done something like this

 fun = function(model1) {

  for(i in 1:length(model1)){

   v=vif(model1)

  ss=any(v[i]=10)

here you select only one item from vif, Why do you use any?


 if(ss==1){update(model1,.~.,-v[i])}

 else{print(no update)}


Why do you change i here?

  i-i+1

 }



 return(model1)

   }


if you want to get rid of all terms bigger than some threshold in once you can 
use

sel - which(vif(model1)10)

and select values for update possibly by

update(model1,.~. - names(vif(model1))[sel])

or if you want to get rid one by one you can use

vmax - which.max(vif(model1))
and check if max vif value is bigger than 10.

vif(model1)[vmax]=10

If it is just update with

- names(vif(model1))[vmax])

if it is not do not update.

All of this untested.

Cheers
Petr

 fun(model1)

 but giving error as

 Error in if (ss == 1) { : missing value where TRUE/FALSE needed.

 please tell me how do i solve this problem.



   [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
 and provide commented, minimal, self-contained, reproducible code.


Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny 
pouze jeho adresátům.
Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
svého systému.
Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
zpožděním přenosu e-mailu.

V případě, že je tento e-mail součástí obchodního jednání:
- vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a 
to z jakéhokoliv důvodu i bez uvedení důvodu.
- a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce 
s dodatkem či odchylkou.
- trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
dosažením shody na všech jejích náležitostech.
- odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně 
pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně 
osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi 
či osobě jím zastoupené známá.

This e-mail and any documents attached to it may be confidential and are 
intended only for its intended recipients.
If you received this e-mail by mistake, please immediately inform its sender. 
Delete the contents of this e-mail with all attachments and its copies from 
your system.
If you are not the intended recipient of this e-mail, you are not authorized to 
use, disseminate, copy or disclose this e-mail in any manner.
The sender of this e-mail shall not be liable for any possible damage caused by 
modifications of the e-mail or by delay with transfer of the email.

In case that this e-mail forms part of business dealings:
- the sender reserves the right to end negotiations about entering into a 
contract in any time, for any reason, and without stating any reasoning.
- if the e-mail contains an offer, the recipient is entitled to immediately 
accept such offer; The sender of this e-mail (offer) excludes any acceptance of 
the offer on the part of the recipient containing any amendment or variation.
- the sender insists on that the respective contract is concluded only upon an 
express mutual agreement

Re: [R] How to calculate vif of each term of model in R?

2015-04-17 Thread PIKAL Petr

Hi

 -Original Message-
 From: Methekar, Pushpa (GE Transportation, Non-GE)
 [mailto:pushpa.methe...@ge.com]
 Sent: Friday, April 17, 2015 1:12 PM
 To: PIKAL Petr
 Subject: RE: How to calculate vif of each term of model in R?

 Hi Petr,
 You got my problem ,the solution which u specified is little good but
 with


 update(model1,.~. - names(vif(model1))[vmax])



 I won't be able to update my model .
 Instead I have to write like


  update(model1,.~. - x8)


maybe something like

fun = function(model1) {

vmax - which.max(vif(model1))

while( vif(model1)[vmax]=10) {

model1 - update(model1,.~. - names(vif(model1))[vmax]))
vmax - which.max(vif(model1))

}

return(model1)

}

I am not sure about cycle (i did not use while in R for a while).

Without some data I cannot check syntax so it is up to you.

Cheers
Petr



 i.e my x which having highest vif value.
 Each time for removing highest x .i  have to explicitly write ..-
 x8) So is there any way to avoid this?




 Thanks,
 Pushpa

 -Original Message-
 From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
 Sent: Friday, April 17, 2015 4:16 PM
 To: Methekar, Pushpa (GE Transportation, Non-GE)
 Subject: RE: How to calculate vif of each term of model in R?

 Comments to your first mail which are among lines and directly related
 to your code.

 I specifically mentioned it

  answers and comments in line

 If you have my first respond you can find them easier then now as they
 are buried within your mail.

 Cheers
 Petr


  -Original Message-
  From: Methekar, Pushpa (GE Transportation, Non-GE)
  [mailto:pushpa.methe...@ge.com]
  Sent: Friday, April 17, 2015 12:10 PM
  To: PIKAL Petr
  Subject: RE: How to calculate vif of each term of model in R?
 
  What comments are you talking about?
 
  -Original Message-
  From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
  Sent: Friday, April 17, 2015 3:38 PM
  To: Methekar, Pushpa (GE Transportation, Non-GE)
  Cc: r-help@r-project.org
  Subject: RE: How to calculate vif of each term of model in R?
 
  Did you follow my advice/comments?
 
  Cheers
  Petr
 
 
   -Original Message-
   From: Methekar, Pushpa (GE Transportation, Non-GE)
   [mailto:pushpa.methe...@ge.com]
   Sent: Friday, April 17, 2015 11:43 AM
   To: PIKAL Petr
   Subject: RE: How to calculate vif of each term of model in R?
  
   Car package
  
   -Original Message-
   From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
   Sent: Friday, April 17, 2015 3:11 PM
   To: Methekar, Pushpa (GE Transportation, Non-GE); r-help@r-
  project.org
   Subject: RE: How to calculate vif of each term of model in R?
  
   Hi
  
   I did not see any answer so I try.
  
   Your question lacks some info:
  
   Which vif - car or HH?
  
   answers and comments in line
  
-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
Methekar, Pushpa (GE Transportation, Non-GE)
Sent: Wednesday, April 08, 2015 10:24 AM
To: r-help@r-project.org
Subject: [R] How to calculate vif of each term of model in R?
   
   
I am beginner in R doing modelling in R, I loaded excel sheet in
R, i have chosen x elements and y elements then fitted model for
  linear
   and
second order regression. Now I have both models. I am bit
 confused
   how
to calculate vif for each term in model like
   
e.g model1-lm(y1~x1+x2+.x9) when I am using rms package then
   it's
giving me like
   
vif(model1)
   
   x1 x2 x3 x4 x5 x6
x7
   
 6.679692   1.520271   1.667125   3.618439   4.931810   2.073879
13.870630
   
x8 x9
   
   220.969628 214.034135
   
now i want to compare each term with std vif as vif=10 and which
   will
satisfy this condition i want to delete that term and update
  model1.
   i
have done something like this
   
fun = function(model1) {
   
 for(i in 1:length(model1)){
   
  v=vif(model1)
   
 ss=any(v[i]=10)
  
   here you select only one item from vif, Why do you use any?
  
   
if(ss==1){update(model1,.~.,-v[i])}
   
else{print(no update)}
   
  
   Why do you change i here?
  
 i-i+1
   
}
   
   
   
return(model1)
   
  }
   
  
   if you want to get rid of all terms bigger than some threshold in
  once
   you can use
  
   sel - which(vif(model1)10)
  
   and select values for update possibly by
  
   update(model1,.~. - names(vif(model1))[sel])
  
   or if you want to get rid one by one you can use
  
   vmax - which.max(vif(model1))
   and check if max vif value is bigger than 10.
  
   vif(model1)[vmax]=10
  
   If it is just update with
  
   - names(vif(model1))[vmax])
  
   if it is not do not update.
  
   All of this untested.
  
   Cheers
   Petr
  
fun(model1)
   
but giving error as
   
Error in if (ss == 1

Re: [R] How to calculate vif of each term of model in R?

2015-04-17 Thread PIKAL Petr

Did you follow my advice/comments?

Cheers
Petr


 -Original Message-
 From: Methekar, Pushpa (GE Transportation, Non-GE)
 [mailto:pushpa.methe...@ge.com]
 Sent: Friday, April 17, 2015 11:43 AM
 To: PIKAL Petr
 Subject: RE: How to calculate vif of each term of model in R?

 Car package

 -Original Message-
 From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
 Sent: Friday, April 17, 2015 3:11 PM
 To: Methekar, Pushpa (GE Transportation, Non-GE); r-help@r-project.org
 Subject: RE: How to calculate vif of each term of model in R?

 Hi

 I did not see any answer so I try.

 Your question lacks some info:

 Which vif - car or HH?

 answers and comments in line

  -Original Message-
  From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
  Methekar, Pushpa (GE Transportation, Non-GE)
  Sent: Wednesday, April 08, 2015 10:24 AM
  To: r-help@r-project.org
  Subject: [R] How to calculate vif of each term of model in R?
 
 
  I am beginner in R doing modelling in R, I loaded excel sheet in R, i
  have chosen x elements and y elements then fitted model for linear
 and
  second order regression. Now I have both models. I am bit confused
 how
  to calculate vif for each term in model like
 
  e.g model1-lm(y1~x1+x2+.x9) when I am using rms package then
 it's
  giving me like
 
  vif(model1)
 
 x1 x2 x3 x4 x5 x6
  x7
 
   6.679692   1.520271   1.667125   3.618439   4.931810   2.073879
  13.870630
 
  x8 x9
 
 220.969628 214.034135
 
  now i want to compare each term with std vif as vif=10 and which
 will
  satisfy this condition i want to delete that term and update model1.
 i
  have done something like this
 
  fun = function(model1) {
 
   for(i in 1:length(model1)){
 
v=vif(model1)
 
   ss=any(v[i]=10)

 here you select only one item from vif, Why do you use any?

 
  if(ss==1){update(model1,.~.,-v[i])}
 
  else{print(no update)}
 

 Why do you change i here?

   i-i+1
 
  }
 
 
 
  return(model1)
 
}
 

 if you want to get rid of all terms bigger than some threshold in once
 you can use

 sel - which(vif(model1)10)

 and select values for update possibly by

 update(model1,.~. - names(vif(model1))[sel])

 or if you want to get rid one by one you can use

 vmax - which.max(vif(model1))
 and check if max vif value is bigger than 10.

 vif(model1)[vmax]=10

 If it is just update with

 - names(vif(model1))[vmax])

 if it is not do not update.

 All of this untested.

 Cheers
 Petr

  fun(model1)
 
  but giving error as
 
  Error in if (ss == 1) { : missing value where TRUE/FALSE needed.
 
  please tell me how do i solve this problem.
 
 
 
[[alternative HTML version deleted]]
 
  __
  R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-
  guide.html and provide commented, minimal, self-contained,
  reproducible code.



Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny 
pouze jeho adresátům.
Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
svého systému.
Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
zpožděním přenosu e-mailu.

V případě, že je tento e-mail součástí obchodního jednání:
- vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a 
to z jakéhokoliv důvodu i bez uvedení důvodu.
- a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce 
s dodatkem či odchylkou.
- trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
dosažením shody na všech jejích náležitostech.
- odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně 
pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně 
osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi 
či osobě jím zastoupené známá.

This e-mail and any documents attached to it may be confidential and are 
intended only for its intended recipients.
If you received this e-mail by mistake, please immediately inform its sender. 
Delete the contents of this e-mail with all attachments and its copies from 
your system.
If you are not the intended recipient of this e-mail, you are not authorized to 
use, disseminate, copy or disclose this e-mail in any manner.
The sender of this e-mail shall not be liable for any

[R] How to calculate vif of each term of model in R?

2015-04-08 Thread Methekar, Pushpa (GE Transportation, Non-GE)




I am beginner in R doing modelling in R, I loaded excel sheet in R, i have 
chosen x elements and y elements then fitted model for linear and second order 
regression. Now I have both models. I am bit confused how to calculate vif for 
each term in model like

e.g model1-lm(y1~x1+x2+.x9) when I am using rms package then it's giving 
me like

vif(model1)

   x1 x2 x3 x4 x5 x6 x7

 6.679692   1.520271   1.667125   3.618439   4.931810   2.073879  13.870630

x8 x9

   220.969628 214.034135

now i want to compare each term with std vif as vif=10 and which will satisfy 
this condition i want to delete that term and update model1. i have done 
something like this

fun = function(model1) {

 for(i in 1:length(model1)){

  v=vif(model1)

 ss=any(v[i]=10)

if(ss==1){update(model1,.~.,-v[i])}

else{print(no update)}

 i-i+1

}



return(model1)

  }

fun(model1)

but giving error as

Error in if (ss == 1) { : missing value where TRUE/FALSE needed.

please tell me how do i solve this problem.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate correlation of a vector in R?

2014-11-02 Thread C W

Hi list,
I have trying to calculate the covariance/correlation of three elements.  I
have vector say,

v - c(700, 800, 1000)

I want to have a 3 by 3 correlation matrix, meaning cor(v1, c2), cor(v1,
c3), cor(v2, v3), etc...

So far I get,
 cor(v)
Error in cor(v) : supply both 'x' and 'y' or a matrix-like 'x'

 vvv - cbind(v, v, v)
 cor(vvv)
  v v v
v 1 1 1
v 1 1 1
v 1 1 1


I am calculating squared exponential kernel as seen here.
http://mlg.eng.cam.ac.uk/duvenaud/cookbook/index.html

Thanks a bunch,

Mike

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate correlation of a vector in R?

2014-11-02 Thread Jeff Newmiller

What is your question? The matrix form is probably what you are looking for, 
but you put the same vector in three times so if course it is all ones. I don't 
know what you expected to happen when you entered cor(v) since there is nothing 
to correlate if you only have one vector.

Please post in plain text as the footer and Posting Guide request.
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

On November 2, 2014 3:30:14 PM PST, C W tmrs...@gmail.com wrote:
Hi list,
I have trying to calculate the covariance/correlation of three
elements.  I
have vector say,

v - c(700, 800, 1000)

I want to have a 3 by 3 correlation matrix, meaning cor(v1, c2),
cor(v1,
c3), cor(v2, v3), etc...

So far I get,
 cor(v)
Error in cor(v) : supply both 'x' and 'y' or a matrix-like 'x'

 vvv - cbind(v, v, v)
 cor(vvv)
  v v v
v 1 1 1
v 1 1 1
v 1 1 1


I am calculating squared exponential kernel as seen here.
http://mlg.eng.cam.ac.uk/duvenaud/cookbook/index.html

Thanks a bunch,

Mike

   [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate correlation of a vector in R?

2014-11-02 Thread C W

Thanks, Jeff.  I had some misunderstanding.

So, I want to calculate the squared exponential of vector v

v = c(700, 800, 1029)

formula is:
k(x_i, x_j)=sigma^2 * exp(-1/(2*l^2) * (x_i - x_j) ^2)

where,
sigma=7, l=100


I used,
 v - c(700, 800, 1029)
 corr.matrix(cbind(v),scales=0.5)
 [,1] [,2] [,3]
[1,]100
[2,]010
[3,]001

the output should be covariance matrix = [49, 29.7, 0.02; 29.7, 49, 3.6;
0.2,  3.6, 49]



On Sun, Nov 2, 2014 at 7:04 PM, Jeff Newmiller jdnew...@dcn.davis.ca.us
wrote:

 What is your question? The matrix form is probably what you are looking
 for, but you put the same vector in three times so if course it is all
 ones. I don't know what you expected to happen when you entered cor(v)
 since there is nothing to correlate if you only have one vector.

 Please post in plain text as the footer and Posting Guide request.
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live
 Go...
   Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 ---
 Sent from my phone. Please excuse my brevity.

 On November 2, 2014 3:30:14 PM PST, C W tmrs...@gmail.com wrote:
 Hi list,
 I have trying to calculate the covariance/correlation of three
 elements.  I
 have vector say,
 
 v - c(700, 800, 1000)
 
 I want to have a 3 by 3 correlation matrix, meaning cor(v1, c2),
 cor(v1,
 c3), cor(v2, v3), etc...
 
 So far I get,
  cor(v)
 Error in cor(v) : supply both 'x' and 'y' or a matrix-like 'x'
 
  vvv - cbind(v, v, v)
  cor(vvv)
   v v v
 v 1 1 1
 v 1 1 1
 v 1 1 1
 
 
 I am calculating squared exponential kernel as seen here.
 http://mlg.eng.cam.ac.uk/duvenaud/cookbook/index.html
 
 Thanks a bunch,
 
 Mike
 
[[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate correlation of a vector in R?

2014-11-02 Thread Jeff Newmiller

k - sigma^2 * exp( -1/(2*l^2) * outer( v,v,FUN=function(x,y){(x-y)^2}))

but perhaps you should look at the e1071 package instead?
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

On November 2, 2014 4:19:56 PM PST, C W tmrs...@gmail.com wrote:
Thanks, Jeff.  I had some misunderstanding.

So, I want to calculate the squared exponential of vector v

v = c(700, 800, 1029)

formula is:
k(x_i, x_j)=sigma^2 * exp(-1/(2*l^2) * (x_i - x_j) ^2)

where,
sigma=7, l=100


I used,
 v - c(700, 800, 1029)
 corr.matrix(cbind(v),scales=0.5)
 [,1] [,2] [,3]
[1,]100
[2,]010
[3,]001

the output should be covariance matrix = [49, 29.7, 0.02; 29.7, 49,
3.6;
0.2,  3.6, 49]



On Sun, Nov 2, 2014 at 7:04 PM, Jeff Newmiller
jdnew...@dcn.davis.ca.us
wrote:

 What is your question? The matrix form is probably what you are
looking
 for, but you put the same vector in three times so if course it is
all
 ones. I don't know what you expected to happen when you entered
cor(v)
 since there is nothing to correlate if you only have one vector.

 Please post in plain text as the footer and Posting Guide request.

---
 Jeff NewmillerThe .   .  Go
Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live
 Go...
   Live:   OO#.. Dead: OO#.. 
Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#. 
rocks...1k

---
 Sent from my phone. Please excuse my brevity.

 On November 2, 2014 3:30:14 PM PST, C W tmrs...@gmail.com wrote:
 Hi list,
 I have trying to calculate the covariance/correlation of three
 elements.  I
 have vector say,
 
 v - c(700, 800, 1000)
 
 I want to have a 3 by 3 correlation matrix, meaning cor(v1, c2),
 cor(v1,
 c3), cor(v2, v3), etc...
 
 So far I get,
  cor(v)
 Error in cor(v) : supply both 'x' and 'y' or a matrix-like 'x'
 
  vvv - cbind(v, v, v)
  cor(vvv)
   v v v
 v 1 1 1
 v 1 1 1
 v 1 1 1
 
 
 I am calculating squared exponential kernel as seen here.
 http://mlg.eng.cam.ac.uk/duvenaud/cookbook/index.html
 
 Thanks a bunch,
 
 Mike
 
[[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate a numeric's digits count?

2014-10-24 Thread PO SU


Ok,  what i want is  find how many numbers after  . in a numeric ,and i don't 
know if there is already exists a function to do it( i wrote one by myself 
which will be showed later).
e.g.
1.234 has 3 numbers after . 
1 has 0 number
1.5342 has 4 numbers 
And i solved the above format using:
find-function(x)
{
str-as.character(x)
if(is.na(strsplit(str,\\.)[[1]][2])) return(0)
else return(nchar(strsplit(str,\\.)[[1]][2]))  
}

But when i  find(1.340)  i get 2 not 3. find(1.3400) will also get 2 not 4.
So,my question is how to implement the above needing? TKS.




--

PO SU
mail: desolato...@163.com 
Majored in Statistics from SJTU




At 2014-10-24 12:04:18, Jeff Newmiller jdnew...@dcn.davis.ca.us wrote:
I am baffled. I think those were English words but they didn't make any sense 
to me. Not was there a reproducible example to turn to. Can you try again?
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

On October 23, 2014 8:35:06 PM PDT, PO SU rhelpmaill...@163.com wrote:

Dear usRers,
  Now i want to cal ,e.g. 
 cal(1.234)  will get 3
 cal(1) will get 0
 cal(1.3045) will get 4
 But the difficult part is cal(1.3450) will get 4 not 3.
So, is there anyone happen to know the solution to this problem, or it
can't be solved in R, because 1.340 will always be transformed autolly
to 1.34?






--

PO SU
mail: desolato...@163.com 
Majored in Statistics from SJTU
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate a numeric's digits count?

2014-10-24 Thread Olivier Crouzet

Le Fri, 24 Oct 2014 14:11:08 +0800 (CST), PO SU a écrit :

 
 Ok,  what i want is  find how many numbers after  . in a numeric ,and
 i don't know if there is already exists a function to do it( i wrote
 one by myself which will be showed later). e.g. 1.234 has 3 numbers
 after . 1 has 0 number
 1.5342 has 4 numbers 
 And i solved the above format using:
 find-function(x)
 {
 str-as.character(x)
 if(is.na(strsplit(str,\\.)[[1]][2])) return(0)
  else return(nchar(strsplit(str,\\.)[[1]][2]))  
 }

It appears that your initial vector (x) is a numeric one that you
(obviously) need to transform to a character vector. So any number
it contains will have been shortened to its minimal representation when
you convert it to string. As far as I can tell (though I'm no expert in
number representation), you should work with x as a string vector from
the very beginning (which, to me, seems rather intuitive as the 0 in
1.230 is, really, only a string isn't it ?)

If you can see zeros when you print the numeric vector by e.g. 
print (x), that's simply because R displays them with a default
precision but they may contain many more digits... Try dput (x) before
str - as.character (x) and you will see what I mean.

Also, it may be usefull to look at what the following lines produce :
sprintf(%.10f, x)
sprintf(%.30f, x)

You will see that actually R stores many more digits that what you
think there are and that your 0 in 1.340 for example is,
really not single... except when it is basically a string.

Olivier.

 
 But when i  find(1.340)  i get 2 not 3. find(1.3400) will also get 2
 not 4. So,my question is how to implement the above needing? TKS.
 
 
 
 
 --
 
 PO SU
 mail: desolato...@163.com 
 Majored in Statistics from SJTU
 
 
 
 
 At 2014-10-24 12:04:18, Jeff Newmiller jdnew...@dcn.davis.ca.us
 wrote:
 I am baffled. I think those were English words but they didn't make
 any sense to me. Not was there a reproducible example to turn to.
 Can you try again?
 ---
 Jeff NewmillerThe .   .  Go
 Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.
 ##.#.  Live Go...
   Live:   OO#.. Dead: OO#..
  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.
 rocks...1k
 --- 
 Sent from my phone. Please excuse my brevity.
 
 On October 23, 2014 8:35:06 PM PDT, PO SU rhelpmaill...@163.com
 wrote:
 
 Dear usRers,
   Now i want to cal ,e.g. 
  cal(1.234)  will get 3
  cal(1) will get 0
  cal(1.3045) will get 4
  But the difficult part is cal(1.3450) will get 4 not 3.
 So, is there anyone happen to know the solution to this problem, or
 it can't be solved in R, because 1.340 will always be transformed
 autolly to 1.34?
 
 
 
 
 
 
 --
 
 PO SU
 mail: desolato...@163.com 
 Majored in Statistics from SJTU
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html and provide commented,
 minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate a numeric's digits count?

2014-10-24 Thread S Ellison

Dear Po Su,

All floating point numbers in R have exactly the same number of binary digits 
(53) and therefore the same number of decimal digits (15.95, as decimals aren't 
represented exactly in binary). If you want to find out how many decimal digits 
are after the decimal point, you could try subtracting floor(log10(x)) from 16.
Anything else is pretty much futile - to the point of nonsensical. You can 
never have exact representation of all decimal fractions in a binary computer, 
and once you understand that you can see that the number of decimal digits you 
get from any character representation depends only on how much you decide to 
round when converting to character format. That is essentially arbitrary, so 
any games you play with conversion to character are just telling you how many 
digits you decided to round each number to, not how many there were to start 
with.


S Ellison




***
This email and any attachments are confidential. Any use...{{dropped:8}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate a numeric's digits count?

2014-10-24 Thread Duncan Murdoch

On 23/10/2014, 11:35 PM, PO SU wrote:
 
 Dear usRers,
   Now i want to cal ,e.g. 
  cal(1.234)  will get 3
  cal(1) will get 0
  cal(1.3045) will get 4
  But the difficult part is cal(1.3450) will get 4 not 3.
 So, is there anyone happen to know the solution to this problem, or it can't 
 be solved in R, because 1.340 will always be transformed autolly to 1.34?
 

No, there's no way to do what you want unless you put quotes around the
number.  R parses 1.345 and 1.3450 as exactly the same thing, whereas
1.345 and 1.3450 are different.

Duncan Murdoch

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate a numeric's digits count?

2014-10-24 Thread David L Carlson

Where do these numbers come from? If they are calculated values, they are 
actually many decimal places longer than your examples. They are represented on 
your terminal with fewer decimals according to the setting of 
options(digits). 

For example:

 sqrt(2)*sqrt(2)
[1] 2
 sqrt(2)*sqrt(2) == 2  
[1] FALSE
# FAQ 7.31 Why doesn’t R think these numbers are equal?
 options(digits)
$digits
[1] 7
 options(digits=22)
 sqrt(2)*sqrt(2)
[1] 2.000444089

If the numbers were read from a plain text file and you are talking about how 
they are represented in the file, analyze them as character strings.

-
David L Carlson
Department of Anthropology
Texas AM University
College Station, TX 77840-4352

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of PO SU
Sent: Thursday, October 23, 2014 10:35 PM
To: R. Help
Subject: [R] how to calculate a numeric's digits count?


Dear usRers,
  Now i want to cal ,e.g. 
 cal(1.234)  will get 3
 cal(1) will get 0
 cal(1.3045) will get 4
 But the difficult part is cal(1.3450) will get 4 not 3.
So, is there anyone happen to know the solution to this problem, or it can't be 
solved in R, because 1.340 will always be transformed autolly to 1.34?






--

PO SU
mail: desolato...@163.com 
Majored in Statistics from SJTU
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate a numeric's digits count?

2014-10-23 Thread PO SU


Dear usRers,
  Now i want to cal ,e.g. 
 cal(1.234)  will get 3
 cal(1) will get 0
 cal(1.3045) will get 4
 But the difficult part is cal(1.3450) will get 4 not 3.
So, is there anyone happen to know the solution to this problem, or it can't be 
solved in R, because 1.340 will always be transformed autolly to 1.34?






--

PO SU
mail: desolato...@163.com 
Majored in Statistics from SJTU
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 a numeric's digits count?

2014-10-23 Thread Jeff Newmiller

I am baffled. I think those were English words but they didn't make any sense 
to me. Not was there a reproducible example to turn to. Can you try again?
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

On October 23, 2014 8:35:06 PM PDT, PO SU rhelpmaill...@163.com wrote:

Dear usRers,
  Now i want to cal ,e.g. 
 cal(1.234)  will get 3
 cal(1) will get 0
 cal(1.3045) will get 4
 But the difficult part is cal(1.3450) will get 4 not 3.
So, is there anyone happen to know the solution to this problem, or it
can't be solved in R, because 1.340 will always be transformed autolly
to 1.34?






--

PO SU
mail: desolato...@163.com 
Majored in Statistics from SJTU
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate moving average without using filter()?

2014-02-18 Thread C W

Thanks everyone.

For 5 point moving average,
filter(x, side=2, filter=rep(1/5, 5)), versus,
filter(x, side=2, filter=rep(1, 5)

Do they have the same effect, since the total needs to be 1.

Gabor  Rui: I am aware of the zoo package, I did not want to install a
package for one function.  Same reason for sos package.

David, thanks, that is what I am looking for.

Mike

On Mon, Feb 17, 2014 at 2:07 PM, Rui Barradas ruipbarra...@sapo.pt wrote:

 Hello,

 Many packages have a movind average function. For instance package
 forecast. Or

 library(sos)
 findFn(moving average)

 In your example, what you compute is not exactly a moving average, but in
 can be computed with something like the following.

 s - (seq_along(dat) - 1) %/% 3
 sapply(split(dat, s), mean)


 Hope this helps,

 Rui Barradas


 Em 17-02-2014 18:45, C W escreveu:

 Hi list,
 How do I calculate a moving average without using filter().  filter() does
 not seem to give weighted averages.

 I am looking into apply(), tapply,... But nothing moves.

 For example,

 dat-c(1:20)
 mean(dat[1:3])
 mean(dat[4:6])
 mean(dat[7:9])
 mean(dat[10:12])

 etc...

 I understand the point of apply is to avoid loops, how should I
 incorporate
 this idea into using an apply()?

 Thanks,
 Mike

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/
 posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate moving average without using filter()?

2014-02-17 Thread C W

Hi list,
How do I calculate a moving average without using filter().  filter() does
not seem to give weighted averages.

I am looking into apply(), tapply,... But nothing moves.

For example,

dat-c(1:20)
mean(dat[1:3])
mean(dat[4:6])
mean(dat[7:9])
mean(dat[10:12])

etc...

I understand the point of apply is to avoid loops, how should I incorporate
this idea into using an apply()?

Thanks,
Mike

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate moving average without using filter()?

2014-02-17 Thread Bert Gunter

There are a zillion answers to this, because your question is really:
How do I smooth a time series? So you can search on appropriate
keywords.

My answer is: don't use moving averages -- that's pathetically
ancient. ?loess is one among the zillions of alternatives you might
consider. Post on CV (stats.stackexchange.com) for other statistical
alternatives for time series smoothing.

Also, the understanding you expressed above is flawed. apply-type
constructs **are** (R-level) loops. So have you done your homework by
reading An Intro to R
(http://cran.r-project.org/doc/manuals/R-intro.pdf) or other web
tutorials? If not, please do so before posting here further.

Cheers,
Bert

-- Bert

Bert Gunter
Genentech Nonclinical Biostatistics
(650) 467-7374

Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom.
H. Gilbert Welch




On Mon, Feb 17, 2014 at 10:45 AM, C W tmrs...@gmail.com wrote:
 Hi list,
 How do I calculate a moving average without using filter().  filter() does
 not seem to give weighted averages.

 I am looking into apply(), tapply,... But nothing moves.

 For example,

 dat-c(1:20)
 mean(dat[1:3])
 mean(dat[4:6])
 mean(dat[7:9])
 mean(dat[10:12])

 etc...

 I understand the point of apply is to avoid loops, how should I incorporate
 this idea into using an apply()?

 Thanks,
 Mike

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate moving average without using filter()?

2014-02-17 Thread David Winsemius


On Feb 17, 2014, at 10:45 AM, C W wrote:

 Hi list,
 How do I calculate a moving average without using filter().  filter() does
 not seem to give weighted averages.
 
 I am looking into apply(), tapply,... But nothing moves.
 
 For example,
 
 dat-c(1:20)
 mean(dat[1:3])
 mean(dat[4:6])
 mean(dat[7:9])
 mean(dat[10:12])
 
 etc...
 
 I understand the point of apply is to avoid loops, how should I incorporate
 this idea into using an apply()?
 

Construct a vector for grouping and use tapply. Modulo division is a common 
method for achieving this. Sometimes the seq-function can be used if you adjust 
the length properly.

 tapply(dat, (0:(length(dat)-1))%/%3, mean)
   0123456 
 2.0  5.0  8.0 11.0 14.0 17.0 19.5 


tapply(dat, round(seq(1, (length(dat)/3),  len=length(dat))), mean)
   1234567 
 1.5  4.5  8.0 11.0 14.5 18.0 20.0 

The comment about weighting dos not seem to be exemplified in your example.

 Thanks,
 Mike
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, 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.

Re: [R] How to calculate moving average without using filter()?

2014-02-17 Thread Rui Barradas


Hello,

Many packages have a movind average function. For instance package 
forecast. Or


library(sos)
findFn(moving average)

In your example, what you compute is not exactly a moving average, but 
in can be computed with something like the following.


s - (seq_along(dat) - 1) %/% 3
sapply(split(dat, s), mean)


Hope this helps,

Rui Barradas


Em 17-02-2014 18:45, C W escreveu:

Hi list,
How do I calculate a moving average without using filter().  filter() does
not seem to give weighted averages.

I am looking into apply(), tapply,... But nothing moves.

For example,

dat-c(1:20)
mean(dat[1:3])
mean(dat[4:6])
mean(dat[7:9])
mean(dat[10:12])

etc...

I understand the point of apply is to avoid loops, how should I incorporate
this idea into using an apply()?

Thanks,
Mike

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate running 8-hour averages of temperature

2013-10-24 Thread arun

Hi,
You may try
?rollmean from library(zoo)
?SMA from library(TTR)
dat1 - structure(...
library(zoo)

res1 - rollmean(dat1[-c(1:5),2],8)
library(TTR)
res2 - SMA(dat1[-c(1:5),2],8)
 all.equal(res1,res2[!is.na(res2)])
#[1] TRUE



A.K.

Hi, 

How can I calculate running 8-hour averages of temperature 
values that start at specific date and time? In the example data set 
shown below I want to start the calculation on May 31, at 5pm (row 6). 
Reproducible data is provided. 

Thanks 
structure(list(date = structure(c(644148000, 644151600, 644155200, 
644158800, 644162400, 644166000, 644169600, 644173200, 644176800, 
644180400, 644184000, 644187600, 644191200, 644194800, 644198400, 
644202000, 644205600, 644209200, 644212800, 644216400, 64422, 
644223600, 644227200, 644230800, 644234400, 644238000, 644241600, 
644245200, 644248800, 644252400, 644256000), class = c(POSIXct, 
POSIXt), tzone = ), temp = c(17.7, 18.3, 18.9, 19.7, 19.5, 
19, 18.8, 18, 16.3, 15.1, 13.9, 12.7, 11.9, 10.6, 10, 8.7, 7.8, 
8.0, 9.2, 12.4, 14.5, 15.9, 17.4, 18.6, 18.3, 18.4, 18.5, 
20.4, 19.7, 18.9, 18.1)), .Names = c(date, temp), row.names = 3613:3643, 
class = data.frame)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate bioclim for table dataset in dismo package

2013-09-02 Thread Kristi Glover

Hi R Experts,
I was trying to develop model (bioclim, Domin) using a table data in 'dismo' 
package. The data I have is at table format instead of images (stack of raster 
images). I followed the procedures of  dismo package to calculate the bioclim 
but I could not figure it out how I can implement these procedures using table 
data. 

I have pasted an example how I did it, but it did not work. if some one has 
done it before, would you mind to suggest me how I can implement bioclim of 
'dismo' using table data?
I really appropriate your help.
Sincerely
KG


#
library(dismo)
dd-structure(list(long = c(-75.747, -106.84, -105.27, -104.64, -103.71, 
-102.72, -101.79, -100.8, -105.69, -104.67, -103.71, -102.69, 
-101.66, -100.71, -99.685, -98.656, -97.627, -111.82, -110.2, 
-109.55, -106.88, -106.84, -105.85, -104.82, -103.84, -102.8, 
-101.83, -100.79, -99.742, -98.695, -97.744, -112.46, -111.8, 
-110.79, -108.41, -107.74), lat = c(19.792, 21.576, 21.347, 21.243, 
21.078, 21.185, 20.99, 21.067, 22.621, 22.763, 22.594, 22.704, 
22.797, 22.581, 22.642, 22.685, 22.711, 24.153, 23.44, 23.385, 
23.4, 23.702, 23.565, 23.71, 23.542, 23.654, 23.454, 23.533, 
23.596, 23.64, 23.377, 24.796, 24.459, 24.399, 24.503, 24.429
), sp1 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L), env1 = c(182.0037689, 163.3301239, 
443.0214233, 1240.155273, 1774.867432, 1909.528809, 2153.244141, 
1878.018433, 909.1315308, 1980.898438, 2271.118896, 2122.672852, 
2033.41626, 1534.658447, 828.4759522, 222.3117523, 23.4761219, 
79.52642822, 322.994751, 273.6637268, 35.5085907, 516.3795776, 
2205.419434, 2118.727539, 2178.901123, 1995.210083, 2048.075928, 
1824.84021, 1043.509888, 184.4526062, 12.6928978, 20.87172508, 
126.5344544, 258.8852844, 3.28964257, 140.3287811), env2 = c(1134L, 
550L, 2111L, 2523L, 2156L, 1209L, 1107L, 2605L, 3176L, 2490L, 
1360L, 801L, 1118L, 1484L, 2730L, 1309L, 104L, 197L, 2033L, 1339L, 
567L, 2694L, 2708L, 1806L, 1344L, 912L, 1377L, 2346L, 3265L, 
989L, 69L, 428L, 764L, 896L, 100L, 1521L), env3 = c(24.533, 24.928, 
24.707, 21.052, 21.318, 18.428, 19.041, 17.743, 24.371, 19.689, 
16.879, 16.528, 16.901, 18.015, 20.648, 25.31, 24.308, 22.528, 
22.912, 22.001, 25.097, 25.391, 19.154, 14.943, 17.143, 16.898, 
16.891, 17.563, 15.63, 24.354, 24.088, 22.527, 22.528, 22.126, 
25.317, 25.643), env4 = c(0.047969, 0.003469, 0.003385, 0.002253, 
0.000791, 0.001834, 0.008016, 0.009262, 0.003934, 0.002322, 0.00061, 
0.000799, 6.4e-05, 0, 0, 0.000107, 0.003151, 0.018915, 0.015077, 
0.004554, 0.003499, 0.002705, 0.003507, 0.001173, 0.000149, 0.000308, 
0, 0, 0, 0.00074, 0.002845, 0.017047, 0.018915, 0.017111, 0.002417, 
0.002668)), .Names = c(long, lat, sp1, env1, env2, 
env3, env4), class = data.frame, row.names = c(NA, -36L))

bioclim.model-bioclim(dd[,4:7],dd[,3])

I got the following message
Error in function (classes, fdef, mtable)  : 
  unable to find an inherited method for function bioclim, for signature 
data.frame, integer





  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate the average values of each row in a matrix

2013-07-19 Thread Elaine Kuo

Hello,

I have a matrix (class matrix) composed of GridCell (row and column).
The matrix value is the beta diversity index value between two grids.

Now I would like to get the average value of each GridCell.
Please kindly advise how to make the calculation.
Thank you.

Elaine

The matrix looks like (cited from Michael Friendly)
I would like to get the average value of each color.


  Obs  stim   RPur   Red   Yel   Gy1   Gy2  Green  Blue  BlP  Pur1
 Pur2

  1  RPur . . . . .  . .
 . . .
  2  Red11.5. . . .  . .
 . . .
  3  Yel13.1   6.0. . .  . .
 . . .
  4  Gy112.6   7.9   6.2. .  . .
 . . .
  5  Gy210.6   8.4   8.4   5.2.  . .
 . . .
  6  Green  10.6   9.4   9.9   6.5   4.1 . .
 . . .
  7  Blue   10.8  10.2  10.3   8.8   7.06.4.
 . . .
   8  BlP 7.3  11.3  12.7  11.2  10.49.9   4.2
  . . .
   9  Pur15.4  11.5  12.9  11.7  10.89.4   8.4
 4.5. .
  10  Pur25.0  11.5  10.7  10.2  10.6   10.1   8.1
 6.43 .

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate the average values of each row in a matrix

2013-07-19 Thread R. Michael Weylandt michael.weyla...@gmail.com



On Jul 19, 2013, at 20:19, Elaine Kuo elaine.kuo...@gmail.com wrote:

 Hello,
 
 I have a matrix (class matrix) composed of GridCell (row and column).
 The matrix value is the beta diversity index value between two grids.
 
 Now I would like to get the average value of each GridCell.
 Please kindly advise how to make the calculation.
 Thank you.
 

Perhaps the rowMeans() function?

MW

 Elaine
 
 The matrix looks like (cited from Michael Friendly)
 I would like to get the average value of each color.
 
 
  Obs  stim   RPur   Red   Yel   Gy1   Gy2  Green  Blue  BlP  Pur1
 Pur2
 
  1  RPur . . . . .  . .
 . . .
  2  Red11.5. . . .  . .
 . . .
  3  Yel13.1   6.0. . .  . .
 . . .
  4  Gy112.6   7.9   6.2. .  . .
 . . .
  5  Gy210.6   8.4   8.4   5.2.  . .
 . . .
  6  Green  10.6   9.4   9.9   6.5   4.1 . .
 . . .
  7  Blue   10.8  10.2  10.3   8.8   7.06.4.
 . . .
   8  BlP 7.3  11.3  12.7  11.2  10.49.9   4.2
  . . .
   9  Pur15.4  11.5  12.9  11.7  10.89.4   8.4
 4.5. .
  10  Pur25.0  11.5  10.7  10.2  10.6   10.1   8.1
 6.43 .
 
[[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate the mean in a period of time?

2013-05-22 Thread arun

HI GG,
I thought you were referring about the solution that I sent today.  Sorry, I 
didn't check it.
TRy this:


dat1- read.table(text=
patient_id  number  responsed_at  delais  scores1   scores2   
scores3
 1  1    2010-05-26 NA   2.6   
0.5    0.7  
 1  2 2010-07-07 42    2.5  
 NA   NA   
 1   3    2010-08-14 38    2.3  
 NA   NA  
 1   4    2010-10-01  48    2.5 
  0.7   0.6
 1   5    2010-12-01  61    2.5 
  NA   NA
2    1    2011-07-19  NA   2.5  
 0.8    0.5
2 2   2011-09-22  65    2.6 
  NA    NA
2    3 2011-10-26  34    2.7
  NA    NA
3    1 2011-07-17 NA   2.8  
0.5 0.6
3    2    2011-10-30  103   2.6 
 NA NA
3 3  2011-12-23    54    2.5
  NA NA
,sep=,header=TRUE,stringsAsFactors=FALSE)
dat1$idx-with(dat1,ifelse(is.na(delais)|delais45 
delais20, 1,ifelse(delais60 
delais=45,2,ifelse(delais=90  delais=60,3,NA

#deleted lines
dat1[as.logical(with(dat1,ave(idx,patient_id,FUN=function(x) 
cumsum(is.na(x),-8]
#   patient_id number responsed_at delais scores1 scores2 scores3
#10  3  2   2011-10-30    103 2.6  NA  NA
#11  3  3   2011-12-23 54 2.5  NA  NA

 nrow(dat1[as.logical(with(dat1,ave(idx,patient_id,FUN=function(x) 
cumsum(is.na(x),-8])
#[1] 2
#or
table(with(dat1,ave(idx,patient_id,FUN=function(x) cumsum(is.na(x)[2]
#1 
#2 

#added lines
datSub-dat1[!as.logical(with(dat1,ave(idx,patient_id,FUN=function(x) 
cumsum(is.na(x),]
#count
 sum(with(datSub,idx[idx1  !is.na(idx)])-1)
#[1] 5

A.K.




From: GUANGUAN LUO guanguan...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Wednesday, May 22, 2013 11:30 AM
Subject: Re: how to calculate the mean in a period of time?



Hi Arun,
I meant to say number of observations you have deleted.
In the table, you have deleted some lines which have a delais  90 days, and 
added some lines copied from the precedent lines. Can I calculate the number of 
lines that you have deleted and the number of lines that you have added?

GG


2013/5/22 arun smartpink...@yahoo.com

Hi GG,
It is not clear what you meant by number of observations you have detected.  
Is it the number of observations in the initial dataset for each patient_id? 
Number of observations added also means what?  In the final result, we are 
taking the means for each 3 observations keeping the t=0  along with the 
first cluster for each Patient_id.






From: GUANGUAN LUO guanguan...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Wednesday, May 22, 2013 11:13 AM

Subject: Re: how to calculate the mean in a period of time?



Hello A.K,
for this reconstructed table, can i calculate the number of the observations 
you have delected and the number of the observations that you have added?
Thanks

GG


2013/5/17 arun smartpink...@yahoo.com



Hi,
Try this:

 dat1$idx-with(dat1,ifelse(is.na(delais)|delais45 
delais20, 1,ifelse(delais60 
delais=45,2,ifelse(delais=90  delais=60,3,NA
dat1$idx1-c(dat1$idx[-head(dat1$idx,1)],1)

library(zoo)
res1-do.call(rbind,lapply(split(dat1,dat1$patient_id),function(x) 
{x$idx[as.logical(cumsum(is.na(x$idx)))]-NA; x1-x[!is.na(x$idx),]; 
x1[,6:8]-na.locf(x1[,6:8]);x1$idx1[is.na(x1$idx1)]-1; 
x2-x1[rep(seq_len(nrow(x1)),x1$idx1),]; 
x2$delais[duplicated(x2$delais,fromLast=FALSE)]-0; 
x2$t-seq(0,nrow(x2)-1,1);x2[,-c(8:9)]}))
 row.names(res1)- 1:nrow(res1)
 res2- res1[,c(1:3,8,4:7)]
 res2

#   patient_id number responsed_at t delais scores1 scores2 scores3
#1   1  1   2010-05-26 0 NA 2.6 0.5 0.7
#2   1  2   2010-07-07 1 42 2.5 0.5 0.7
#3   1  3   2010-08-14 2 38 2.3 0.5 0.7
#4   1  3   2010-08-14 3  0 2.3 0.5 0.7

#5   1  4   2010-10-01 4 48 2.5 0.7 0.6
#6   1  4   2010-10-01 5  0 2.5 0.7 0.6
#7   1  4   2010-10-01 6  0 2.5 0.7 0.6

#8   1  5   2010-12-01 7 61 2.5 0.7 0.6
#9   2  1   2011-07-19 0 NA 2.5 0.8 0.5
#10  2  1   2011-07-19 1  0 2.5 0.8 0.5
#11  2  1   2011-07-19 2  0 2.5 0.8 0.5

#12  2  2   2011-09-22 3 65 2.6 0.8 0.5
#13  2  3

Re: [R] how to calculate the mean in a period of time?

2013-05-17 Thread arun



Hi,
Try this:
 dat1$idx-with(dat1,ifelse(is.na(delais)|delais45  
delais20, 1,ifelse(delais60  
delais=45,2,ifelse(delais=90  delais=60,3,NA
dat1$idx1-c(dat1$idx[-head(dat1$idx,1)],1)

library(zoo)
res1-do.call(rbind,lapply(split(dat1,dat1$patient_id),function(x) 
{x$idx[as.logical(cumsum(is.na(x$idx)))]-NA; x1-x[!is.na(x$idx),]; 
x1[,6:8]-na.locf(x1[,6:8]);x1$idx1[is.na(x1$idx1)]-1; 
x2-x1[rep(seq_len(nrow(x1)),x1$idx1),]; 
x2$delais[duplicated(x2$delais,fromLast=FALSE)]-0; 
x2$t-seq(0,nrow(x2)-1,1);x2[,-c(8:9)]}))
 row.names(res1)- 1:nrow(res1)
 res2- res1[,c(1:3,8,4:7)]
 res2
#   patient_id number responsed_at t delais scores1 scores2 scores3
#1   1  1   2010-05-26 0 NA 2.6 0.5 0.7
#2   1  2   2010-07-07 1 42 2.5 0.5 0.7
#3   1  3   2010-08-14 2 38 2.3 0.5 0.7
#4   1  3   2010-08-14 3  0 2.3 0.5 0.7
#5   1  4   2010-10-01 4 48 2.5 0.7 0.6
#6   1  4   2010-10-01 5  0 2.5 0.7 0.6
#7   1  4   2010-10-01 6  0 2.5 0.7 0.6
#8   1  5   2010-12-01 7 61 2.5 0.7 0.6
#9   2  1   2011-07-19 0 NA 2.5 0.8 0.5
#10  2  1   2011-07-19 1  0 2.5 0.8 0.5
#11  2  1   2011-07-19 2  0 2.5 0.8 0.5
#12  2  2   2011-09-22 3 65 2.6 0.8 0.5
#13  2  3   2011-10-26 4 34 2.7 0.8 0.5
#14  3  1   2011-07-17 0 NA 2.8 0.5 0.6
A.K.


From: GUANGUAN LUO guanguan...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Friday, May 17, 2013 9:33 AM
Subject: Re: how to calculate the mean in a period of time?



Hello， 
Thank you for your help
the lines added to the tables are the precedent lines but not the followed 
lines, if i just change x2-x1[rep(seq_len(nrow(x1)-1), is that ok? and so the 
delais should be changed too, isn't it?


GG


2013/5/17 arun smartpink...@yahoo.com

Hi,
No problem.

Arun







From: GUANGUAN LUO guanguan...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Friday, May 17, 2013 4:19 AM

Subject: Re: how to calculate the mean in a period of time?



Ah, yes, that is the wrong thing i have written. Thank you so much. the output 
which you have got is right.
Thanks a lot.
GG


2013/5/16 arun smartpink...@yahoo.com

Hi,
The output you showed is not clear especially the for the scores3,

  2                         2           2011-09-22    3     65            2.6 
      0.8            0.8
2                        3             2011-10-26   4      34            2.7  
   0.8            0.8
3                        1             2011-07-17    0     NA           2.8   
   0.5             0.6
In the input data, the scores3 column didn't had 0.8.


This is what I got:
dat1- read.table(text=

patient_id  number  responsed_at  delais  scores1   scores2   
scores3
 1  1    2010-05-26 NA   2.6  
 0.5    0.7  
 1  2 2010-07-07 42    2.5
   NA   NA   
 1   3    2010-08-14 38    2.3
   NA   NA  
 1   4    2010-10-01  48    2.5   
    0.7   0.6
 1   5    2010-12-01  61    2.5   
    NA   NA
2    1    2011-07-19  NA   2.5
   0.8    0.5
2 2   2011-09-22  65    2.6   
    NA    NA
2    3 2011-10-26  34    2.7  
    NA    NA
3    1 2011-07-17 NA   2.8
  0.5 0.6
3    2    2011-10-30  103   2.6   
   NA NA
3 3  2011-12-23    54    2.5  
    NA NA
,sep=,header=TRUE,stringsAsFactors=FALSE)

 dat1$idx-with(dat1,ifelse(is.na(delais)|delais45  delais20, 
1,ifelse(delais60  delais=45,2,ifelse(delais=90  delais=60,3,NA
library(zoo)
res-do.call(rbind,lapply(split(dat1,dat1$patient_id),function(x) 
{x$idx[as.logical(cumsum(is.na(x$idx)))]-NA; x1-x[!is.na(x$idx),]; 
x1[,6:8]-na.locf(x1[,6:8]);x2-x1[rep(seq_len(nrow(x1)),x1$idx),]; 
x2$delais[duplicated(x2$delais,fromLast=TRUE)]-0; 
x2$t-seq(0,nrow(x2)-1,1);x2}))


row.names(res)- 1:nrow(res)
 res1- res[,c(1:3,9,4:7)]
res1
#   patient_id number responsed_at t delais scores1 scores2 scores3
#1   1  1   2010-05-26 0 NA 2.6 0.5 0.7
#2   1  2   2010-07-07 1 42 2.5 0.5 0.7
#3   1  3   2010-08-14 2 38 2.3 0.5 0.7
#4   1  4

Re: [R] how to calculate the mean in a period of time?

2013-05-16 Thread arun

Hi,
The output you showed is not clear especially the for the scores3, 
  2                         2           2011-09-22    3     65            2.6   
    0.8            0.8
2                        3             2011-10-26   4      34            2.7    
 0.8            0.8
3                        1             2011-07-17    0     NA           2.8     
 0.5             0.6
In the input data, the scores3 column didn't had 0.8.


This is what I got:
dat1- read.table(text=
patient_id  number  responsed_at  delais  scores1   scores2   
scores3
 1  1    2010-05-26 NA   2.6   
0.5    0.7   
 1  2 2010-07-07 42    2.5  
 NA   NA    
 1   3    2010-08-14 38    2.3  
 NA   NA   
 1   4    2010-10-01  48    2.5 
  0.7   0.6
 1   5    2010-12-01  61    2.5 
  NA   NA
2    1    2011-07-19  NA   2.5  
 0.8    0.5
2 2   2011-09-22  65    2.6 
  NA    NA
2    3 2011-10-26  34    2.7
  NA    NA
3    1 2011-07-17 NA   2.8  
0.5 0.6
3    2    2011-10-30  103   2.6 
 NA NA
3 3  2011-12-23    54    2.5
  NA NA
,sep=,header=TRUE,stringsAsFactors=FALSE)

 dat1$idx-with(dat1,ifelse(is.na(delais)|delais45  delais20, 
1,ifelse(delais60  delais=45,2,ifelse(delais=90  delais=60,3,NA
library(zoo)
res-do.call(rbind,lapply(split(dat1,dat1$patient_id),function(x) 
{x$idx[as.logical(cumsum(is.na(x$idx)))]-NA; x1-x[!is.na(x$idx),]; 
x1[,6:8]-na.locf(x1[,6:8]);x2-x1[rep(seq_len(nrow(x1)),x1$idx),]; 
x2$delais[duplicated(x2$delais,fromLast=TRUE)]-0; 
x2$t-seq(0,nrow(x2)-1,1);x2}))


row.names(res)- 1:nrow(res)
 res1- res[,c(1:3,9,4:7)]
res1
#   patient_id number responsed_at t delais scores1 scores2 scores3
#1   1  1   2010-05-26 0 NA 2.6 0.5 0.7
#2   1  2   2010-07-07 1 42 2.5 0.5 0.7
#3   1  3   2010-08-14 2 38 2.3 0.5 0.7
#4   1  4   2010-10-01 3  0 2.5 0.7 0.6
#5   1  4   2010-10-01 4 48 2.5 0.7 0.6
#6   1  5   2010-12-01 5  0 2.5 0.7 0.6
#7   1  5   2010-12-01 6  0 2.5 0.7 0.6
#8   1  5   2010-12-01 7 61 2.5 0.7 0.6
#9   2  1   2011-07-19 0 NA 2.5 0.8 0.5
#10  2  2   2011-09-22 1  0 2.6 0.8 0.5
#11  2  2   2011-09-22 2  0 2.6 0.8 0.5
#12  2  2   2011-09-22 3 65 2.6 0.8 0.5
#13  2  3   2011-10-26 4 34 2.7 0.8 0.5
#14  3  1   2011-07-17 0 NA 2.8 0.5 0.6
A.K.


From: GUANGUAN LUO guanguan...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Thursday, May 16, 2013 12:05 PM
Subject: Re: how to calculate the mean in a period of time?



Hello, AK, 
now i have a problem really complicated for me,

Now my table is like this:

patient_id  number  responsed_at  delais  scores1   scores2   
scores3
 1  1    2010-05-26     NA   2.6   
0.5        0.7    
 1  2             2010-07-07         42            2.5  
 NA           NA         
 1                       3    2010-08-14 38    2.3  
 NA       NA    
 1   4    2010-10-01          48            2.5     
  0.7           0.6
 1                       5            2010-12-01          61            2.5     
  NA           NA
2                        1            2011-07-19          NA           2.5      
 0.8            0.5
2                         2           2011-09-22          65            2.6     
  NA            NA
2                        3             2011-10-26          34            2.7    
  NA            NA
3                        1             2011-07-17         NA           2.8      
0.5             0.6
3                        2            2011-10-30          103           2.6     
 NA             NA
3                         3              2011-12-23        54            2.5    
  NA             NA

explications: delais = the date of responsed_at  - the date of precedent 
responsed_at
scores1 is measured every month
scores2 and 3 are measured every three months

first thing is :   if the 20delais 45, this count one month of delais
if the 45=delais 60, this

Re: [R] How to calculate Hightest Posterior Density (HPD) of coeficients in a simple regression (lm) in R?

2013-05-08 Thread Ben Bolker

Richard Asturia richard.asturia at gmail.com writes:

 
 Hi!
 
 I am trying to calculate HPD for the coeficients of regression models
 fitted with lm or lmrob in R, pretty much in the same way that can be
 accomplished by the association of mcmcsamp and HPDinterval functions for
 multilevel models fitted with lmer. Can anyone point me in the right
 direction on which packages/how to implement this?
 
 Thanks for your time!
 
 R.
 

 Hmmm. 
  At least for lm(), if the assumptions of the model are met
then the sampling distribution of the parameters should be
multivariate normal, so with a flat prior the posterior distributions
should be symmetric and equivalent to the sampling distributions of
the parameters -- so I think that the highest 95% posterior density
interval should be equivalent to classical frequentist confidence
intervals [see confint()].

  You might be interested in the bayeslm() function from the arm
package.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 Hightest Posterior Density (HPD) of coeficients in a simple regression (lm) in R?

2013-05-08 Thread Richard Asturia

Thanks a lot for your reply!

That is interesting function. And you are completely right. Of course it
may be an overkill to use hpd with lm... But anyway, I am precisely
interested in comparing CI and HPD for a model fitted with lm and then do
the same comparison for the same model fitted with lmrob. I played a little
bit with the bayeslm() you've suggested and I think it pretty much solves
my issue for lm. But what about doing something similar to lmrob? In this
case, as far as I know, there is no similar function. If this is right,
some manual solution should be needed, i.e. some code-implementation
using MCMC.

Any ideas on which funcionts/code could help me doing so for lmrob fitted
objects?

Thanks again!

R.


2013/5/8 Ben Bolker bbol...@gmail.com

 Richard Asturia richard.asturia at gmail.com writes:

 
  Hi!
 
  I am trying to calculate HPD for the coeficients of regression models
  fitted with lm or lmrob in R, pretty much in the same way that can be
  accomplished by the association of mcmcsamp and HPDinterval functions for
  multilevel models fitted with lmer. Can anyone point me in the right
  direction on which packages/how to implement this?
 
  Thanks for your time!
 
  R.
 

  Hmmm.
   At least for lm(), if the assumptions of the model are met
 then the sampling distribution of the parameters should be
 multivariate normal, so with a flat prior the posterior distributions
 should be symmetric and equivalent to the sampling distributions of
 the parameters -- so I think that the highest 95% posterior density
 interval should be equivalent to classical frequentist confidence
 intervals [see confint()].

   You might be interested in the bayeslm() function from the arm
 package.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate the mean in a period of time?

2013-05-07 Thread arun

Hi,
Your question is still not clear.
May be this helps:

dat2- read.table(text=
patient_id  t scores
1  0    1.6
1  1    2.6
1  2 2.2
1  3 1.8
2  0  2.3
2   2 2.5
2  4  2.6
2   5 1.5
,sep=,header=TRUE)

library(plyr)
 dat2New-ddply(dat2,.(patient_id),summarize,t=seq(min(t),max(t)))
 res-join(dat2New,dat2,type=full)
res1-do.call(rbind,lapply(split(res,res$patient_id),function(x) 
{x1-x[x$t!=0,];do.call(rbind,lapply(split(x1,((x1$t-1)%/%3)+1),function(y) 
{y1-if(any(y$t==1)) rbind(x[x$t==0,],y) else y; 
data.frame(patient_id=unique(y1$patient_id),scores=mean(y1$scores,na.rm=TRUE))})
 ) }))
 row.names(res1)-1:nrow(res1)
res1$period-with(res1,ave(patient_id,patient_id,FUN=seq))
 res1
#  patient_id scores period
#1  1   2.05  1
#2  2   2.40  1
#3  2   2.05  2


A.K.




From: GUANGUAN LUO guanguan...@gmail.com
To: arun smartpink...@yahoo.com 
Sent: Tuesday, May 7, 2013 11:29 AM
Subject: Re: how to calculate the mean in a period of time?



Yes ， as you have said, probably , it's not continuous.


2013/5/7 arun smartpink...@yahoo.com

Hi,
Your question is not clear.  You mentioned to calculate the mean of 3 months, 
but infact you added the scores for t=0,1,2,3 as first 3 months, then possibly 
4,5,6 as the next.  So, it is not exactly three months.  Isn't it?


Dear R experts,
sorry to trouble you again.
My data is like this now :
patient_id      t         scores
1                      0                1.6
1                      1                2.6
1                      2                 2.2
1                      3                 1.8
2                      0                  2.3
2                       2                 2.5
2                      4                  2.6
2                       5                 1.5

I want to calculate the mean of period of 3 months, just get a table like this

patient_id     period     scores
1                            1           2.05                      
(1.6+2.6+2.2+1.8)/4
2                            1               2.4                     
(2.3+2.5)/2
2                            2               2.05                    
(2.6+1.5)/2

thank you in avance


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate Hightest Posterior Density (HPD) of coeficients in a simple regression (lm) in R?

2013-05-07 Thread Richard Asturia

Hi!

I am trying to calculate HPD for the coeficients of regression models
fitted with lm or lmrob in R, pretty much in the same way that can be
accomplished by the association of mcmcsamp and HPDinterval functions for
multilevel models fitted with lmer. Can anyone point me in the right
direction on which packages/how to implement this?

Thanks for your time!

R.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate average of each column

2013-04-10 Thread Ye Lin

Hey All,

I have a large dataset and I want to calculate the average of each column
then return a new dataset.

Here is my question: I dont know if there is a function that can allow me
to calculate the average every 60 records of data in the whole dataset, and
return a new data frame. Not sure if I have to divide the dataset first for
every 60, then do the mean or can i directly do that.

thanks for your help!


cici

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate average of each column

2013-04-10 Thread arun

Hi,
TRy this:
set.seed(52)
dat1- as.data.frame(matrix(sample(1:40,100*60,replace=TRUE), nrow=600))

lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),nrow)
#$`1`
#[1] 60

#$`2`
#[1] 60

#$`3`
#[1] 60

res-lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),colMeans)
 res[1:2]
#$`1`
#  V1   V2   V3   V4   V5   V6   V7   V8 
#20.95000 20.5 20.55000 21.1 20.1 18.3 21.1 20.5 
#  V9  V10 
#20.86667 18.7 

#$`2`
#  V1   V2   V3   V4   V5   V6   V7   V8 
#22.16667 19.85000 19.7 22.26667 18.8 19.9 18.85000 20.46667 
#  V9  V10 
#17.81667 18.51667 

A.K.



- Original Message -
From: Ye Lin ye...@lbl.gov
To: r-help@r-project.org
Cc: 
Sent: Wednesday, April 10, 2013 1:46 PM
Subject: [R] how to calculate average of each column

Hey All,

I have a large dataset and I want to calculate the average of each column
then return a new dataset.

Here is my question: I dont know if there is a function that can allow me
to calculate the average every 60 records of data in the whole dataset, and
return a new data frame. Not sure if I have to divide the dataset first for
every 60, then do the mean or can i directly do that.

thanks for your help!


cici

    [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate average of each column

2013-04-10 Thread Ye Lin

Thanks so much! it works fine!

But when I convert the result to a data frame, it mess up.

so i have:

res - data.frame(res)

then i get:

  X1   X2   X3   X4   X5   X6
V1  20.95000 22.16667 21.9 21.1 20.35000 19.36667
V2  20.5 19.85000 20.15000 17.5 19.3 20.88333
V3  20.55000 19.7 19.36667 22.4 20.06667 20.86667
V4  21.1 22.26667 20.2 18.16667 19.25000 19.7
V5  20.1 18.8 19.68333 21.45000 22.08333 22.7
V6  18.3 19.9 21.36667 21.2 23.0 21.5
V7  21.1 18.85000 19.55000 19.51667 22.25000 21.86667
V8  20.5 20.46667 21.6 20.3 19.15000 17.5
V9  20.86667 17.81667 20.5 19.65000 21.21667 21.96667
V10 18.7 18.51667 21.98333 20.1 22.15000 21.85000
  X7   X8   X9  X10
V1  19.08333 20.56667 18.9 22.4
V2  22.06667 23.28333 22.38333 19.56667
V3  20.75000 20.81667 19.41667 21.05000
V4  20.41667 20.36667 18.5 20.16667
V5  22.76667 21.16667 21.7 17.41667
V6  20.8 20.1 22.8 19.85000
V7  21.65000 21.0 19.2 18.95000
V8  19.85000 19.91667 20.68333 19.86667
V9  20.81667 21.08333 22.3 22.7
V10 16.81667 21.26667 19.25000 21.5

how can i return a final output like this:

with original column name

10 rows of each average record

Thanks!




On Wed, Apr 10, 2013 at 1:13 PM, arun smartpink...@yahoo.com wrote:

 Hi,
 TRy this:
 set.seed(52)
 dat1- as.data.frame(matrix(sample(1:40,100*60,replace=TRUE), nrow=600))

 lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),nrow)
 #$`1`
 #[1] 60

 #$`2`
 #[1] 60

 #$`3`
 #[1] 60

 res-lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),colMeans)
  res[1:2]
 #$`1`
 #  V1   V2   V3   V4   V5   V6   V7   V8
 #20.95000 20.5 20.55000 21.1 20.1 18.3 21.1 20.5
 #  V9  V10
 #20.86667 18.7

 #$`2`
 #  V1   V2   V3   V4   V5   V6   V7   V8
 #22.16667 19.85000 19.7 22.26667 18.8 19.9 18.85000 20.46667
 #  V9  V10
 #17.81667 18.51667

 A.K.



 - Original Message -
 From: Ye Lin ye...@lbl.gov
 To: r-help@r-project.org
 Cc:
 Sent: Wednesday, April 10, 2013 1:46 PM
 Subject: [R] how to calculate average of each column

 Hey All,

 I have a large dataset and I want to calculate the average of each column
 then return a new dataset.

 Here is my question: I dont know if there is a function that can allow me
 to calculate the average every 60 records of data in the whole dataset, and
 return a new data frame. Not sure if I have to divide the dataset first for
 every 60, then do the mean or can i directly do that.

 thanks for your help!


 cici

 [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate average of each column

2013-04-10 Thread Rui Barradas


Hello,

You should provide us with a data example, like the posting guide asks 
you to.

Without one, you could adapt the following example to your case


# Make up some data
dat - data.frame(X = rnorm(200), Y = rnorm(200))

# Divide into subsets of 60 rows each and compute the col means
grp - rep(1:(1 + nrow(dat) / 60), each = 60)[seq_len(nrow(dat))]
do.call(rbind, lapply(split(dat, grp), colMeans))


Hope this helps,

Rui Barradas

Em 10-04-2013 18:46, Ye Lin escreveu:

Hey All,

I have a large dataset and I want to calculate the average of each column
then return a new dataset.

Here is my question: I dont know if there is a function that can allow me
to calculate the average every 60 records of data in the whole dataset, and
return a new data frame. Not sure if I have to divide the dataset first for
every 60, then do the mean or can i directly do that.

thanks for your help!


cici

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate average of each column

2013-04-10 Thread Ye Lin

That is very helpful Rui, thanks so much!


On Wed, Apr 10, 2013 at 1:54 PM, Rui Barradas ruipbarra...@sapo.pt wrote:

 Hello,

 You should provide us with a data example, like the posting guide asks you
 to.
 Without one, you could adapt the following example to your case


 # Make up some data
 dat - data.frame(X = rnorm(200), Y = rnorm(200))

 # Divide into subsets of 60 rows each and compute the col means
 grp - rep(1:(1 + nrow(dat) / 60), each = 60)[seq_len(nrow(dat))]
 do.call(rbind, lapply(split(dat, grp), colMeans))


 Hope this helps,

 Rui Barradas

 Em 10-04-2013 18:46, Ye Lin escreveu:

 Hey All,

 I have a large dataset and I want to calculate the average of each column
 then return a new dataset.

 Here is my question: I dont know if there is a function that can allow me
 to calculate the average every 60 records of data in the whole dataset,
 and
 return a new data frame. Not sure if I have to divide the dataset first
 for
 every 60, then do the mean or can i directly do that.

 thanks for your help!


 cici

 [[alternative HTML version deleted]]

 __**
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/**listinfo/r-helphttps://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/**
 posting-guide.html http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate average of each column

2013-04-10 Thread arun

Hi,

Is this what you need?
as.data.frame(do.call(rbind,res))
# V1   V2   V3   V4   V5   V6   V7   V8
#1  20.95000 20.5 20.55000 21.1 20.1 18.3 21.1 20.5
#2  22.16667 19.85000 19.7 22.26667 18.8 19.9 18.85000 20.46667
#3  21.9 20.15000 19.36667 20.2 19.68333 21.36667 19.55000 21.6
#4  21.1 17.5 22.4 18.16667 21.45000 21.2 19.51667 20.3
#5  20.35000 19.3 20.06667 19.25000 22.08333 23.0 22.25000 19.15000
#6  19.36667 20.88333 20.86667 19.7 22.7 21.5 21.86667 17.5
#7  19.08333 22.06667 20.75000 20.41667 22.76667 20.8 21.65000 19.85000
#8  20.56667 23.28333 20.81667 20.36667 21.16667 20.1 21.0 19.91667
#9  18.9 22.38333 19.41667 18.5 21.7 22.8 19.2 20.68333
#10 22.4 19.56667 21.05000 20.16667 17.41667 19.85000 18.95000 19.86667
  #   V9  V10
#1  20.86667 18.7
#2  17.81667 18.51667
#3  20.5 21.98333
#4  19.65000 20.1
#5  21.21667 22.15000
#6  21.96667 21.85000
#7  20.81667 16.81667
#8  21.08333 21.26667
#9  22.3 19.25000
#10 22.7 21.5
A.K.




 From: Ye Lin ye...@lbl.gov
To: arun smartpink...@yahoo.com 
Cc: R help r-help@r-project.org 
Sent: Wednesday, April 10, 2013 4:50 PM
Subject: Re: [R] how to calculate average of each column
 

Thanks so much! it works fine!


But when I convert the result to a data frame, it mess up.


so i have:


res - data.frame(res)


then i get:

X1   X2   X3   X4   X5   X6
V1  20.95000 22.16667 21.9 21.1 20.35000 19.36667
V2  20.5 19.85000 20.15000 17.5 19.3 20.88333
V3  20.55000 19.7 19.36667 22.4 20.06667 20.86667
V4  21.1 22.26667 20.2 18.16667 19.25000 19.7
V5  20.1 18.8 19.68333 21.45000 22.08333 22.7
V6  18.3 19.9 21.36667 21.2 23.0 21.5
V7  21.1 18.85000 19.55000 19.51667 22.25000 21.86667
V8  20.5 20.46667 21.6 20.3 19.15000 17.5
V9  20.86667 17.81667 20.5 19.65000 21.21667 21.96667
V10 18.7 18.51667 21.98333 20.1 22.15000 21.85000 X7   X8   X9  
X10
V1  19.08333 20.56667 18.9 22.4
V2  22.06667 23.28333 22.38333 19.56667
V3  20.75000 20.81667 19.41667 21.05000
V4  20.41667 20.36667 18.5 20.16667
V5  22.76667 21.16667 21.7 17.41667
V6  20.8 20.1 22.8 19.85000
V7  21.65000 21.0 19.2 18.95000
V8  19.85000 19.91667 20.68333 19.86667
V9  20.81667 21.08333 22.3 22.7
V10 16.81667 21.26667 19.25000 21.5


how can i return a final output like this:


with original column name

10 rows of each average record


Thanks!





On Wed, Apr 10, 2013 at 1:13 PM, arun smartpink...@yahoo.com wrote:

Hi,
TRy this:
set.seed(52)
dat1- as.data.frame(matrix(sample(1:40,100*60,replace=TRUE), nrow=600))

lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),nrow)
#$`1`
#[1] 60

#$`2`
#[1] 60

#$`3`
#[1] 60

res-lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),colMeans)
 res[1:2]
#$`1`
#  V1   V2   V3   V4   V5   V6   V7   V8
#20.95000 20.5 20.55000 21.1 20.1 18.3 21.1 20.5
#  V9  V10
#20.86667 18.7

#$`2`
#  V1   V2   V3   V4   V5   V6   V7   V8
#22.16667 19.85000 19.7 22.26667 18.8 19.9 18.85000 20.46667
#  V9  V10
#17.81667 18.51667

A.K.




- Original Message -
From: Ye Lin ye...@lbl.gov
To: r-help@r-project.org
Cc:
Sent: Wednesday, April 10, 2013 1:46 PM
Subject: [R] how to calculate average of each column

Hey All,

I have a large dataset and I want to calculate the average of each column
then return a new dataset.

Here is my question: I dont know if there is a function that can allow me
to calculate the average every 60 records of data in the whole dataset, and
return a new data frame. Not sure if I have to divide the dataset first for
every 60, then do the mean or can i directly do that.

thanks for your help!


cici

    [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] how to calculate average of each column

2013-04-10 Thread arun

Hi,
Try this:
set.seed(52)
dat1- as.data.frame(matrix(sample(c(1:40,NA),100*60,replace=TRUE), nrow=600))
 
res1-as.data.frame(do.call(rbind,lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 
60)+1),function(x) colMeans(x,na.rm=TRUE
 res1
# V1   V2   V3   V4   V5   V6   V7   V8
#1  21.20339 21.1 20.64407 20.94828 20.22034 17.91379 21.38983 20.0
#2  21.48214 19.5 19.41379 22.13793 18.53448 20.4 18.94915 19.77193
#3  22.11864 20.03448 19.55932 20.61667 19.41379 21.49153 20.08333 21.44828
#4  21.7 17.9 22.32759 18.7 21.61017 20.94828 19.13793 20.32203
#5  20.91667 19.37288 20.16949 18.12500 22.05172 23.01724 22.17241 19.22034
#6  19.50847 21.05085 20.70690 20.16667 22.22807 21.36207 21.63793 17.13793
#7  19.5 22.7 20.98305 20.96667 23.06780 20.98305 21.83051 19.91525
#8  20.36207 23.55932 20.94915 20.47458 21.25424 19.94828 19.31481 20.01695
#9  18.62069 22.03509 19.50847 18.95000 21.19298 23.01695 19.6 20.44828
#10 21.96491 20.1 21.61667 20.65000 17.76667 20.25000 18.28070 19.68966
 #    V9  V10
#1  21.01695 18.84746
#2  17.46552 18.9
#3  20.69492 22.6
#4  19.05263 20.30508
#5  21.7 22.40678
#6  21.86207 21.3
#7  20.81034 17.25000
#8  21.5 21.45763
#9  22.18966 19.7
#10 22.31579 20.58929
A.K.




 From: Ye Lin ye...@lbl.gov
To: arun smartpink...@yahoo.com 
Sent: Wednesday, April 10, 2013 6:02 PM
Subject: Re: [R] how to calculate average of each column
 

Hey A.K,


I want to exclude the missing values in the table when do the col mean, and I 
code like this:

res-lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 
60)+1),colMeans(dat1,na.rm=TRUE)


then i get this error message:

colMeans(dat1, na.rm = TRUE)' is not a function, character or symbol


How can I tell R to omit the NAs automatically and do the mean? For example, if 
there is 4 out of 10 NAs in one column, it will calculate mean as sum of the 
remaining 6 values and divide by 6.


Thanks!






On Wed, Apr 10, 2013 at 1:13 PM, arun smartpink...@yahoo.com wrote:

Hi,
TRy this:
set.seed(52)
dat1- as.data.frame(matrix(sample(1:40,100*60,replace=TRUE), nrow=600))

lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),nrow)
#$`1`
#[1] 60

#$`2`
#[1] 60

#$`3`
#[1] 60

res-lapply(split(dat1,((seq_len(nrow(dat1))-1)%/% 60)+1),colMeans)
 res[1:2]
#$`1`
#  V1   V2   V3   V4   V5   V6   V7   V8
#20.95000 20.5 20.55000 21.1 20.1 18.3 21.1 20.5
#  V9  V10
#20.86667 18.7

#$`2`
#  V1   V2   V3   V4   V5   V6   V7   V8
#22.16667 19.85000 19.7 22.26667 18.8 19.9 18.85000 20.46667
#  V9  V10
#17.81667 18.51667

A.K.




- Original Message -
From: Ye Lin ye...@lbl.gov
To: r-help@r-project.org
Cc:
Sent: Wednesday, April 10, 2013 1:46 PM
Subject: [R] how to calculate average of each column

Hey All,

I have a large dataset and I want to calculate the average of each column
then return a new dataset.

Here is my question: I dont know if there is a function that can allow me
to calculate the average every 60 records of data in the whole dataset, and
return a new data frame. Not sure if I have to divide the dataset first for
every 60, then do the mean or can i directly do that.

thanks for your help!


cici

    [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] how to calculate left kronecker product?

2013-02-23 Thread Michael Friendly

For an application, I have formulas defined in terms of a left Kronecker 
product of matrices,

A,B, meaning
A \otimes_L B = {A * B[i,j]}  -- matrix on the left multiplies each 
element on the right.


The standard kronecker() function is the right Kronecker product,
A \otimes_R B = {A[i,j] * B}  -- matrix on the right multiplies each 
element on the left.


The example below shows the result of kronecker() and what I want, but
kronecker() is now defined in generic S4 methods, and  I can't see
how to use more basic functions to get the result I want.  Or, alternatively
how to transform the result of kronecker() to give my wanted.

 test code ---
A - matrix(1:4, 2,2, dimnames=list(c(a1, a2), c(a1,a2)))
B - diag(2)
dimnames(B) - list(c(b1, b2), c(b1,b2))
# standard, right kronecker product: each A[i,j] * B
kronecker(A, B, make.dimnames=TRUE)

# left kronecker product: A * each B[i,j]
wanted - rbind(
cbind(A * B[1,1], A*B[1,2]),
cbind(A * B[2,1], A*B[2,2]))
rownames(wanted) - colnames(wanted) - paste(rep(c(b1, b2), 
each=2), rownames(wanted), sep=:)

wanted

 R output 
 A - matrix(1:4, 2,2, dimnames=list(c(a1, a2), c(a1,a2)))
 B - diag(2)
 dimnames(B) - list(c(b1, b2), c(b1,b2))
 # standard, right kronecker product: each A[i,j] * B
 kronecker(A, B, make.dimnames=TRUE)
  a1:b1 a1:b2 a2:b1 a2:b2
a1:b1 1 0 3 0
a1:b2 0 1 0 3
a2:b1 2 0 4 0
a2:b2 0 2 0 4

 # left kronecker product: A * each B[i,j]
 wanted - rbind(
+ cbind(A * B[1,1], A*B[1,2]),
+ cbind(A * B[2,1], A*B[2,2]))

 rownames(wanted) - colnames(wanted) - paste(rep(c(b1, b2), 
each=2), rownames(wanted), sep=:)

 wanted
  b1:a1 b1:a2 b2:a1 b2:a2
b1:a1 1 3 0 0
b1:a2 2 4 0 0
b2:a1 0 0 1 3
b2:a2 0 0 2 4




--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.  Chair, Quantitative Methods
York University  Voice: 416 736-2100 x66249 Fax: 416 736-5814
4700 Keele StreetWeb:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 left kronecker product?

2013-02-23 Thread Berend Hasselman


On 23-02-2013, at 20:46, Michael Friendly frien...@yorku.ca wrote:

 For an application, I have formulas defined in terms of a left Kronecker 
 product of matrices,
 A,B, meaning
 A \otimes_L B = {A * B[i,j]}  -- matrix on the left multiplies each element 
 on the right.
 
 The standard kronecker() function is the right Kronecker product,
 A \otimes_R B = {A[i,j] * B}  -- matrix on the right multiplies each element 
 on the left.
 
 The example below shows the result of kronecker() and what I want, but
 kronecker() is now defined in generic S4 methods, and  I can't see
 how to use more basic functions to get the result I want.  Or, alternatively
 how to transform the result of kronecker() to give my wanted.
 
  test code ---
 A - matrix(1:4, 2,2, dimnames=list(c(a1, a2), c(a1,a2)))
 B - diag(2)
 dimnames(B) - list(c(b1, b2), c(b1,b2))
 # standard, right kronecker product: each A[i,j] * B
 kronecker(A, B, make.dimnames=TRUE)
 
 # left kronecker product: A * each B[i,j]
 wanted - rbind(
 cbind(A * B[1,1], A*B[1,2]),
 cbind(A * B[2,1], A*B[2,2]))
 rownames(wanted) - colnames(wanted) - paste(rep(c(b1, b2), each=2), 
 rownames(wanted), sep=:)
 wanted
 
  R output 
  A - matrix(1:4, 2,2, dimnames=list(c(a1, a2), c(a1,a2)))
  B - diag(2)
  dimnames(B) - list(c(b1, b2), c(b1,b2))
  # standard, right kronecker product: each A[i,j] * B
  kronecker(A, B, make.dimnames=TRUE)
  a1:b1 a1:b2 a2:b1 a2:b2
 a1:b1 1 0 3 0
 a1:b2 0 1 0 3
 a2:b1 2 0 4 0
 a2:b2 0 2 0 4
 
  # left kronecker product: A * each B[i,j]
  wanted - rbind(
 + cbind(A * B[1,1], A*B[1,2]),
 + cbind(A * B[2,1], A*B[2,2]))
 
  rownames(wanted) - colnames(wanted) - paste(rep(c(b1, b2), each=2), 
  rownames(wanted), sep=:)
  wanted
  b1:a1 b1:a2 b2:a1 b2:a2
 b1:a1 1 3 0 0
 b1:a2 2 4 0 0
 b2:a1 0 0 1 3
 b2:a2 0 0 2 4
 

How about

kronecker(B, A, make.dimnames=TRUE)

Berend

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 left kronecker product?

2013-02-23 Thread Michael Friendly


Ugh. So simple!  Thanks, Berend.

On 2/23/2013 2:58 PM, Berend Hasselman wrote:

On 23-02-2013, at 20:46, Michael Friendly frien...@yorku.ca wrote:


For an application, I have formulas defined in terms of a left Kronecker 
product of matrices,
A,B, meaning
A \otimes_L B = {A * B[i,j]}  -- matrix on the left multiplies each element on 
the right.

The standard kronecker() function is the right Kronecker product,
A \otimes_R B = {A[i,j] * B}  -- matrix on the right multiplies each element on 
the left.

The example below shows the result of kronecker() and what I want, but
kronecker() is now defined in generic S4 methods, and  I can't see
how to use more basic functions to get the result I want.  Or, alternatively
how to transform the result of kronecker() to give my wanted.

 test code ---
A - matrix(1:4, 2,2, dimnames=list(c(a1, a2), c(a1,a2)))
B - diag(2)
dimnames(B) - list(c(b1, b2), c(b1,b2))
# standard, right kronecker product: each A[i,j] * B
kronecker(A, B, make.dimnames=TRUE)

# left kronecker product: A * each B[i,j]
wanted - rbind(
cbind(A * B[1,1], A*B[1,2]),
cbind(A * B[2,1], A*B[2,2]))
rownames(wanted) - colnames(wanted) - paste(rep(c(b1, b2), each=2), 
rownames(wanted), sep=:)
wanted

 R output 

A - matrix(1:4, 2,2, dimnames=list(c(a1, a2), c(a1,a2)))
B - diag(2)
dimnames(B) - list(c(b1, b2), c(b1,b2))
# standard, right kronecker product: each A[i,j] * B
kronecker(A, B, make.dimnames=TRUE)

  a1:b1 a1:b2 a2:b1 a2:b2
a1:b1 1 0 3 0
a1:b2 0 1 0 3
a2:b1 2 0 4 0
a2:b2 0 2 0 4

# left kronecker product: A * each B[i,j]
wanted - rbind(

+ cbind(A * B[1,1], A*B[1,2]),
+ cbind(A * B[2,1], A*B[2,2]))

rownames(wanted) - colnames(wanted) - paste(rep(c(b1, b2), each=2), 
rownames(wanted), sep=:)
wanted

  b1:a1 b1:a2 b2:a1 b2:a2
b1:a1 1 3 0 0
b1:a2 2 4 0 0
b2:a1 0 0 1 3
b2:a2 0 0 2 4

How about

kronecker(B, A, make.dimnames=TRUE)

Berend





--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.  Chair, Quantitative Methods
York University  Voice: 416 736-2100 x66249 Fax: 416 736-5814
4700 Keele StreetWeb:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] How to calculate the moving average for binary files?

2013-02-18 Thread Jonsson

I have 12 binary (raster) files   
https://echange-fichiers.inra.fr/get?k=k3M2jatJyHy65Cs99G4 .
I would like to calculate the moving average for the 12 values for each
pixel in the 12 files.

For a simple vector we can get a moving average by using this :
  
 x - c(1,2,3,NA,NA,4,6,5,6,4,2,5)
movingmean - rollapply(x, 3, FUN = mean, na.rm = T,fill=NA)
now I want to do the same but with binary files and I tried:

files   - list.files(C:final-2010, *.envi, full.names = TRUE)
results - list()
for (.files in files) {
# read in the 12 files as a vector of numbers 
# that we take the average of
x - do.call(rbind,(lapply(.files, readBin  , double() , 
 size = 4 ,n =1440 * 720 , signed = T)))
# take the moving average across the 12 values 
# from the 12 files for each pixel
results[[length(results) + 1L]] - rollapply(x, 3, FUN = mean,na.rm
= T)
}
 
But got this error:
 
Error in seq.default(start.at, NROW(data), by = by) : 
 wrong sign in 'by' argument



--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-calculate-the-moving-average-for-binary-files-tp4658986.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] How to calculate the moving average for binary files?

2013-02-18 Thread David Winsemius


On Feb 18, 2013, at 12:28 PM, Jonsson wrote:

 I have 12 binary (raster) files   
 https://echange-fichiers.inra.fr/get?k=k3M2jatJyHy65Cs99G4 .
 I would like to calculate the moving average for the 12 values for each
 pixel in the 12 files.
 
 For a simple vector we can get a moving average by using this :
 
 x - c(1,2,3,NA,NA,4,6,5,6,4,2,5)
movingmean - rollapply(x, 3, FUN = mean, na.rm = T,fill=NA)
 now I want to do the same but with binary files and I tried:
 
files   - list.files(C:final-2010, *.envi, full.names = TRUE)
results - list()
for (.files in files) {
# read in the 12 files as a vector of numbers 
# that we take the average of
x - do.call(rbind,(lapply(.files, readBin  , double() , 
 size = 4 ,n =1440 * 720 , signed = T)))
# take the moving average across the 12 values 
# from the 12 files for each pixel
results[[length(results) + 1L]] - rollapply(x, 3, FUN = mean,na.rm
 = T)
}
 
 But got this error:
 
Error in seq.default(start.at, NROW(data), by = by) : 
 wrong sign in 'by' argument

Cross-posting on StackOverflow and Rhelp is deprecated in the Posting Guide.

http://stackoverflow.com/questions/14944824/how-to-calculate-the-moving-averag-for-several-binary-files-in-r

-- 

David Winsemius
Alameda, CA, 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.

Re: [R] How to calculate monthly average from daily files in R?

2013-01-18 Thread Pascal Oettli


Hello,

I guess that you are working on Windows. So, I don't know whether the 
following is pertinent or not.


On Linux, the first file is 'ET100.bin', not 'ET1.bin'. Thus, the first 
'monthly' mean is calculated between April, 10 and May, 10.


You can try the following script. When I read your data, the code for 
the missing values is -999, not -32765, maybe due to the machine. 
You can manage your own missing value code by modifying missval.


#-

missval - -999

files - paste0(C:\\New folder (4)\\New folder\\ET,1:365,.bin)

mm - factor(rep(1:12, c(31,28,31,30,31,30,31,31,30,31,30,31)))
files.group - split(files, mm)

monthly - NULL
for(ff in files.group){
  x - do.call(rbind,(lapply(ff, readBin, double(), size=4, n=360*720, 
signed=T)))

  x[x==missval] - NA
  x - apply(x,2,mean)
  monthly - rbind(monthly, x)
  rm(x)
}
dim(monthly)

#-

HTH,
Pascal

Le 18/01/2013 02:37, Jonsson a écrit :

I have 365 binary files:
https://echange-fichiers.inra.fr/get?k=oy3CN1yV1Um7ouRWm2U   ,I want to
calculate the monthly average. So from the 365 files, I will get 12 files.I
would like also to tell R not to take into account the no-data value
(-32765).for example, for the first month, there are 31 records: 3 of these
records has the value -32765,I want R to take the average of the rest
records(28  records) and so on with all months.

This code will take the average of every 30 files(any idea on how to make it
according to number of days in a month?and not to take into account the
no-data values)

files- list.files(C:\\New folder (4)\\New folder,
*.bin,full.names=TRUE)
# assume that we want to take the average of every 30 files
 files.group- split(files , rep(seq_along(files), each = 30,
length =length(files)))
   results- list()
  for (.files in files.group){
# read in the 30 files as a vector of numbers that you take
the average of
  x- do.call(rbind,(lapply(.files, readBin  , double() , size =
4 ,n =360 * 720 , signed =T)))
   ## take the means across the 30 files
  results[[length(results) + 1L]]- colMeans(x)}
   close(x)
  for (i in seq_along(results)){
 fileName - sprintf(C:/New folder/glo_%d.flt, i)
 writeBin(as.double(results[[i]]), fileName, size = 4)}



--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-calculate-monthly-average-from-daily-files-in-R-tp4655869.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

1 2 3 4 >

1 - 100 of 346 matches

Mail list logo