Re: [R] Extract time and state of charge (Start and End) and Count
Hi Roslina, I think you have changed the code as "bc_start" in my code is "BCStartTime" in yours. When I run the attached code, I get a data frame "hourly_SoC" that looks right, and a matrix "result" (hour by SoC) that checks against the data frame. I have tried to comment the code so that you an see what I am doing. Jim On Wed, Jul 20, 2022 at 1:18 AM roslinazairimah zakaria wrote: > > Hi Jim, > > I tried to run your code and got this error. > > > # get the temporal order of observations > > obs_order <- order(c(as.numeric(dt$BCStartTime),as.numeric(dt$BCStopTime))) > Warning messages: > 1: In order(c(as.numeric(dt$BCStartTime), as.numeric(dt$BCStopTime))) : > NAs introduced by coercion > 2: In order(c(as.numeric(dt$BCStartTime), as.numeric(dt$BCStopTime))) : > NAs introduced by coercion > > numeric_time<-c(as.numeric(dt$BCStartTime),as.numeric(dt$BCStopTime))[obs_order] > Warning messages: > 1: NAs introduced by coercion > 2: NAs introduced by coercion > > nobs<-diff(range(numeric_time))/3600 > > # find the linear approximation of charge state by hours > > hourly_SoC <- approx(numeric_time, > + > c(dt$Starting_SoC_of_12,dt$Ending_SoC_of_12)[obs_order],n=nobs) > Error in approx(numeric_time, c(dt$Starting_SoC_of_12, > dt$Ending_SoC_of_12)[obs_order], : > need at least two non-NA values to interpolate > __ 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] Need panel plots by each state, comprising three columns in R
Hello, I need to draw a panel plot by each importing the attached csv data (sample panel plot attached as well). Each state should contain three columns- wt, avg, and observed, may be represented as legend or x-axis, and y-axis should show values. I tried doing it in ggplot in R but couldn't succeed. Please help plotting this panel plot. Thanks, Ranjeet __ 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] Extract time and state of charge (Start and End) and Count
Hello,
There was a bug in the way I copied your data to my R session,
hence the NA's.
Here is a tidyverse way of doing what you want. Its output matches the
expected output in your last post. The column names don't start at zero
because there was no Starting_SoC_of_12 equal to 0.
library(tidyverse)
result <- dt_2014 %>%
mutate(Hour = lubridate::hour(BatteryChargeStartDate)) %>%
group_by(Hour, Starting_SoC_of_12, Ending_SoC_of_12) %>%
mutate(Count = n()) %>%
ungroup() %>%
arrange(Hour, Starting_SoC_of_12, Ending_SoC_of_12) %>%
pivot_wider(
id_cols = Hour,
names_from = c(Starting_SoC_of_12, Ending_SoC_of_12),
names_sep = "-",
values_from = Count,
values_fill = 0L
)
i <- str_order(names(result)[-1], numeric = TRUE)
result <- cbind(result[1], result[-1][i])
result
# Hour 1-11 2-10 2-11 4-4 4-12 5-8 8-12
#17000 00 10
#28000 01 01
#3 15000 10 00
#4 16110 00 00
#5 18001 00 01
#6 21000 01 01
Hope this helps,
Rui Barradas
Às 16:10 de 19/07/2022, roslinazairimah zakaria escreveu:
Hi Rui,
I try to run your code, but all data became NA. Not sure why...
# these columns need to be fixed
cols <- c("BatteryChargeStartDate", "BatteryChargeStopDate")
dt[cols] <- lapply(dt[cols], \(x) sub("\n", " ", x))
# use package lubridate to coerce to a datetime class
library(lubridate)
dt <- lapply(dt, lubridate::dmy_hm)
dt
dt[cols] <- lapply(dt[cols], lubridate::dmy_hm)
h <- lubridate::hour(dt[["BatteryChargeStartDate"]])
aggregate(Starting_SoC_of_12 ~ h, dt, length)
$BCStartTime
[1] 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
[33] 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
[65] 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
[97] 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
[129] 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
[161] 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
[193] 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
[225] 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
[257] 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
[289] NA NA NA NA NA NA NA NA
[929] 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
[961] 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
[993] NA NA NA NA NA NA NA NA
[ reached 'max' / getOption("max.print") -- omitted 34418 entries ]
dt[cols] <- lapply(dt[cols], lubridate::dmy_hm)
Warning messages:
1: All formats failed to parse. No formats found.
2: All formats failed to parse. No formats found.
h <- lubridate::hour(dt[["BatteryChargeStartDate"]])
aggregate(Starting_SoC_of_12 ~ h, dt, length)
Error in aggregate.data.frame(lhs, mf[-1L], FUN = FUN, ...) :
no rows to aggregate
On Tue, Jul 19, 2022 at 4:53 PM roslinazairimah zakaria <
roslina...@gmail.com> wrote:
Hi Rui,
Yes, I would like to count for each hour, how many in the state of charge
start 0 and SOC 12, then SOC 1 and SOC 12 and so on.
Thank you for your help.
On Tue, Jul 19, 2022 at 1:11 AM Rui Barradas wrote:
Hello,
I'm not sure I understand the problem. Do you want counts of how many
rows are there per hour?
# these columns need to be fixed
cols <- c("BatteryChargeStartDate", "BatteryChargeStopDate")
dt_2014[cols] <- lapply(dt_2014[cols], \(x) sub("\n", " ", x))
# use package lubridate to coerce to a datetime class
dt_2014[cols] <- lapply(dt_2014[cols], lubridate::dmy_hm)
h <- lubridate::hour(dt_2014[["BatteryChargeStartDate"]])
aggregate(Starting_SoC_of_12 ~ h, dt_2014, length)
It would be better if you post the expected output corresponding to the
posted data set.
Hope this helps,
Rui Barradas
Às 05:04 de 18/07/2022, roslinazairimah zakaria escreveu:
Dear all,
I have data of Battery Electric vehicle (BEV). I would like to extract
data
from every hour starting from 0.00 to 0.59, 1:00-1:59 for SOC(state of
charge) start to end.
Some examples:
I can extract data from SOC=0 and SOC=12
dt_2014[which(dt_2014$Starting_SoC_of_12==0 &
dt_2014$Ending_SoC_of_12==12),]
I can extract data from SOC=1 and SOC=12
dt_2014[which(dt_2014$Starting_SoC_of_12==1 &
dt_2014$Ending_SoC_of_12==12),]
and I would like to further categorise the data by hour and count how
many
cars from 0 state charge to 12 state charge at in that particular hour.
Thank you so much for any help given.
Some data
dput(dt_2014[1:10,])
structure(list(ï..CarID =
Re: [R] Extract time and state of charge (Start and End) and Count
Hi Jim,
I tried to run your code and got this error.
> # get the temporal order of observations
> obs_order <-
order(c(as.numeric(dt$BCStartTime),as.numeric(dt$BCStopTime)))
Warning messages:
1: In order(c(as.numeric(dt$BCStartTime), as.numeric(dt$BCStopTime))) :
NAs introduced by coercion
2: In order(c(as.numeric(dt$BCStartTime), as.numeric(dt$BCStopTime))) :
NAs introduced by coercion
>
numeric_time<-c(as.numeric(dt$BCStartTime),as.numeric(dt$BCStopTime))[obs_order]
Warning messages:
1: NAs introduced by coercion
2: NAs introduced by coercion
> nobs<-diff(range(numeric_time))/3600
> # find the linear approximation of charge state by hours
> hourly_SoC <- approx(numeric_time,
+
c(dt$Starting_SoC_of_12,dt$Ending_SoC_of_12)[obs_order],n=nobs)
Error in approx(numeric_time, c(dt$Starting_SoC_of_12,
dt$Ending_SoC_of_12)[obs_order], :
need at least two non-NA values to interpolate
On Tue, Jul 19, 2022 at 3:22 PM roslinazairimah zakaria <
roslina...@gmail.com> wrote:
> Jim,
>
> Thank you so much!
>
> On Tue, Jul 19, 2022 at 3:05 PM Jim Lemon wrote:
>
>> Hi Roslina,
>> The following gives you the state of charge for the vehicle in your
>> example data for each hour. This is approximate as your times are not
>> on even hours.
>>
>> # get the temporal order of observations
>>
>> obs_order<-order(c(as.numeric(dt_2014$bc_start),as.numeric(dt_2014$bc_stop)))
>>
>> numeric_time<-c(as.numeric(dt_2014$bc_start),as.numeric(dt_2014$bc_stop))[obs_order]
>> nobs<-diff(range(numeric_time))/3600
>> # find the linear approximation of charge state by hours
>> hourly_SoC<-approx(numeric_time,
>> c(dt_2014$Starting_SoC_of_12,dt_2014$Ending_SoC_of_12)[obs_order],n=nobs)
>>
>> To get the POSIX times:
>>
>> hourly_POSIX<-seq(dt_2014$bc_start[1],dt_2014$bc_stop,length.out=nobs)
>>
>> That will fill part of your table. If you want the state of charge by
>> hour regardless of day, you'll have to create a new "hour" variable
>> from, hourly_POSIX, then:
>>
>> mean_charge_by_hour<-by(hourly_SoC$hour,hourly_SoC$y,mean) #untested
>> since I don't know whether you want this
>>
>> Jim
>>
>> On Mon, Jul 18, 2022 at 2:04 PM roslinazairimah zakaria
>> wrote:
>> >
>> > Dear all,
>> >
>> > I have data of Battery Electric vehicle (BEV). I would like to extract
>> data
>> > from every hour starting from 0.00 to 0.59, 1:00-1:59 for SOC(state of
>> > charge) start to end.
>> >
>> > Some examples:
>> > I can extract data from SOC=0 and SOC=12
>> > dt_2014[which(dt_2014$Starting_SoC_of_12==0 &
>> > dt_2014$Ending_SoC_of_12==12),]
>> >
>> > I can extract data from SOC=1 and SOC=12
>> > dt_2014[which(dt_2014$Starting_SoC_of_12==1 &
>> > dt_2014$Ending_SoC_of_12==12),]
>> >
>> > and I would like to further categorise the data by hour and count how
>> many
>> > cars from 0 state charge to 12 state charge at in that particular hour.
>> >
>> > Thank you so much for any help given.
>> >
>> > Some data
>> > > dput(dt_2014[1:10,])
>> > structure(list(ï..CarID = c("GC10", "GC10", "GC10", "GC10", "GC10",
>> > "GC10", "GC10", "GC10", "GC10", "GC10"), BatteryChargeStartDate =
>> > c("16/2/2014 16:05",
>> > "16/2/2014 18:20", "17/2/2014 8:10", "18/2/2014 7:41", "18/2/2014
>> 15:36",
>> > "18/2/2014 16:36", "18/2/2014 21:26", "19/2/2014 8:57", "19/2/2014
>> 21:08",
>> > "20/2/2014 18:11"), BCStartTime = c("16:05", "18:20", "8:10",
>> > "7:41", "15:36", "16:36", "21:26", "8:57", "21:08", "18:11"),
>> > Year = c(2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
>> > 2014L, 2014L, 2014L), Month = c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
>> > 2L, 2L, 2L), Day = c(16L, 16L, 17L, 18L, 18L, 18L, 18L, 19L,
>> > 19L, 20L), BatteryChargeStopDate = c("16/2/2014 17:05", "16/2/2014
>> > 19:00",
>> > "17/2/2014 15:57", "18/2/2014 9:52", "18/2/2014 15:39", "18/2/2014
>> > 17:36",
>> > "19/2/2014 1:55", "19/2/2014 14:25", "20/2/2014 5:17", "20/2/2014
>> 23:20"
>> > ), BCStopTime = c("17:05", "19:00", "15:57", "9:52", "15:39",
>> > "17:36", "1:55", "14:25", "5:17", "23:20"), Year2 = c(2014L,
>> > 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L
>> > ), Month2 = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Day2 = c(16L,
>> > 16L, 17L, 18L, 18L, 18L, 19L, 19L, 20L, 20L), Starting_SoC_of_12 =
>> > c(1L,
>> > 2L, 4L, 5L, 4L, 2L, 8L, 8L, 4L, 8L), Ending_SoC_of_12 = c(11L,
>> > 11L, 12L, 8L, 4L, 10L, 12L, 12L, 12L, 12L)), row.names = c(NA,
>> > 10L), class = "data.frame")
>> >
>> >
>> > --
>> > *Roslinazairimah Zakaria*
>> > *Tel: +609-5492370; Fax. No.+609-5492766*
>> >
>> > *Email: roslina...@gmail.com *
>> > University Malaysia Pahang
>> > Lebuhraya Tun Razak, 26300 Gambang, Pahang, Malaysia
>> >
>> > [[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
[R] mgcv: relative risk from GAM with distributed lag
Dear list members,
Does anyone know how to obtain a relative risk/ risk ratio from a GAM
with a distributed lag model implemented in mgcv? I have a GAM
predicting daily deaths from time series data consisting of daily
temperature, humidity and rainfall. The GAM includes a distributed lag
model because deaths may occur over several days following a high heat
day.
What I'd like to do is compute (and plot) the relative risk
(accumulated across all lags) for a given temperature vs the
temperature at which the risk is lowest, with corresponding confidence
intervals. I am aware of the predict.gam function but am not sure if
and how it should be used in this case. (Additionally, I'd also like
to plot the relative risk for different lags separately).
I apologise if this seems trivial to some. (Actually, I hope it is,
because that might mean I get a solution!) I've been looking for
examples on how to do this, but found nothing so far. Suggestions
would be very much appreciated!
Below is a reproducible example with the GAM:
library(mgcv)
set.seed(3) # make reproducible example
simdat <- gamSim(1,400) # simulate data
g <- exp(simdat$f/5)
simdat$y <- rnbinom(g,size=3,mu=g) # negative binomial response var
simdat$time <- 1:400 # create time series
names(simdat) <- c("deaths", "temp", "humidity", "rain", "x3", "f",
"f0", "f1", "f2", "f3", "time")
# lag function based on Simon Wood (book 2017, p.349 and gamair
package documentation p.54
# https://cran.rstudio.com/web/packages/gamair/gamair.pdf)
lagard <- function(x,n.lag=7) {
n <- length(x); X <- matrix(NA,n,n.lag)
for (i in 1:n.lag) X[i:n,i] <- x[i:n-i+1]
X
}
# set up lag, temp, rain and humidity as 7-column matrices
# to create lagged variables - based on Simon Wood's example
dat <- list(lag=matrix(0:6,nrow(simdat),7,byrow=TRUE),
deaths=simdat$deaths, time = simdat$time)
dat$temp <- lagard(simdat$temp)
dat$rain <- lagard(simdat$rain)
dat$humidity <- lagard(simdat$humidity)
mod <- gam(deaths~s(time, k=70) + te(temp, lag, k=c(12, 4)) +
te(humidity, lag, k=c(12, 4)) + te(rain, lag, k=c(12, 4)), data = dat,
family = nb, method = 'REML', select = TRUE)
summary(mod)
plot(mod, scheme = 1)
plot(mod, scheme = 2)
Thanks for any suggestions you may have,
Jade
__
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.
