[R] Hypothesis Testing using Wald Criterion for two regression models with dummy variables
I have two models, controlled by dummy variables to see if the models can be combined into one model with similar intercepts and slopes. Has anyone tried to conduct this type of test in R. I am utilizing the econometric idea of hypothesis testing through the hypothesis of coincidence. I have tried to run an anova with test of Chisq, but I am not sure what the results are telling. In addition, I used the rms package with a lrm model in an anova test, again I am not sure what the results are telling me: Try 1 anova(H,Ha,test="Chi") Analysis of Variance Table Model 1: logload ~ logflow Model 2: logload ~ dummy + logflow Res.DfRSS Df Sum of Sq Pr(>Chi) 1 17 4.6742 2 16 2.6314 12.0428 0.0004245 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Try 2: fit<-lrm(logload~logflow+dummy) > anova(fit) Wald Statistics Response: logload Factor Chi-Square d.f. P logflow17.56 1<.0001 dummy 5.22 10.0224 TOTAL 18.03 20.0001 Can anyone help me with this? Thanks! -- View this message in context: http://r.789695.n4.nabble.com/Hypothesis-Testing-using-Wald-Criterion-for-two-regression-models-with-dummy-variables-tp4601582.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] Hypothesis Testing using Wald Criterion for two regression models with dummy variables
Peter- Maybe I have not articulately my problem clearly, I have had local help with the statistical part just trying to figure out how to correctly program this test. For clarity's sake, I have months worth of data, I want to potentially combine those months into four, shall we say seasons, that have close to the same behaviour. Therefore to do this, I am trying a monthly moving window to categorize these seasons. After talking to a couple water resources statisician's we decided the way to test if the months are different is through the use of hypothesis testing and a dummy variable. So I have one regression, Model A, that includes a zero in the dummy spot with the two months of data combined, then I have another regression, Model B, that includes the interaction term for the changes between the months (the intercept changes, using a 0 or 1 dummy variable). Now we discussed running a Wald testing, Chi squared, to test to see if the interaction term is of importance probability wise, can I do this utilizing anova? Does this make more sense? Then I will run another set of restricted and unrestricted models to account for potential differences in the mean (i.e the slope). Does this explain my problem better? Meredith -- View this message in context: http://r.789695.n4.nabble.com/Hypothesis-Testing-using-Wald-Criterion-for-two-regression-models-with-dummy-variables-tp4601582p4603260.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] Subtracting a matrix 1x28 from a scalar
Afternoon- I am trying to subtract a matrix, basically a vector of 28 values, each by the same number to account for differences in regression fitting. I am taking the 1x28 and minus it by the mean value of the matrix, each number. The result I receive is a 1X29 matrix. Does anyone know why the result has an extra value? > lffeb_march [1] 6.588926 7.663877 6.917706 6.824374 7.029973 6.549651 6.517671 6.070738 5.916202 [10] 6.993933 6.091310 7.313220 5.135798 6.762730 6.381816 7.999679 6.851185 5.799093 [19] 7.774856 6.956545 8.218787 8.218787 7.549609 7.222566 7.170120 6.173786 5.075174 [28] 4.700480 6.672033 5.198497 4.867534 7.170120 6.131226 6.802395 > meanf_fm<-mean(lffeb_march) > X3_fm<-(lffeb_march-meanf_fm) > X3_fm [1] -0.03784990 1.03710088 0.29092923 0.19759729 0.40319653 -0.07712564 -0.10910511 [8] -0.55603865 -0.71057432 0.36715660 -0.53546650 0.68644401 -1.49097794 0.13595313 [15] -0.24496036 1.37290220 0.22440855 -0.82768372 1.14807939 0.32976906 1.59201078 [22] 1.59201078 0.92283279 0.59578964 0.54334317 -0.45299027 -1.55160256 -1.92629601 [29] 0.04525657 -1.42827935 -1.75924193 0.54334317 -0.49554989 0.17561839 Thanks for the help Meredith -- View this message in context: http://r.789695.n4.nabble.com/Subtracting-a-matrix-1x28-from-a-scalar-tp4615590.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] Subtracting a matrix 1x28 from a scalar
* Disregard message please Afternoon- I am trying to subtract a matrix, basically a vector of 28 values, each by the same number to account for differences in regression fitting. I am taking the 1x28 and minus it by the mean value of the matrix, each number. The result I receive is a 1X29 matrix. Does anyone know why the result has an extra value? > lffeb_march [1] 6.588926 7.663877 6.917706 6.824374 7.029973 6.549651 6.517671 6.070738 5.916202 [10] 6.993933 6.091310 7.313220 5.135798 6.762730 6.381816 7.999679 6.851185 5.799093 [19] 7.774856 6.956545 8.218787 8.218787 7.549609 7.222566 7.170120 6.173786 5.075174 [28] 4.700480 6.672033 5.198497 4.867534 7.170120 6.131226 6.802395 > meanf_fm<-mean(lffeb_march) > X3_fm<-(lffeb_march-meanf_fm) > X3_fm [1] -0.03784990 1.03710088 0.29092923 0.19759729 0.40319653 -0.07712564 -0.10910511 [8] -0.55603865 -0.71057432 0.36715660 -0.53546650 0.68644401 -1.49097794 0.13595313 [15] -0.24496036 1.37290220 0.22440855 -0.82768372 1.14807939 0.32976906 1.59201078 [22] 1.59201078 0.92283279 0.59578964 0.54334317 -0.45299027 -1.55160256 -1.92629601 [29] 0.04525657 -1.42827935 -1.75924193 0.54334317 -0.49554989 0.17561839 Thanks for the help Meredith -- View this message in context: http://r.789695.n4.nabble.com/Subtracting-a-matrix-1x28-from-a-scalar-tp4615590p4615599.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] Automating R for Hypothesis Testing
R Users- I have been trying to automate a manual code that I have developed for calling in a .csv file, isolating certain rows and columns that correspond to specified months: something to the effect i=name.csv N=length(i$month) iphos1=0 iphos2=0 isphos3=0 for i=1,N if month=1 iphos1=iphos+1 iphos1(iphos1)=i an so on to call out the months into there own arrays (unless there is a way I can wrap it into the next automation) Next: I would like to run a simple linear regression combining each of the months 1 by 1: for instance I want to run a regression on a combined model from months 1 and 2 and a dummy model for 1 and 2, compare them using a Chi-sq distribution, if Chi-sq is less than the Critical value, we accept and go on to test another set of models with both 1 and 2. If it rejects, then we proceed to months 2 and 3. If we move on to the second set on months 1 and 2, and the critical value is accepted, I want to print an accept or reject and move on to months 2 and 3, until finally comparing months 12-1 at the end. Is there a way to loop or automate this in R? Thanks Meredith -- View this message in context: http://r.789695.n4.nabble.com/Automating-R-for-Hypothesis-Testing-tp4618653.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] Automating R for Hypothesis Testing
dput: http://r.789695.n4.nabble.com/file/n4620188/milwaukeephos.csv
milwaukeephos.csv
# Feb-march
> modelH_febmarch<-lm(llfeb_march~lffeb_march)
>modelHa_febmarch<-lm(llfeb_march~X1feb_mar+lffeb_march)
> anova(modelHa_febmarch)
> coefficients(modelH_febmarch)
(Intercept) lffeb_march
-2.4298901.172821
> coefficients(modelHa_febmarch)
(Intercept) X1feb_mar lffeb_march
-2.8957776 -0.5272793 1.3016303
> bres_fm<-matrix(c(-2.429890,0,1.172821),nrow=3)
> bunres_fm<-matrix(c(-2.8957776,-0.5272793,1.3016303),nrow=3)
>bfm<-t(bunres_fm-bres_fm)
> fmvect<-seq(1,1,length=34)
> X1a_febmar<-seq(0,0,length=9) # dummy variable step 1
> X1b_febmar<-seq(1,1,length=25) # dummy variable step 2
> X1feb_mar<-c(X1a_febmar,X1b_febmar) #dummy variable creation
# Test Stat Equation for Chisq
> fmxx<-cbind(fmvect,X1feb_mar,lffeb_march)
> tfmx<-t(fmxx)
> xcom_fm<-(tfmx %*% fmxx)
> xinv_fm<-ginv(xcom_fm)
> var_fm<-xinv_fm*0.307
> chi_fm<-bfm %*% var_fm %*% (bunres_fm-bres_fm)
> chi_fm # chisq value for recording
if less than CV move onto to slope modification
> modelH2_febmarch<-lm(llfeb_march~X3feb_march)
> modelH2a_febmarch<-lm(llfeb_march~X3feb_march+X4feb_march)
> anova(modelH2a_febmarch)
> coefficients(modelH2_febmarch) # get coefficients to make beta vectors for
> test
(Intercept) X3feb_march
5.3421301.172821
> coefficients(modelH2a_febmarch)
(Intercept) X3feb_march X4feb_march
5.2936263 1.0353752 0.2407557
# Test Stat
> bsres_fm<-matrix(c(5.342130,1.172821,0),nrow=3)
> bsunres_fm<-matrix(c(5.2936263,1.0353752,0.2407557),nrow=3)
> bsfm<-t(bsunres_fm-bsres_fm)
> #X matrix
> fmxs<-cbind(fmvect,X3feb_march,X4feb_march)
> tfmxs<-t(fmxs)
> xcoms_fm<-(tfmxs %*% fmxs)
> xinvs_fm<-ginv(xcoms_fm)
> var_fms<-xinvs_fm*0.341
> chi_fms<-bsfm %*% var_fms %*% (bsunres_fm-bsres_fm)
> chi_fms
# Record Chisq value
Does this help?
Here lffeb_march is the combination of Feb and March log flows
and llfeb_march is the combination of Feb and March log loads
X3: lffeb_march-mean(feb_march)
X4: X1*X3
Thanks
Rui Barradas wrote
>
> Hello,
>
> I'm not at all sure if I understand your problem. Does this describe it?
>
>
> test first model for months 1 and 2
> if test statistic less than critical value{
> test second model for months 1 and 2
> print results of the first and second tests? just one of them?
> }
> move on to months 2 and 3
> etc, until months 12 and 1
>
>
> Please post example data using dput(dataset).
> Just copy it's output and paste it in your post.
>
> And example code, what you're already doing.
> (Possibly simplified)
>
> Rui Barradas
>
>
> meredith wrote
>>
>> R Users-
>> I have been trying to automate a manual code that I have developed for
>> calling in a .csv file, isolating certain rows and columns that
>> correspond to specified months:
>> something to the effect
>> i=name.csv
>> N=length(i$month)
>> iphos1=0
>> iphos2=0
>> isphos3=0
>> for i=1,N
>> if month=1
>> iphos1=iphos+1
>> iphos1(iphos1)=i
>>
>> an so on to call out the months into there own arrays (unless there is a
>> way I can wrap it into the next automation)
>>
>> Next: I would like to run a simple linear regression combining each of
>> the months 1 by 1:
>> for instance I want to run a regression on a combined model from months 1
>> and 2 and a dummy model for 1 and 2, compare them using a Chi-sq
>> distribution, if Chi-sq is less than the Critical value, we accept and go
>> on to test another set of models with both 1 and 2. If it rejects, then
>> we proceed to months 2 and 3. If we move on to the second set on months
>> 1 and 2, and the critical value is accepted, I want to print an accept or
>> reject and move on to months 2 and 3, until finally comparing months 12-1
>> at the end.
>> Is there a way to loop or automate this in R?
>>
>> Thanks
>> Meredith
>>
>
Rui Barradas wrote
>
> Hello,
>
> I'm not at all sure if I understand your problem. Does this describe it?
>
>
> test first model for months 1 and 2
> if test statistic less than critical value{
> test second model for months 1 and 2
> print results of the first and second tests? just one of them?
> }
> move on to months 2 and 3
> etc, until months 12 and 1
>
>
> Please post example data using dput(dataset).
> Just copy it's output and paste it in your post.
>
> And example code, what you're already doing.
> (Possibly simplified)
>
> Rui Barradas
>
>
Re: [R] Automating R for Hypothesis Testing
Rui-
Thanks this definitely helps, just one quick question. How would you code
the values of chi-fm and chi-fms to change based on the degrees of freedom
of each model H(i)?
Meredith
Rui Barradas wrote
>
> Hello,
>
> Yes, it does help. Now we can see your data and what you're doing.
> What follows is a suggestion on what you could do, not full solution.
> (You forgot to say what X1 is, but I don't think it's important to
> understand the suggestion.)
> (If I'm wrong, say something.)
>
>
> milwaukeephos <- read.csv("milwaukeephos.csv", header=TRUE,
> stringsAsFactors=FALSE)
> # list of data.frames, one per month
> ls1 <- split(milwaukeephos, milwaukeephos$month)
>
> #- if you want to keep the models, not needed if you don't.
> # (yoy probably don't)
> modelH <- vector("list", 12)
> modelHa <- vector("list", 12)
> modelH2 <- vector("list", 12)
> modelH2a <- vector("list", 12)
> #- values to record, these are needed, create them beforehand.
> chi_fm <- numeric(12)
> chi_fms <- numeric(12)
> #
> seq_months <- c(1:12, 1) # wrap months around.
> for(i in 1:12){
> month_this <- seq_months[i]
> month_next <- seq_months[i + 1]
>
> lload <- c(ls1[[month_this]]$load_kg, ls1[[month_next]]$load_kg)
> lflow <- c(ls1[[month_this]]$flow, ls1[[month_next]]$flow)
> modelH[[i]] <- lm(lload ~ lflow)
> # If you don't want to keep the models, use modelH only
> # ( without [[i]] )
> # and do the same with X1
>
> # rest of your code for first test goes here
> chi_fm[i] <- bfm %*% var_fm %*% (bunres_fm - bres_fm)
>
> # and the same for the second test
> chi_fms[i] <- ...etc...
> }
>
>
> Hope this helps,
>
> Rui Barradas
>
>
> meredith wrote
>>
>> dput: http://r.789695.n4.nabble.com/file/n4620188/milwaukeephos.csv
>> milwaukeephos.csv
>>
>> # Feb-march
>>> modelH_febmarch<-lm(llfeb_march~lffeb_march)
>>>modelHa_febmarch<-lm(llfeb_march~X1feb_mar+lffeb_march)
>>> anova(modelHa_febmarch)
>>> coefficients(modelH_febmarch)
>> (Intercept) lffeb_march
>> -2.4298901.172821
>>> coefficients(modelHa_febmarch)
>> (Intercept) X1feb_mar lffeb_march
>> -2.8957776 -0.5272793 1.3016303
>>> bres_fm<-matrix(c(-2.429890,0,1.172821),nrow=3)
>>> bunres_fm<-matrix(c(-2.8957776,-0.5272793,1.3016303),nrow=3)
>>>bfm<-t(bunres_fm-bres_fm)
>>> fmvect<-seq(1,1,length=34)
>>> X1a_febmar<-seq(0,0,length=9) # dummy variable step 1
>>> X1b_febmar<-seq(1,1,length=25) # dummy variable step 2
>>> X1feb_mar<-c(X1a_febmar,X1b_febmar) #dummy variable creation
>> # Test Stat Equation for Chisq
>>> fmxx<-cbind(fmvect,X1feb_mar,lffeb_march)
>>> tfmx<-t(fmxx)
>>> xcom_fm<-(tfmx %*% fmxx)
>>> xinv_fm<-ginv(xcom_fm)
>>> var_fm<-xinv_fm*0.307
>>> chi_fm<-bfm %*% var_fm %*% (bunres_fm-bres_fm)
>>> chi_fm # chisq value for recording
>> if less than CV move onto to slope modification
>>> modelH2_febmarch<-lm(llfeb_march~X3feb_march)
>>> modelH2a_febmarch<-lm(llfeb_march~X3feb_march+X4feb_march)
>>> anova(modelH2a_febmarch)
>>> coefficients(modelH2_febmarch) # get coefficients to make beta vectors
>>> for test
>> (Intercept) X3feb_march
>>5.3421301.172821
>>> coefficients(modelH2a_febmarch)
>> (Intercept) X3feb_march X4feb_march
>> 5.2936263 1.0353752 0.2407557
>> # Test Stat
>>> bsres_fm<-matrix(c(5.342130,1.172821,0),nrow=3)
>>> bsunres_fm<-matrix(c(5.2936263,1.0353752,0.2407557),nrow=3)
>>> bsfm<-t(bsunres_fm-bsres_fm)
>>> #X matrix
>>> fmxs<-cbind(fmvect,X3feb_march,X4feb_march)
>>> tfmxs<-t(fmxs)
>>> xcoms_fm<-(tfmxs %*% fmxs)
>>> xinvs_fm<-ginv(xcoms_fm)
>>> var_fms<-xinvs_fm*0.341
>>> chi_fms<-bsfm %*% var_fms %*% (bsunres_fm-bsres_fm)
>>> chi_fms
>> # Record Chisq value
>>
>> Does this help?
>> Here lffeb_march is the combination of Feb and March log flows
>> and llfeb_march is the combination of Feb and March log loads
>> X3: lffeb_march-mean(feb_march)
>> X4: X1*X3
>>
>> Thanks
>>
>> Rui Barradas wrote
>>>
>>> Hello,
>>>
>>> I'm not at all sure if I understand your problem. Does this describe it?
>>>
>>>
&g
Re: [R] Automating R for Hypothesis Testing
Rui-
Just a quick question. I understand your comment on using ANOVA, but
doesn't this only test for similarities of the mean. We are trying to see if
a the same model can fit for two or three months, therefore have the similar
slope and intercept. The ANOVA would only do one part of this correct, with
the F-Test?
Thanks!
Meredith
Rui Barradas wrote
>
> Hello,
>
> I'm glad it helped.
>
> As for your second question, I don't know, but I'm not very comfortable
> with the way you're doing things.
> Why subtract the coefficients of model 1 from model 2?
> And why the dummy? Why set model 1 to zero?
>
> Isn't it better to use anova's F? After all, it's designed for it, for the
> linear model...
> And if you really want/need the dummy, wouldn't a nested anova do it? (F
> statistic, once again.)
>
> anova(model1, model2)
>
> is simple and statistically speaking seems to me much better. (I specially
> don't like the subtraction bit.)
>
> Rui Barradas
>
> meredith wrote
>>
>> Rui-
>> Thanks this definitely helps, just one quick question. How would you
>> code the values of chi-fm and chi-fms to change based on the degrees of
>> freedom of each model H(i)?
>>
>> Meredith
>>
>>
>> Rui Barradas wrote
>>>
>>> Hello,
>>>
>>> Yes, it does help. Now we can see your data and what you're doing.
>>> What follows is a suggestion on what you could do, not full solution.
>>> (You forgot to say what X1 is, but I don't think it's important to
>>> understand the suggestion.)
>>> (If I'm wrong, say something.)
>>>
>>>
>>> milwaukeephos <- read.csv("milwaukeephos.csv", header=TRUE,
>>> stringsAsFactors=FALSE)
>>> # list of data.frames, one per month
>>> ls1 <- split(milwaukeephos, milwaukeephos$month)
>>>
>>> #- if you want to keep the models, not needed if you don't.
>>> # (yoy probably don't)
>>> modelH <- vector("list", 12)
>>> modelHa <- vector("list", 12)
>>> modelH2 <- vector("list", 12)
>>> modelH2a <- vector("list", 12)
>>> #- values to record, these are needed, create them beforehand.
>>> chi_fm <- numeric(12)
>>> chi_fms <- numeric(12)
>>> #
>>> seq_months <- c(1:12, 1) # wrap months around.
>>> for(i in 1:12){
>>> month_this <- seq_months[i]
>>> month_next <- seq_months[i + 1]
>>>
>>> lload <- c(ls1[[month_this]]$load_kg, ls1[[month_next]]$load_kg)
>>> lflow <- c(ls1[[month_this]]$flow, ls1[[month_next]]$flow)
>>> modelH[[i]] <- lm(lload ~ lflow)
>>> # If you don't want to keep the models, use modelH only
>>> # ( without [[i]] )
>>> # and do the same with X1
>>>
>>> # rest of your code for first test goes here
>>> chi_fm[i] <- bfm %*% var_fm %*% (bunres_fm - bres_fm)
>>>
>>> # and the same for the second test
>>> chi_fms[i] <- ...etc...
>>> }
>>>
>>>
>>> Hope this helps,
>>>
>>> Rui Barradas
>>>
>>>
>>> meredith wrote
>>>>
>>>> dput: http://r.789695.n4.nabble.com/file/n4620188/milwaukeephos.csv
>>>> milwaukeephos.csv
>>>>
>>>> # Feb-march
>>>>> modelH_febmarch<-lm(llfeb_march~lffeb_march)
>>>>>modelHa_febmarch<-lm(llfeb_march~X1feb_mar+lffeb_march)
>>>>> anova(modelHa_febmarch)
>>>>> coefficients(modelH_febmarch)
>>>> (Intercept) lffeb_march
>>>> -2.4298901.172821
>>>>> coefficients(modelHa_febmarch)
>>>> (Intercept) X1feb_mar lffeb_march
>>>> -2.8957776 -0.5272793 1.3016303
>>>>> bres_fm<-matrix(c(-2.429890,0,1.172821),nrow=3)
>>>>> bunres_fm<-matrix(c(-2.8957776,-0.5272793,1.3016303),nrow=3)
>>>>>bfm<-t(bunres_fm-bres_fm)
>>>>> fmvect<-seq(1,1,length=34)
>>>>> X1a_febmar<-seq(0,0,length=9) # dummy variable step 1
>>>>> X1b_febmar<-seq(1,1,length=25) # dummy variable step 2
>>>>> X1feb_mar<-c(X1a_febmar,X1b_febmar) #dummy variable creation
>>>> # Test Stat Equation for Chisq
>>>>> fmxx<-cbind(fmvect,X1feb_mar,lffeb_march)
>>>>> tfmx<-t(fmxx)
>>>>> xcom_fm<-(tfmx %*% f
[R] Event-triggered change in value with a time-delay
Hello,
I am solving a set of ODEs using deSolve and have run into a problem I
would appreciate some advice on. One of the parameters (m) in the ODEs
changes between two states when one of the variables (D) crosses a
threshold (D_T) for the first time in either direction. Additionally, when
the variable crosses the threshold (either by increasing or decreasing),
there is a time delay (delay) before the parameter shifts values.
Toy model:
# define parameters
P <- c(m = 10,
R_1 = 0.5,
R_2 = 0.1,
D_0 = 50,
D_T = 5,
delay = 3
)
# set initial values
y_0 <- c(D = 0, Y = 1)
# set timesteps
times <- seq(0, 100, length = 101)
# define function
bit <- function(t,y,parms){
with(as.list(c(parms,y)),{
# How I would implement the change if there was no time delay
# m <- ifelse(D>=D_T, m*.68, m)
# Option 1: nested if statements with time delay
# if (D >= D_T) {
# t_start1 <- t + delay
# if (t >= t_start1){
# m <- 6.8
# }
# }
# if (D < D_T) {
# t_start2<- t + delay
# if (t >= t_start2){
# m <- 10
# }
# }
dD <- D_0 * (R_1 / (R_1 - R_2)) * (-R_2 * exp(-R_2 * t) + R_1 *
exp(-R_1 * t))
dY <- Y + m;
res <- c(dD, dY)
list(res, m = m)
})
}
# Solve model
library(deSolve)
out <- ode(y = y_0, times = times, func = bit, parms = P)
out.df <- as.data.frame(out)
~~
I have come up with a number of ideas, but I either haven't been able to
get them to work or don't know if R has a command for it:
1. Nested "if" statements with something to track time (as shown above).
One issue here is* t-start *is redefined with each time step in the ode
solver so *t* never reaches *t_start +delay*.
2. Event triggered changes. From the examples I have found (link
p.
24-25), the time at which the event happens is either pre-set or occurs at
a root of the ODE. Also- would this strategy require multiple events? One
event in which the variable crossing the threshold triggers a timer,
another for the timer triggering the parameter change?
3. Lagvalue. This might take care of the time lag component, but I am still
struggling to understand (1) how lagvalue works and (2) how to integrate
both the threshold and timer components.
if ( D < D_T )
m<- ifelse(lagvalue(t-delay)<=0, 6.8, 10)
end
4. Sigmoidal function. Someone suggested this to me earlier, but I am not
sure how this would work. Is it possible to have a sigmoidal function that
switches repeatedly between two states (as opposed to the characteristic S
curve that transitions once from one state to another)?
5. Mathmematica has a WhenEvent function- is there something similar in R?
Thank you for your time and any tips you might have!
[[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] Event-triggered change in value with a time-delay
Hi David,
Thank you for your suggestion- I will look into it some more.
Thanks again,
Hannah
On Wed, Apr 18, 2018 at 5:09 PM, David Winsemius
wrote:
>
> > On Apr 18, 2018, at 1:04 AM, Hannah Meredith
> wrote:
> >
> > Hello,
> >
> > I am solving a set of ODEs using deSolve and have run into a problem I
> > would appreciate some advice on. One of the parameters (m) in the ODEs
> > changes between two states when one of the variables (D) crosses a
> > threshold (D_T) for the first time in either direction. Additionally,
> when
> > the variable crosses the threshold (either by increasing or decreasing),
> > there is a time delay (delay) before the parameter shifts values.
>
> I think you migh not have fully reviewed the available functions in that
> package. Doing an sos::findFn-search on ODE's with time delay brings up
> help pages for the `dede` function in deSolve.
>
> For the conditional value of the within regime-fixed parameter `m`,
> couldn't you just use this:
>
> ... + m*(6.8+ 3.2*(D < D_T))
>
> Caveat: very lightly tested in your toy example, which did generate a
> result where D cross that threshold without blowing up.
>
> --
>
> David.
>
>
> >
> > Toy model:
> > # define parameters
> > P <- c(m = 10,
> > R_1 = 0.5,
> > R_2 = 0.1,
> > D_0 = 50,
> > D_T = 5,
> > delay = 3
> > )
> > # set initial values
> > y_0 <- c(D = 0, Y = 1)
> >
> > # set timesteps
> > times <- seq(0, 100, length = 101)
> >
> > # define function
> > bit <- function(t,y,parms){
> > with(as.list(c(parms,y)),{
> >
> >
> ># How I would implement the change if there was no time delay
> ># m <- ifelse(D>=D_T, m*.68, m)
> >
> ># Option 1: nested if statements with time delay
> ># if (D >= D_T) {
> ># t_start1 <- t + delay
> ># if (t >= t_start1){
> ># m <- 6.8
> ># }
> ># }
> ># if (D < D_T) {
> ># t_start2<- t + delay
> ># if (t >= t_start2){
> ># m <- 10
> ># }
> ># }
> >
> >
> >dD <- D_0 * (R_1 / (R_1 - R_2)) * (-R_2 * exp(-R_2 * t) + R_1 *
> > exp(-R_1 * t))
> >dY <- Y + m;
> >res <- c(dD, dY)
> >list(res, m = m)
> > })
> > }
> >
> >
> > # Solve model
> > library(deSolve)
> > out <- ode(y = y_0, times = times, func = bit, parms = P)
> > out.df <- as.data.frame(out)
> > ~~
> > I have come up with a number of ideas, but I either haven't been able to
> > get them to work or don't know if R has a command for it:
> >
> > 1. Nested "if" statements with something to track time (as shown above).
> > One issue here is* t-start *is redefined with each time step in the ode
> > solver so *t* never reaches *t_start +delay*.
> >
> > 2. Event triggered changes. From the examples I have found (link
> > <https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf>
> p.
> > 24-25), the time at which the event happens is either pre-set or occurs
> at
> > a root of the ODE. Also- would this strategy require multiple events?
> One
> > event in which the variable crossing the threshold triggers a timer,
> > another for the timer triggering the parameter change?
> >
> > 3. Lagvalue. This might take care of the time lag component, but I am
> still
> > struggling to understand (1) how lagvalue works and (2) how to integrate
> > both the threshold and timer components.
> >
> > if ( D < D_T )
> > m<- ifelse(lagvalue(t-delay)<=0, 6.8, 10)
> > end
> >
> > 4. Sigmoidal function. Someone suggested this to me earlier, but I am not
> > sure how this would work. Is it possible to have a sigmoidal function
> that
> > switches repeatedly between two states (as opposed to the characteristic
> S
> > curve that transitions once from one state to another)?
> >
> > 5. Mathmematica has a WhenEvent function- is there something similar in
> R?
> >
> >
> > Thank you for your time and any tips you might have!
> >
> > [[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.
>
> David Winsemius
> Alameda, CA, USA
>
> 'Any technology distinguishable from magic is insufficiently advanced.'
> -Gehm's Corollary to Clarke's Third Law
>
>
>
>
>
>
--
Hannah Meredith | PhD
Whitaker Scholar
London School of Hygiene and Tropical Medicine
[[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] negative binomial GAMM with variance structures
Hello, I am having some difficulty converting my gam code to a correct gamm code, and I'm really hoping someone will be able to help me. I was previously using this script for my overdispersed gam data: M30 <-gam(efuscus~s(mic, k=7) +temp +s(date)+s(For3k, k=7) + pressure+ humidity, family=negbin(c(1,10)), data=efuscus) My gam.check gave me the attached result. In order to deal with my heterogeneity, I need to switch over to a gamm structure and use at least one, but possibly multiple, variance structures, and I am starting by applying varPower to my temperature covariate. (Efuscus is my square root transformed response variable). Here is the code I have for the gamm: K1 <-(efuscus~s(mic, k=7) +temp +s(date)+s(For3k, k=7) + pressure+ humidity+ windspeed + year) M17.4A <-gamm(K1, method="REML", family=negbin(c(1,10), data=efuscus, weights=varPower(form=~temp)) With my various versions of the script, the two error messages I have been getting repeatedly are: Error in eval(expr, envir, enclos) : object 'temp' not found, and Error: unexpected symbol in: "M17.4A <-gamm(K1, method="REML", family=negbin(c(1,10), data=efuscus, weights=varPower(form=~temp))" I know I must be missing something obvious that is causing my script not to work...I've been looking at both Wood (2006) and Zuur's Mixed Models book for examples but none of them are completely the same as my situation, which is causing me to get tripped up. Am I defining my negative binomial family correctly for the situation? Do I need to somehow define the variance term before I apply it? Thank you very much. I appreciate your time and patience. Meredith __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/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] normalizing a negative binomial distribution and/or incorporating variance structures in a GAMM
Hello everyone, Apologies in advance, as this is partially a stats question and partially an R question. I have been using a GAM to model the activity level of bats going into and coming out from a forested edge. I had eight microphones set up in a line transect at each of eight sites, and I am hoping to construct a model for each of 7 species. My count data has a reverse J-shaped skew and is overdispersed with a fair amount of zeros, and I haven't found any transformations that will completely normalize it (I've tried square roots and logs). Meanwhile, the variance in call numbers varies between sites and between microphones. I wanted to use a GAMM to incorporate varComb and varIdent, but these can only be applied on data with a gaussian distribution. Are there any packages I should be looking into that I don't know about that will apply a variance structure on a negative binomial distribution? Or is there some transformation that I should be using that will solve my normality issues? I've been searching the R-help boards, everything in Zuur and Woods, but I haven't found an answer yet. Thanks very much, I am very appreciative of any help I can get on this matter. Sincerely, Meredith Jantzen M.Sc. Candidate Department of Biology University of Western Ontario [[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] problem with glsm.krige: trendd and trend l must have similar specifications error
Hello, I used glsm.mcmc and likfit.glsm to create model. Now I want to predict at different locations, but I can't get glsm.krige to work. I keep getting the error that trend.d and trend.l must have similar specifications. I have tried numerous ways to include the covariates in the glsm.krige model, and I keep getting the same error message. The bolded part is the part that doesn't work. Any help would be appreciated! reddmcmc=glsm.mcmc(geodata, coords=geodata$coords, data=geodata$data,model=list(family="poisson", cov.pars=c(50,8000),beta=c(-26,-.04,6.94),cov.model="spherical",nugget=5,trend=~powerb1+temp1),mcmc.input = mcmc.control(S.scale =.01 ,thin=10)) prepareredd1=prepare.likfit.glsm(reddmcmc) lik.redd1.1 <- likfit.glsm(prepareredd1, ini.phi = 10,cov.model = "spherical",trend=~powerb1 + temp1) # predict at new locations trendpred=trend.spatial(~habitat2$powerb + habitat2$avgtemp) test2 <- glsm.krige(reddmcmc,coords,trend.l=trendpred) # trend.d and trend.l must have similar specifications ? test2 <- glsm.krige(reddmcmc,coords,trend.l=habitat2$powerb + habitat2$avgtemp) # same error message... [[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] bug in R2WinBUGS
Problem with the change is that none of our old scripts work! If the model.file is in R's current working directory we either have to use working.directory=getwd() or specify a full path with model.file=file.path(getwd(), "mymodel.bug"). It would be really nice if 'bugs' looked in R's c.w.d. for input files -- in particular model.file -- while still putting its temporary files in a tempdir. Cheers, Mike. Uwe Ligges-3 wrote: > > > > John Smith wrote: >> In newest version of R2WinBUGS, the default directory is changed to >> working.directory, but never changed back once finished bugs call. > > Now finally fixed, thanks again for the report and sorry for the delay. > > Best, > Uwe > > > > >> [[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. > > -- View this message in context: http://www.nabble.com/bug-in-R2WinBUGS-tp21377525p21761975.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 find "p" in binomial(n,p)
I think the function you need is 'help.search'; try:
help.search("binomial")
and look for something obvious in the 'stats' package. A good deal quicker
and easier than posting to an internet forum!
Cheers, Mike.
cathelf wrote:
>
> Dear all,
>
> I am trying a find the value "p" in binomial.
>
> X ~ Bin(n,p)
>
> I want to find the value "p", so that Pr(X <= k) <= alpha
>
> Here, n, k and alpha are known. n, k are integers. alpha is between (0,1).
>
> Thanks a lot!
>
> Catherine
>
--
View this message in context:
http://www.nabble.com/how-to-find-%22p%22-in-binomial%28n%2Cp%29-tf4484227.html#a12787900
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] Comparing density plots using same axes or same axes scale
Good Evening-
I have a set of nine scenarios I want to plot to see how the distribution
is changing, if one tail is getting larger in certain scenario, currently I
am using this code
colnames<-dimnames(sag_pdfs)[[2]]
par(mfrow=c(3,3))
for(i in 1:9) {
d<-density(sag[,i])
plot(d,type="n", main=colnames[i])
polygon(d,col="red",border="grey")}
where sag is a 7305x9 double matrix and 9 different scenarios. I want to be
able to compare the distribution using the same axes scale.
Can anyone help?
Thanks
Meredith LaBeau
--
Doctoral Candidate
Department of Civil and Environmental Engineering
Michigan Technological University
[[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] Merging multiple columns into one column
Good Evening- I have a dataframe that has 10 columns that has a header and 7306 rows in each column, I want to combine these columns into one. I utilized the stack function but it only returned 3/4 of the data...my code is: where nfcuy_bw is the dataframe with 7305 obs. and 10 variables Once I apply this code I only receive a data frame with 58440 obs. of 2 variables, of which there should be 73,050 obs. of 2 variables, just wondering what is happening here? View(nfcuy_bw) attach(nfcuy_bw) cuyahoga_nf<-data.frame(s5,s10,s25,s27,s33,s41,s51,his_c) cuy_nf<-stack(cuyahoga_nf) Thanks Meredith -- Doctoral Candidate Department of Civil and Environmental Engineering Michigan Technological University [[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] Removing lower whisker in boxplot to see the effects of the high values
Good Afternoon- I was wanting to alter the boxplot to remove the lower whisker, both the whisker line and staple just on the lower end. Is there a way to do this? As my code is currently: boxplot(log_loads~ind,data=nfmaum, horizontal=TRUE, notch=T, outline=FALSE, whisker=0, main="Maumee River Near Future Climate Scenarios", ylab="Log Load",xlab="Climate Scenarios") I just want to better see the medians and high end tail. Thanks Meredith -- Doctoral Candidate Department of Civil and Environmental Engineering Michigan Technological University [[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] using correlation compound correlation structure with nlme; how to incorporate multple random effects?
Hello, I have developed this model to test change in PTFines6 over time. I have random effects of watershed (HUC3) and management type (mgmt3), and then I have the YrC/SiteID random effect which is the longitudinal time effect. But I recently found out that I need to incorporate a compound correlation structure. I found examples of this where there is just one random effect, but not my case where I have multiple. I am not sure if I just include YrC|SiteID or the entire random effect. I have tried both, and get error messages: Original code: model5.15_e=lme(PTFines6~HUC3 + YrC*mgmt3 + Bf* mgmt3 + LnGrad * mgmt3+ Precip ,random=list( ~1|HUC3,~1|mgmt3,1~(YrC)|SiteID),na.action=na.omit, data=habitat2,method="REML",control=control1) #22835 no precip interaction I have tried: model5.15_e=lme(PTFines6~HUC3 + YrC*mgmt3 + Bf* mgmt3 + LnGrad * mgmt3+ Precip ,random=list( ~1|HUC3,~1|mgmt3,1~(YrC)|SiteID),na.action=na.omit, data=habitat2, correlation=corCompSymm(form=1~(YrC)|SiteID),method="REML",control=control1) #22835 no precip interaction I get the message: incompatible formulas for groups in "random" and "correlation" I have tried: model5.15_e=lme(PTFines6~HUC3 + YrC*mgmt3 + Bf* mgmt3 + LnGrad * mgmt3+ Precip ,random=list( ~1|HUC3,~1|mgmt3,1~(YrC)|SiteID),na.action=na.omit, data=habitat2, correlation=corCompSymm(form=~1|HUC3,~1|mgmt3,1~(YrC)|SiteID),method="REML",control=control1) #22835 no precip interaction I get the error message "non-numeric argument to mathematical function Thanks for any help. Christy Meredith This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately. [[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] randomly select another observation with same grouping factor and year value, do for every record in dataframe
Hello,
I am trying to create a function that will move through each record of a data
frame, find the value in the "HUC" column, then
randomly select another observation from the dataframe with the same value in
"HUC" column, as well as the same value in "Yr" column as the first
observation. I want the function to produce a list of the RchID of the first
observation, the RchID of the second randomly chosen observation, and several
other characteristics of the randomly chosen observation. Below is the code I
have written, but it doesn't work.
Thanks for any help.
Christy
test
roads=read.csv("streamland23.csv")
for (i in 1:nrow (roads)){
Sitetype= roads$Sitetype
yr=roads$REACH_Yr
initRchid=roads$RchID
huc1=roads$HUC
sample.df <- function(df, n) df[sample(nrow(df), n), , drop = FALSE]
selected=sample.df(roads[roads$HUC == "huc1"& roads$REACH_Yr =="yr" , ], 1)
output=cbind (initRchid,selected$RchID,selected$Sitetype,selected$REACH_Yr)
}
Christy Meredith
USDA Forest Service
Rocky Mountain Research Station
PIBO Monitoring
Data Analyst
Voice: 435-755-3573
Fax: 435-755-3563
This electronic message contains information generated by the USDA solely for
the intended recipients. Any unauthorized interception of this message or the
use or disclosure of the information it contains may violate the law and
subject the violator to civil or criminal penalties. If you believe you have
received this message in error, please notify the sender and delete the email
immediately.
[[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 with for loop: new column giving count of observation for each SITEID
Hello,
I think this is easy, but I can't seem to find a good way to do this in the R
help. I have a list of sites, with multiple years of data for each site id. I
want to create a new column that gives a number describing whether it is the
1st year ("1" ) the data was collected for the site, the second year ("2"),
etc. I have different years for each siteid, but I don't care which year it was
collected, just the order that it is in for that siteid. This is what I have
so far, but it doesn't do the analysis separately for each SiteID.
indexi<-indexg[order(indexg$SiteID,indexg$Yr),]
obs=0
indexi=na.omit(indexi)
for(i in 1:length(indexi$SiteID)){
obs=obs+1
indexi$obs[i]=obs
}
Thanks for any help you can give.
Christy Meredith
USDA Forest Service
Rocky Mountain Research Station
PIBO Monitoring
Data Analyst
Voice: 435-755-3573
Fax: 435-755-3563
This electronic message contains information generated by the USDA solely for
the intended recipients. Any unauthorized interception of this message or the
use or disclosure of the information it contains may violate the law and
subject the violator to civil or criminal penalties. If you believe you have
received this message in error, please notify the sender and delete the email
immediately.
[[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] help with for loop: new column giving count of observation for each SITEID
Not quite,
I need it like this, a new number for each ordered year in the sequence within
each site, regardless of what the years are, and to retain the RchID column.
RchID siteyearindex
1 A 20021
2 A 20042
3 A 20053
4 B 20031
5 B 20062
6 B 20083
7 C 20021
8 C 20032
9 C 20043
Thanks so much for you help!
-Original Message-
From: William Dunlap [mailto:wdun...@tibco.com]
Sent: Tuesday, October 30, 2012 1:07 PM
To: Meredith, Christy S -FS; r-help@R-project.org
Subject: RE: [R] help with for loop: new column giving count of observation for
each SITEID
Is this what you want?
> withinGroupIndex <- function(group, ...) ave(integer(length(group)), group,
..., FUN=seq_along)
> site <- c("A","A","C","D","C","A","B")
> data.frame(site, index=withinGroupIndex(site))
site index
1A 1
2A 2
3C 1
4D 1
5C 2
6A 3
7B 1
You can add more arguments if the groups depend on more than one value:
> year <- rep(c(1985, 2012), c(4,3))
> data.frame(site, year, index=withinGroupIndex(site, year))
site year index
1A 1985 1
2A 1985 2
3C 1985 1
4D 1985 1
5C 2012 1
6A 2012 1
7B 2012 1
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -Original Message-
> From: r-help-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org] On Behalf Of Meredith, Christy S
> -FS
> Sent: Tuesday, October 30, 2012 11:17 AM
> To: r-help@R-project.org
> Subject: [R] help with for loop: new column giving count of
> observation for each SITEID
>
>
> Hello,
> I think this is easy, but I can't seem to find a good way to do this
> in the R help. I have a list of sites, with multiple years of data for
> each site id. I want to create a new column that gives a number
> describing whether it is the 1st year ("1" ) the data was collected
> for the site, the second year ("2"), etc. I have different years for
> each siteid, but I don't care which year it was collected, just the order
> that it is in for that siteid. This is what I have so far, but it doesn't do
> the analysis separately for each SiteID.
>
> indexi<-indexg[order(indexg$SiteID,indexg$Yr),]
>
> obs=0
> indexi=na.omit(indexi)
> for(i in 1:length(indexi$SiteID)){
> obs=obs+1
> indexi$obs[i]=obs
> }
>
>
> Thanks for any help you can give.
>
> Christy Meredith
> USDA Forest Service
> Rocky Mountain Research Station
> PIBO Monitoring
> Data Analyst
> Voice: 435-755-3573
> Fax: 435-755-3563
>
>
>
>
>
> This electronic message contains information generated by the USDA
> solely for the intended recipients. Any unauthorized interception of
> this message or the use or disclosure of the information it contains
> may violate the law and subject the violator to civil or criminal
> penalties. If you believe you have received this message in error, please
> notify the sender and delete the email immediately.
>
> [[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] help with for loop: new column giving count of observation for each SITEID
Thanks to you all for your help with this code Basically, the purpose of this
was to create a column detailing whether it is the 1st, 2nd, 3rd, etc time that
the site was sampled. I need this for creating a graph which shows this new
variable on the x axis versus a variable of interest. The variable of interest
can change but refers to aspects of fish habitat.
Indexg$newindex refers to this new variable (perhaps not the best names, but
they are indices of habitat quality and time). Ultimately, I want to create a
ggplot graph with the new value on the x axis and the variable of interest in
the y axis (code also shown below). I have been able to make the graph, but not
all of my series are showing on the graph (only the first 5). I am not sure why
this is the case. So now I have another question. Do any of you know why all
the series are not showing on my graph and how to fix this?
Thanks
Christy
#Code for new column nth time sampled
indexg=read.csv("indexg.csv")
indexg=data.frame(indexg)
indexi<-indexg[order(indexg$SiteID,indexg$Yr),]
res<-do.call(rbind,lapply(split(indexi,indexi$SiteID),function(x)
data.frame(x,newindex=1:nrow(x
rownames(res)<-1:nrow(res)
res
res$SiteID=as.factor(res$SiteID)
q=ggplot(data=res,aes(x=newindex,y=Bf,group=SiteID,shape=SiteID))+geom_line()
+ geom_point(size=5,colour="black")
Data:
SiteID Yr Bf newindex
11015 2001 3.771
21015 2006 4.942
31015 2011 5.203
41035 2003 5.841
51035 2008 6.182
61039 2003 4.411
71039 2008 5.242
81047 2001 NA1
91047 2003 4.762
10 1047 2004 4.103
11 1047 2006 5.854
12 1047 2008 4.875
13 1047 2009 5.596
14 1047 2010 4.697
15 1047 2011 4.948
16 1088 2003 3.751
17 1088 2008 5.862
18 1104 2004 7.321
19 1104 2009 7.212
20 1106 2001 6.921
21 1106 2002 4.252
22 1106 2003 4.753
23 1106 2004 6.674
24 1106 2005 4.505
25 1106 2006 6.626
26 1106 2008 6.327
27 1106 2009 6.308
28 1106 2010 6.659
29 1106 2011 5.51 10
30 1110 2004 6.871
31 1110 2009 5.532
32 2702 2009 1.801
33 2944 2010 4.361
34 2946 2010 2.251
-Original Message-
From: arun [mailto:smartpink...@yahoo.com]
Sent: Tuesday, October 30, 2012 1:57 PM
To: Meredith, Christy S -FS
Cc: R help; William Dunlap
Subject: Re: [R] help with for loop: new column giving count of observation for
each SITEID
HI,
You can also use this:res<-do.call(rbind,lapply(split(d,d$site),function(x)
data.frame(x,newindex=1:nrow(x
rownames(res)<-1:nrow(res)
res
# RchID site year index newindex
#1 1 A 2002 1 1
#2 2 A 2004 2 2
#3 3 A 2005 3 3
#4 4 B 2003 1 1
#5 5 B 2006 2 2
#6 6 B 2008 3 3
#7 7 C 2002 1 1
#8 8 C 2003 2 2
#9 9 C 2004 3 3
A.K.
- Original Message -
From: William Dunlap
To: "Meredith, Christy S -FS"
Cc: "r-help@r-project.org"
Sent: Tuesday, October 30, 2012 3:43 PM
Subject: Re: [R] help with for loop: new column giving count of observation for
each SITEID
Your data was, in R-readable format (from dput())
d <- data.frame(
RchID = 1:9,
site = factor(c("A", "A", "A", "B", "B", "B", "C",
"C", "C"), levels = c("A", "B", "C")),
year = c(2002L, 2004L, 2005L, 2003L, 2006L, 2008L,
2002L, 2003L, 2004L),
index = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L)) and I am assuming that
'index' is the desired result. You can use withinGroupIndex to make a new
column identical to 'index'. There are a variety of ways to add that column to
an existing data.frame, one of which is within():
> within(d, newIndex <- withinGroupIndex(site))
RchID site year index newIndex
1 1 A 2002 1 1
2 2 A 2004 2 2
3 3 A 2005 3 3
4 4 B 2003 1 1
5 5 B 2006 2 2
6 6 B 2008 3 3
7 7 C 2002 1 1
8 8 C 2003 2 2
9 9 C 2004 3 3
Or is 'index' not the desired result?
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -Original Message-
> From: Meredith, Christy S -FS [mailto:csmered...@fs.fed.us]
> Sent: Tuesday, October 30, 2012 12:20 PM
> To: William Dunlap
> Subject: RE: [R] help with for loop: new column giving count of
> observation for each SITEID
>
> Not quite,
> I need it like th
