Hey I'm not really sure what I should put on here, but I am having trouble
with my R code. I am trying to get the p-values, R^2s etc for a number of
different groups of variables that are all in one dataset.
This is the code:
#Stand counter
st<-1
#Collections
stands<-numeric(67)
slopes<-numeric(67)
intercepts<-numeric(67)
mses<-numeric(67)
rsquares<-numeric(67)
pValues<-numeric(67)
#Start lists for X and Y values within each stand
xi<-numeric(0)
yi<-numeric(0)
#Set the first element to the starting X and Y values
xi[1]=X[1]
yi[1]=Y[1]
#Start looping working through your data, record by record
for (i in 2:length(X)) {
#If you are in the same stand as on the last record, continue to
#collect X and Y values
if(Stand[i]==Stand[i-1]) {
xi=cbind(xi,X[i])
yi=cbind(yi,Y[i])
} else {
#If a new stand is encountered make your linear model and
#collect statistics
model<-lm(yi~xi)
stands[st]<-Stand[i-1]
intercepts[st]<-model$coefficients[1]
slopes[st]<-model$coefficients[2]
mses[st]<-sum(resid(model)^2)/(length(xi)-2)
ssr<-var(yi)*(length(xi)-1)-sum(resid(model)^2)
rsquares[st]<-ssr/(var(yi)*(length(xi)-1))
fRatio<-ssr/mses[st]
pValues[st]<-1-pf(fRatio,1,length(xi)-2)
#Increment the stand number, zero the within stand collections,
#and start again
st<-st+1
xi<-numeric(0)
yi<-numeric(0)
xi[1]=X[i]
yi[1]=Y[i]
}
}
#Make your data set
standEstimates<-data.frame(standID=stands,intercept=intercepts,slop=slopes,mse=mses,rSquare=rsquares,pValue=pValues)
The standEstimate outputs look like this:
standID intercept slop mse rSquare pValue
1 6833 319.2470 NA 0 NA NA
2 756 708.7470 NA 0 NA NA
3 795 508.2290 NA 0 NA NA
4 1249 503.1460 NA 0 NA NA
5 1331 703.0620 NA 0 NA NA
6 1417 747.7620 NA 0 NA NA
7 4715 549.3400 NA 0 NA NA
8 4850 603.9940 NA 0 NA NA
9 2105 573.3270 NA 0 NA NA
Etc etc
and I get these warnings:
1: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
2: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
number of items to replace is not a multiple of replacement length
3: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
4: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
number of items to replace is not a multiple of replacement length
5: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
6: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
number of items to replace is not a multiple of replacement length
7: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
8: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
number of items to replace is not a multiple of replacement length
9: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
10: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
number of items to replace is not a multiple of replacement length
11: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
12: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
number of items to replace is not a multiple of replacement length
13: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
number of items to replace is not a multiple of replacement length
14: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
etc etc
Sorry if I haven't set this post out right, or haven't provided enough
information. But can anyone see why it is not giving me any returns for R^2
etc?
--
View this message in context:
http://r.789695.n4.nabble.com/R-looping-help-tp4667180.html
Sent from the R help mailing list archive at Nabble.com.
You can use basic debugging techniques to help understand what is going
wrong. For example, insert a cat() statement right before the statement
which is giving the warning. Like this:
## your code
mses[st]<-sum(resid(model)^2)/(length(xi)-2)
ssr<-var(yi)*(length(xi)-1)-sum(resid(model)^2)
## insert two lines
tmp <- (var(yi)*(length(xi)-1))
cat('i',i,'st',st,'ssr',ssr,'denominator',tmp,'\n')
## continue with your code
rsquares[st]<-ssr/(var(yi)*(length(xi)-1))
fRatio<-ssr/mses[st]
This example should help you understand the warning message
> foo <- numeric(5)
> foo[2] <- c(1,3)
Warning message:
In foo[2] <- c(1, 3) :
number of items to replace is not a multiple of replacement length
--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
On 5/16/13 3:17 AM, "rishard" <orc15 at uclive.ac.nz> wrote:
>Hey I'm not really sure what I should put on here, but I am having
trouble
>with my R code. I am trying to get the p-values, R^2s etc for a number of
>different groups of variables that are all in one dataset.
>
>This is the code:
>
>#Stand counter
>st<-1
>#Collections
>stands<-numeric(67)
>slopes<-numeric(67)
>intercepts<-numeric(67)
>mses<-numeric(67)
>rsquares<-numeric(67)
>pValues<-numeric(67)
>#Start lists for X and Y values within each stand
>xi<-numeric(0)
>yi<-numeric(0)
>#Set the first element to the starting X and Y values
>xi[1]=X[1]
>yi[1]=Y[1]
>#Start looping working through your data, record by record
>for (i in 2:length(X)) {
> #If you are in the same stand as on the last record, continue to
> #collect X and Y values
> if(Stand[i]==Stand[i-1]) {
> xi=cbind(xi,X[i])
> yi=cbind(yi,Y[i])
> } else {
> #If a new stand is encountered make your linear model and
> #collect statistics
> model<-lm(yi~xi)
> stands[st]<-Stand[i-1]
> intercepts[st]<-model$coefficients[1]
> slopes[st]<-model$coefficients[2]
> mses[st]<-sum(resid(model)^2)/(length(xi)-2)
> ssr<-var(yi)*(length(xi)-1)-sum(resid(model)^2)
> rsquares[st]<-ssr/(var(yi)*(length(xi)-1))
> fRatio<-ssr/mses[st]
> pValues[st]<-1-pf(fRatio,1,length(xi)-2)
> #Increment the stand number, zero the within stand collections,
> #and start again
> st<-st+1
> xi<-numeric(0)
> yi<-numeric(0)
> xi[1]=X[i]
> yi[1]=Y[i]
> }
>}
>#Make your data set
>standEstimates<-data.frame(standID=stands,intercept=intercepts,slop=slopes
>,mse=mses,rSquare=rsquares,pValue=pValues)
>
>The standEstimate outputs look like this:
>
>standID intercept slop mse rSquare pValue
>1 6833 319.2470 NA 0 NA NA
>2 756 708.7470 NA 0 NA NA
>3 795 508.2290 NA 0 NA NA
>4 1249 503.1460 NA 0 NA NA
>5 1331 703.0620 NA 0 NA NA
>6 1417 747.7620 NA 0 NA NA
>7 4715 549.3400 NA 0 NA NA
>8 4850 603.9940 NA 0 NA NA
>9 2105 573.3270 NA 0 NA NA
>Etc etc
>
>and I get these warnings:
>
>1: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>2: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
>3: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>4: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
>5: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>6: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
>7: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>8: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
>9: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>10: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
>11: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>12: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
>13: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
>14: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>etc etc
>
>Sorry if I haven't set this post out right, or haven't provided
enough
>information. But can anyone see why it is not giving me any returns for
>R^2
>etc?
>
>
>
>--
>View this message in context:
>http://r.789695.n4.nabble.com/R-looping-help-tp4667180.html
>Sent from the R help mailing list archive at Nabble.com.
>
>______________________________________________
>R-help at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.
It's difficult to see where your error is when you provide no accompanying
data for us to test your code out on (X? Y? Stand?).
However, it looks like you are making your code more complex than it needs
to be. Are you simply trying to a fit a separate linear regression to each
subset of your data identified by a unique standID? If so, you can do so
like this:
# a fake version of your data
mydata <- data.frame(standID=rep(1:5, rep(20, 5)), X=rnorm(100),
Y=rnorm(100))
# split the data up by stand ID
standlist <- split(mydata, mydata$standID)
# fit a regression to each stand and save some results
standEst1 <- lapply(standlist, function(df) {
fit <- summary(lm(Y ~ X, data=df))
co <- coef(fit)
intercept <- co[1, 1]
slope <- co[2, 1]
pValue <- co[2, 4]
mse <- fit$sigma^2
rSquare <- fit$r.squared
c(standID=df$standID[1], intercept=intercept, slope=slope, mse=mse,
rSquare=rSquare, pValue=pValue)
})
# convert the results to a data frame
standEstimates <- data.frame(do.call(rbind, standEst1))
Jean
On Thu, May 16, 2013 at 5:17 AM, rishard <orc15@uclive.ac.nz> wrote:
> Hey I'm not really sure what I should put on here, but I am having
trouble
> with my R code. I am trying to get the p-values, R^2s etc for a number of
> different groups of variables that are all in one dataset.
>
> This is the code:
>
> #Stand counter
> st<-1
> #Collections
> stands<-numeric(67)
> slopes<-numeric(67)
> intercepts<-numeric(67)
> mses<-numeric(67)
> rsquares<-numeric(67)
> pValues<-numeric(67)
> #Start lists for X and Y values within each stand
> xi<-numeric(0)
> yi<-numeric(0)
> #Set the first element to the starting X and Y values
> xi[1]=X[1]
> yi[1]=Y[1]
> #Start looping working through your data, record by record
> for (i in 2:length(X)) {
> #If you are in the same stand as on the last record, continue to
> #collect X and Y values
> if(Stand[i]==Stand[i-1]) {
> xi=cbind(xi,X[i])
> yi=cbind(yi,Y[i])
> } else {
> #If a new stand is encountered make your linear model and
> #collect statistics
> model<-lm(yi~xi)
> stands[st]<-Stand[i-1]
> intercepts[st]<-model$coefficients[1]
> slopes[st]<-model$coefficients[2]
> mses[st]<-sum(resid(model)^2)/(length(xi)-2)
> ssr<-var(yi)*(length(xi)-1)-sum(resid(model)^2)
> rsquares[st]<-ssr/(var(yi)*(length(xi)-1))
> fRatio<-ssr/mses[st]
> pValues[st]<-1-pf(fRatio,1,length(xi)-2)
> #Increment the stand number, zero the within stand collections,
> #and start again
> st<-st+1
> xi<-numeric(0)
> yi<-numeric(0)
> xi[1]=X[i]
> yi[1]=Y[i]
> }
> }
> #Make your data set
>
>
standEstimates<-data.frame(standID=stands,intercept=intercepts,slop=slopes,mse=mses,rSquare=rsquares,pValue=pValues)
>
> The standEstimate outputs look like this:
>
> standID intercept slop mse rSquare pValue
> 1 6833 319.2470 NA 0 NA NA
> 2 756 708.7470 NA 0 NA NA
> 3 795 508.2290 NA 0 NA NA
> 4 1249 503.1460 NA 0 NA NA
> 5 1331 703.0620 NA 0 NA NA
> 6 1417 747.7620 NA 0 NA NA
> 7 4715 549.3400 NA 0 NA NA
> 8 4850 603.9940 NA 0 NA NA
> 9 2105 573.3270 NA 0 NA NA
> Etc etc
>
> and I get these warnings:
>
> 1: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 2: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
> 3: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 4: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
> 5: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 6: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
> 7: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 8: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
> 9: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 10: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
> 11: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 12: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> number of items to replace is not a multiple of replacement length
> 13: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
> number of items to replace is not a multiple of replacement length
> 14: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> etc etc
>
> Sorry if I haven't set this post out right, or haven't provided
enough
> information. But can anyone see why it is not giving me any returns for
> R^2
> etc?
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/R-looping-help-tp4667180.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.
>
[[alternative HTML version deleted]]