title: Choose canonical values from clusters of erroneous ones

text: I suspect this is a statistical problem, but it may be just an algorithmic one. So, more formally: Given:   * a set C of unknown ‘canonical' values   * an error window e   * a set V of known values, s.t. every v in V has the property that   `c - e <= v <= c + e for some c in C` Find all c in C for which there are associated v values, and find a (non- unique) c for every v. I have a set of values representing spar lengths on a geodesic dome. Each of these values ought to be one of, say, 12 canonical lengths: all Spar A’s would be 5.2 units long, for example. But, in the models I’ve been given, some of the spars were hand-drawn and there may be small errors in the lengths. I want a way to automatically identify the categories (i.e. “There’s a cluster around 5.2 units. Let’s call that A”) within some error range and assign each [non-anomalous] value to a category. This is easy to do by hand with a histogram, but I’d like to do it programmatically. This seems like the kind of problem statisticians have tools for, but I’m not sure what they are. So I thought I’d ask. I can go through and brute-force for every v, counting how many other values fall within a window of `+/- e`, then use that count to determine the most likely canonical values, but I suspect there’s some conceptual formalism that would resolve this more cleanly. Any ideas? Are there some stats concepts that would help me answer this kind of question better? Anybody who wants to turn this into a more generalized teaching moment would be welcome.


title: Advantages of bootstrapping when calculating mean

text: I'm a statistically incompetent social scientist who is very much out of my depth... I'm using R to analyse a data set, which is a 1% sample of a population showing cancer rates across 20 regions. There are parameters for male and female rates. If I have to calculate 95% credible intervals for the mean based on the two sexes, am I better off using an arithmetic calculation or bootstrapping? I've searched high and low for something comprehensible but I'm struggling. I think I'm right in saying that bootstrapping is better when the data isn't normally distributed (this, however, very much is) and sample size is small (which this is). On the basis of the sample size, am I justified in going down the bootstrapping route? Any advice much appreciated!


title: Right-censored survival fit with JAGS

text: **Update: I got the JAGS model running and this eliminates the distracting part of my question. It's really about the proper preparation of data for dinterval() and inits.** I can't find a concrete example of preparing right-censored survival data for use in JAGS when starting with a vector of times and a vector of censored status. The examples I've seen all have a constant censoring time. My two questions are:   1. What is the correct way to prepare dinterval(t[i], c[i]) when you are starting with a vector of event times and a matching vector of censored/failed status?   2. What is the reason for supplying a initial value vector of event times where _uncensored_ times are set to NA and _censored_ times are set to a value > than the observed censoring time? I'm really struggling to understand that. **Update: I think I have the answer to (2). Apparently JAGS complains if you try to overwrite values in the init vector, but if the vector is initialized to NA it works.** Note that my understanding of these requirements comes from looking at existing online examples. If the requirements are unfounded I'd love to learn that. I created a short test using the aml data set from the survival package which does a survival fit then tries to fit a Weibull distribution using JAGS and add those points for comparison. Originally I thought you had to pass zeros in for the unscensored times in c[i] for dinterval(t[i], c[i]), but that causes an error: > Error in node is.censored[1] Observed node inconsistent with observed > parents at initialization. Try setting appropriate initial values. Here is my sample R script. Some lines seem to be indented too much, not sure why. I converted all tabs to spaces and it looks clean in my edit window:               # Right-censoring test     require(survival)     require(runjags)     graphics.off()     rm(fit, S, results)          # SURVIVAL FIT AND PLOT     fit <- survfit(Surv(aml$time, aml$status) ~ 1, conf.type = "log-log")     plot(fit, xlab = "time", ylab = "survival", main = "AML survival")          # JAGS WEIBULL FIT WITH CENSORING     # Model uses three different versions of observed times.     # t.orig is the original set of observed times including censored times.     # t.mod has the censored times set to NA.     # t.cens is the set I have questions about. Currently the same as t.orig     # but is that the proper choice?     #     model <- "model {         for(i in 1:N) {             is.censored[i] ~ dinterval(t.mod[i], t.cens[i])             t.mod[i] ~ dweib(a, b)             S[i] <- 1 - pweib(t.orig[i], a, b)         }         a ~ dgamma(3.5, 2)         b ~ dgamma(1.5, 100)     }"          censored <- !aml$status         t.orig <- aml$time     t.mod <- aml$time         t.mod[censored] <- NA         t.cens <- aml$time     # Online examples seem to indicate special treatment of the c[i]     # passed to dinterval(t[i], c[i]), changing the uncensored times.     #t.cens[!censored] <- 0 # THIS CAUSES AN ERROR          data <- dump.format(list(         t.orig = t.orig,         t.mod = t.mod,         N = length(t.mod),         t.cens = t.cens,         is.censored = as.numeric(censored)       ))          # Apparently you must provide initial values of the time vector     # with uncensored times of NA and censored times > the observed time.     # I don't understand this requirement at all.     t.init <- rep(NA, length(t.mod))     t.init[censored] <- t.cens[censored] + 1          inits <- c(         dump.format(list(a = 1.1, b = 3.0e-4, t.mod = t.init)),         dump.format(list(a = 2.3, b = 0.02, t.mod = t.init))     )     results = run.jags(model, monitor=c("S", "a", "b"), data=data, inits=inits)     S <- print(results, vars="S")     points(aml$time, S[,4], col="blue", pch=19)         segments(aml$time, S[,1], aml$time, S[,3], col="blue", lwd=2)      Thanks for taking the time to look at this.


title: When can I use STL decomposition?

text: Let us say I extract all tweets containing a certain hashtag (main subject of the tweet) and bin them by the hour i.e. count the number of tweets that occurred during a certain hour. I will then be obtaining a time series for this particular hashtag. I have recently come across the STL decomposition which decomposes a given timeline into its seasonal, trend and error components and `R` provides a simple function `stl` to achieve this. I am trying to understand if using this in my context (timeline of tweet counts) is appropriate. What I would like to examine is whether two given hashtags have a similar timeline and I was wondering what I can understand by using time series decomposition. Any suggestions?


title: Can I compute cross correlation after decomposition?

text: Given two time series within the same time period, I can use the `ccf` function in `R` to compute a notion of similarity of these time series. What I am wondering is this - can I compute the cross correlation on the trend line obtained after using time decomposition function such as `stl`? Time series decomposition is known to decompose a given series into its trend line, seasonal component and the error so by removing the error, would we get a better sense of cross correlation? So in short, my question is this: **_Is cross correlation of the trend line better than cross correlating raw timeseries?_** I tried explaining what I had done in the picture that follows. In the cross correlation graph that I obtained, it indicates a high value at t=0. ![enter image description here](http://i.stack.imgur.com/R82xE.png)


title: Repeated measurement analysis advice in R lme and / or aov and a couple of more questions

text: I want to analyse the visual performance, operationalized via contrast threshold, depending on adaptation luminance and spectrum, gathered for 29 subjects. I'm currently kind of confused of how to do that. Some of my data:                   ID C_measured Subject  LB Spectrum SpectrumLB       1  1 0.1339795   AHI11 0.1       HS     HS.0.1       2  2 0.1040440   AIC19 0.1       HS     HS.0.1       3  3 0.1363313   AUO13 0.1       HS     HS.0.1       4  4 0.1134103   BAR01 0.1       HS     HS.0.1       5  5 0.1117670   BAR02 0.1       HS     HS.0.1       6  6 0.1166350   BCL10 0.1       HS     HS.0.1      LB can be 0.1, 0.21, 0.3 and 1.0, Spectrum HS and LED I know that I want to do 8 planned comparisons of Spectrum at each LB (4x) and the effect of reducing LB from 0.3 to 0.21 for both spectra with and without exchanging the spectrum (4x). Some said I should do a two-way ANOVA, e.g. like                  aov(C_measured ~ LB * Spectrum + Error(Subject / (LB * Spectrum)), data = anovaFrameWithGlareFoveal)      and then a post-hoc test on the interaction variable SpectrumLB, e.g. like               anovaFrameWithGlareFoveal.lme <- lme(C_measured ~ LB, random = ~1 | Subject / LB, data = anovaFrameWithGlareFoveal)     anovaFrameWithGlareFoveal.glht <- glht(am2.subject, linfct = mcp(LB = "Tukey"))     summary(anovaFrameWithGlareFoveal.glht, test = adjusted(type = "none"))      then manually bonferroni-adjusting the p-value to the number of my planned comparisons. I do that for a couple of other parameters (with glare, without glare, old reference group, young group), which I don't want to include in the statistical analysis, because it is common knowledge that this influences visual performance, I just want to analyse whether the planned comparisons differ for those parameters. Till here the question: is this how to do it? DO I need the two-way ANOVA at all? Then I observed some things: for some of the parameters I observed significant values in the two-way ANOVA for the main effect of Spectrum p<.05, but none of the uncorrected multiple comparisons between the two spectra at all four LBs was significant (ok one was <.1, but I'm testing against .05), which some people commented with "impossible". Is this possible? Then people said: "ok paired t.tests should come up with the same results" so I did this:               df <- anovaFrameWithGlareFoveal     df.led <- subset(df, Spectrum == "LED")     df.hs <- subset(df, Spectrum == "HS")     t.test(df.led$C_measured[df.led$LB==.1], df.hs$C_measured[df.hs$LB==.1], paired=T)     t.test(df.led$C_measured[df.led$LB==.21], df.hs$C_measured[df.hs$LB==.21], paired=T)     t.test(df.led$C_measured[df.led$LB==.3], df.hs$C_measured[df.hs$LB==.3], paired=T)     t.test(df.led$C_measured[df.led$LB==1], df.hs$C_measured[df.hs$LB==1], paired=T)      then all of the t.tests were significant (<.05) but only one of the uncorrected multiple comparison was significant (<.05). I'm not really that much into statistics that I can definitively argue for or against one method or combined methods or whether the lme + glht alone is sufficient. Had a tough time on that the last week and am really curiously looking forward to your comments!


title: Ridge regression parameter estimation in R

text: I have seen about 16 parameter estimation methods in different research papers. How can these methods be implemented in R language using Monte Carlo simulation technique.


title: Interpreting coefficients of first differences of logarithms

text: My problem is interpreting coefficients of such time series model: \begin{equation} \ln Y_t - \ln Y_{t-1} =b_1 \cdot \left(X_{t}-X_{t-1}\right)+b_2 \cdot Z_t.\end{equation} I don't know how to interpret coefficients $b_1$ and $b_2$. Hope that someone can help.


title: Finding the best two predictor variables used conjointly, and levels of each

text: Twenty possible predictor variables in data set. One outcome variable. Some of the predictor variables are not linear. So a standard linear multiple regression approach probably won't do. (And I do not want to transform the variables using some arbitrary "cubic," "quadratic" approach -- leave the underlying variables as they are.) Example: Time of day. This oscillates regularly, over 24 hours. (Like a sine wave.) There is considerable random chaos inherent in some variables. Example: Highway traffic frequency. Some curve fitting may be advisable to dampen or smooth out the chaos. I'm looking for possible statistical packages (and specific procedures) that can accomplish the goal -- determining the two best predictor variables taken together of the twenty, and levels of each, that appear to best predict the outcome variable. Considering all twenty variables, there are 190 possible combinations. Doing an analysis of each possible pair would take a great deal of effort. Best to find a statistical package that can run tests on each possible combination in one long go. In summary: **Software Package** and **Specific Procedure** please. No need for a bunch of code. (If I want the code I'll ask in a follow-up.) Thanks. Nicholas Kormanik


title: k-subset with maximal variance

text: I have two versions of the same question:   1. Given a list of numbers (with possible duplicates), how to find a k-subset (with possible duplicates) that maximize the variance? is there a more efficient way than the obvious "check-all-k-subsets"?   2. Given a set of numbers, how do I select from that set a list of k numbers that maximize the variance?


title: Categorical variables and continuous response?

text: Within three distinct groups (three soil types), I have three categorical predictor variables (pH, ionic strength, and SAR). My response variable is continuous (% optical transmission). What statistical analysis do I need to check for interactions between the predictor variables and significant differences in the response? I also want to compare between soil types, and how the effects of the predictor variables on the response variable vary among soil types. Thanks.


title: Finding the Bayesian classifier for a bivariate Gaussian distribution

text: > Very close to: Joint Gaussian of two Gaussians I am trying to find the Bayesian classifier for two classes given by the following bivariate Gaussian distributions: $$p(x|c_1) = N(\mu_1, \Sigma_1)$$ $$p(x|c_2) = N(\mu_2, \Sigma_2)$$ What is the formula and the methodology I should use to do this? I am confused with out a proper understanding on a formula to do that at the moment. Prior probabilities of $c_1$ and $c_2$ are 0.5 respectively. FYI:               μ1 = ( 2        Σ1 = (2 0                        2 )            0 1)               μ2 = ( 2        Σ1 = (4 0                        4 )            0 2)


title: If only a two-way ANOVA main effect is significant, what is the appropriate follow-up?

text: What (with justification) is a valid post-hoc for a two-way ANOVA main effect, if no interaction is present? Example two-way fixed effect ANOVA:               > aov.example <- aov(Response ~ IV1 * IV2, data=data)     > summary(aov.example)                        Df Sum Sq Mean Sq F value   Pr(>F)        IV1            1  13.10  13.099  0.7222 0.40547       IV2            4 315.56  78.891  4.3498 0.01081 *     IV1:IV2        4 141.00  35.251  1.9436 0.14240       Residuals     20 362.74  18.137               ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I see a few general options, but I'm unsure which is most appropriate:   1. **Post-hoc with a multiple comparison test using two-way model**              > TukeyHSD(aov.example, which="IV2")          Tukey multiple comparisons of means        95% family-wise confidence level          Fit: aov(formula = Response ~ IV1 * IV2, data = data)          $IV2                    diff        lwr        upr     p adj     B-A      -2.1485711  -9.506162  5.2090200 0.9031415     C-A      -2.3382727  -9.695864  5.0193184 0.8733517     D-A       1.4732257  -5.884365  8.8308168 0.9735981     E-A      -8.0515205 -15.409112 -0.6939294 0.0277241     C-B      -0.1897016  -7.547293  7.1678895 0.9999912     D-B       3.6217968  -3.735794 10.9793879 0.5905067     E-B      -5.9029494 -13.260541  1.4546416 0.1559460     D-C       3.8114984  -3.546093 11.1690895 0.5439632     E-C      -5.7132478 -13.070839  1.6443432 0.1785266     E-D      -9.5247462 -16.882337 -2.1671552 0.0074740        2. **Post-hoc with a multiple comparison test using one-way model**              > TukeyHSD(aov(Response ~ IV2, data=data))          Tukey multiple comparisons of means        95% family-wise confidence level          Fit: aov(formula = Response ~ IV2, data = data)          $IV2                    diff        lwr        upr     p adj     B-A      -2.1485711  -9.858179  5.5610365 0.9224401     C-A      -2.3382727 -10.047880  5.3713349 0.8976378     D-A       1.4732257  -6.236382  9.1828333 0.9794638     E-A      -8.0515205 -15.761128 -0.3419130 0.0375667     C-B      -0.1897016  -7.899309  7.5199060 0.9999933     D-B       3.6217968  -4.087811 11.3314044 0.6456949     E-B      -5.9029494 -13.612557  1.8066581 0.1951992     D-C       3.8114984  -3.898109 11.5211060 0.6014671     E-C      -5.7132478 -13.422855  1.9963597 0.2212053     E-D      -9.5247462 -17.234354 -1.8151387 0.0102178      Overall, in this example the significant comparisons do not change, but the adjusted p-values are lower using the two-way ANOVA model. Which is most appropriate for a two-way main effect post-hoc?   3. **Post-hoc with subset levels of other independent variable** I'm pretty sure this is not valid, given the lack of any effect across the other independent variable. However, when one level of IV1 has larger effects between levels in IV2 compared to other IV1 levels, wouldn't this strongly skew the overall analysis? In other words, even though overall there is a main effect, could this effect not be significant in the the other level of IV1?                   > aov.level1 <- aov(Response ~ IV2, data=data[1:15,])         > summary(aov.level1)                          Df Sum Sq Mean Sq F value  Pr(>F)           IV2          4 334.55  83.639  3.8413 0.03837 *         Residuals   10 217.74  21.774                                > aov.level2 <- aov(Response ~ IV2, data=data[16:30,])         > summary(aov.level2)                          Df Sum Sq Mean Sq F value Pr(>F)           IV2          4 122.01  30.503  2.1037 0.1551         Residuals   10 145.00  14.500       Level 1 of IV1 has a significant IV2 main effect, but Level 2 does not. There are also differences in significant multiple comparisons. My concern is that I'm committing a Type I error using the previous two methods.                   > TukeyHSD(aov.level1)              Tukey multiple comparisons of means             95% family-wise confidence level              Fit: aov(formula = Response ~ IV2, data = data[1:15, ])              $IV2                          diff        lwr         upr     p adj         B-A       -7.86237771 -20.401217  4.67646157 0.3054028         C-A       -7.90634218 -20.445181  4.63249709 0.3008150         D-A       -0.94962690 -13.488466 11.58921237 0.9989922         E-A      -12.55848654 -25.097326 -0.01964727 0.0496014         C-B       -0.04396448 -12.582804 12.49487479 1.0000000         D-B        6.91275080  -5.626088 19.45159008 0.4167559         E-B       -4.69610883 -17.234948  7.84273044 0.7342625         D-C        6.95671528  -5.582124 19.49555455 0.4111062         E-C       -4.65214436 -17.190984  7.88669492 0.7404793         E-D      -11.60885964 -24.147699  0.92997964 0.0729541              > TukeyHSD(aov.level2)              Tukey multiple comparisons of means             95% family-wise confidence level              Fit: aov(formula = Response ~ IV2, data = data[16:30, ])              $IV2                        diff        lwr       upr     p adj         B-A       3.5652355  -6.667185 13.797656 0.7795224         C-A       3.2297968  -7.002624 13.462218 0.8321507         D-A       3.8960783  -6.336343 14.128499 0.7231230         E-A      -3.5445545 -13.776975  6.687866 0.7829174         C-B      -0.3354387 -10.567860  9.896982 0.9999634         D-B       0.3308428  -9.901578 10.563264 0.9999654         E-B      -7.1097900 -17.342211  3.122631 0.2257310         D-C       0.6662815  -9.566139 10.898702 0.9994435         E-C      -6.7743513 -17.006772  3.458070 0.2618999         E-D      -7.4406328 -17.673054  2.791788 0.1942037      This post has some relevant information.


title: Testing the significance of differences in frequencies across treatments using $\chi^2$ test

text: I've run a small experiment and need to test the data by comparing the frequencies of different treatments, specifically, whether they are significantly different. I think the chi-squared test is appropriate, but I am not sure how to execute it because it's somewhat outside the range of my course. Simplifying, the data looks something like this:                              A Choices (%)       B Choices (%)       Observations     Treatment 1:       40                   60                  n=30          Treatment 2:       65                   35                  n=30          Treatment 3:       70                   30                  n=30      How exactly would I conduct a chi-sq test? Is it the best test for investigating differences in the frequencies of 'A' between treatment pairs? Say, testing H0: Treatment1(% A Choices) = Treatment3 (% A Choices) Is it also possible to do a test like H0: Treatment2(% A Choices) = Treatment3(% A Choices) > Treatment1(% A Choices) Thanks!


title: What method is suitable for short-term forecast for a trendless, oscillatory, bounded time series?

text: I am new to time series analysis and I would appreciate if anyone could provide me some insight on it. I am trying to analyse a past series of numbers that fluctuates between 107 & 210 with a normal frequency bell curve distribution of mean 162. What is a suitable approach to forecast short-term future range for a trendless but oscillatory, range bound type of time series?


title: Online fitting for normal distributions

text: I was wondering if there exist efficient online or dynamic algorithm for fitting a normal distribution to data as it comes in. I am interested in two variants:   1. The algorithm is fed data points one at a time, and has to update its previous best-fit to account for the new point at each step.   2. The algorithm is fed $n$ data points initially. At each time step one data point is removed (the oldest one of the $n$) and one data point is added. The algorithm has to efficiently update the best-fit. Bonus points if you know an implementation of this algorithm in Matlab. Note, that the default Matlab normfit is not online, and it would be too computationally intensive to have to refit the data with it at every timestep.


title: Combining ratings from multiple raters of different accuracy

text: I have a data set of J raters who each give ratings to I objects (on say a scale of 1 to 5). My goal would be to construct some kind of overall rating for each i based on all of the raters' scores. My first approach, which I have seen used before, is to construct a standardized score specific to each rater-object pair so that each of j's ratings are standardized and j's rating for object i is: (j's score on i - j's avg score on all I)/ (standard deviation of j's scores on all I). If all raters were equally trusted, I think this would be an acceptable approach and has a nice interpretation. However, I suspect that there are lazy or poorly calibrated raters whose ratings are not spread out (e.g., they give 95% 4s and an occasional 5 when the distribution should really be uniform). It seems like I would want to somehow give greater weight to those raters with more spread out ratings. I also have an idea that a small subset of raters are "better" raters. Is there a way to perhaps weight each rater's score in line with their inter-rater reliability with the more accurate raters? For simplicity, my preference would be to use some kind of generalized linear model, if possible.


title: How are margins of error related to confidence Intervals?

text: Can somebody tell me the difference between margins of error and confidence intervals? On the Internet I see these two meanings getting used interchangeably. Is it right to say, "Confidence intervals are shown as 1.96 and displayed on the graphs as error margins"?


title: Multinomial logit model

text: Why do we need to use a base category (normalised to zero) when working with multinomial logit models? I would also like to ask why do we need to report conditional logit coefficients? Wouldn't marginal effects give us a better idea of the probability of belonging to a certain class?


title: Correct equation for weighted unbiased sample covariance

text: I'm looking for the correct equation to compute the weighted unbiased sample covariance. Internet sources are quite rare on this theme and they all use different equations. The most likely equation I've found is this one: $q_{jk}=\frac{\sum_{i=1}^{N}w_i}{\left(\sum_{i=1}^{N}w_i\right)^2-\sum_{i=1}^{N}w_i^2} \sum_{i=1}^N w_i \left( x_{ij}-\bar{x}_j \right) \left( x_{ik}-\bar{x}_k \right) .$ From: https://en.wikipedia.org/wiki/Sample_mean_and_sample_covariance#Weighted_samples Of course, you have to compute the weighted (unbiased) sample mean beforehand. However, I have found several other formulas like: $q_{jk}= \frac{1}{\sum_{i=1}^N w_i)-1}\sum_{i=1}^N w_i \left( x_{ij}-\bar{x}_j \right) \left( x_{ik}-\bar{x}_k \right) .$ Or I've even seen some source codes and academic papers just using the standard covariance formula but with the weighted sample mean instead of the sample mean... Can someone help me and shed some light? /EDIT: my weights are simply the number of observations for a sample in the dataset, thus weights.sum() = n


title: Confusion about hidden Markov model

text: I've gone through Hidden Markov models (HMM) for the past few months. However there are a few things that are confusing. The set up is simple: I have to model some human gestures such as walking, jumping, and falling. The observed data have been obtained via an accelerometer while the person was doing the movements. I trained theses observations using the famous Baum-Welch algorithm to get the parameters of an HMM for some states. Further, using the Forward and Backward procedures, the likelihood of the observation sequences given the model (i.e., the parameters) were found. Using a model selection criteria such as Akaike information criterion (AIC), I got the optimum states that represented the data:               a)Walking: 2 states     b)jumping: 2 states     c)Falling: 4 states      All these HMMs are then stored in a directory. Lastly, Viterbi decoding is used to get the most likely sequence of hidden states that produced the data. My questions are:   1. Suppose I performed the experiment again and I just get the data without knowing what kind of movement has been done. After getting the data trained, I got 2 states. How will the machine differentiate which kind of movement has been done, especially if walking and jumping are represented by 2 states?   2. Suppose the person has performed a different kind of gesture, e.g., sliding, what is the expected output after training? Will the machine be able to detect that or generate a false negative result?


title: How to compute an accuracy measure based on RMSE? Is my large dataset normally distributed?

text: I have several datasets on the order of thousands of points. The values in each dataset are X,Y,Z referring to a coordinate in space. The Z-value represents a difference in elevation at coordinate pair (x,y). Typically in my field of GIS, elevation error is referenced in RMSE by subtracting the ground-truth point to a measure point (LiDAR data point). Usually a minimum of 20 ground-truthing check points are used. Using this RMSE value, according to NDEP (National Digital Elevation Guidelines) and FEMA guidelines, a measure of accuracy can be computed: Accuracy = 1.96*RMSE. This Accuracy is stated as: "The fundamental vertical accuracy is the value by which vertical accuracy can be equitably assessed and compared among datasets. Fundamental accuracy is calculated at the 95-percent confidence level as a function of vertical RMSE." I understand that 95% of the area under a normal distribution curve lies within 1.96*std.deviation, however that does not relate to RMSE. Generally I am asking this question: Using RMSE computed from 2-datasets, how can I relate RMSE to some sort of accuracy (i.e. 95-percent of my data points are within +/- X cm)? Also, how can I determine if my dataset is normally distributed using a test that works well with such a large dataset? What is "good enough" for a normal distribution? Should p<0.05 for all tests, or should it match the shape of a normal distribution? * * * I found some very good information on this topic in the following paper: http://paulzandbergen.com/PUBLICATIONS_files/Zandbergen_TGIS_2008.pdf


title: How to get started on multivariate predictive modeling?

text: I'm used to doing descriptive modeling where I'm very concerned with parsimony and interpretability, but now I have a project that calls for the opposite and I'd like advice on techniques I should look into. I have a tool (an R function) that essentially creates random graphs based on 11 parameters. There are 4 statistics about the resulting graph that I care about--the number of nodes, the mean and standard deviation of the degree distribution, and the coefficient of variation of the area distribution (areas are assigned to each node as part of the R function). I'd like to input desired ranges for the statistics and find a set (or sets) of parameter values parameters that has relatively high probability of generating an output graph with those characteristics. I've run simulations (and running more is no problem, but they do take 1-2 seconds each) to generate lots of data, and now I need to model. I don't care what model parameters mean, I don't care how parsimonious the model is, I just need predictive power. My questions: (1) Where do I start? Should I be looking at GAMs or other non- parametric regression techniques? (2) I've been assuming that I would fit the 4 graph statistics as dependent on the parameter values. It occurred to me recently that I could do it the other way around. Would that be reasonable?


title: Plot of copula (based on data set) - R

text: I have to do an empirical analysis for a statistics paper. For this I want to show the differences of dependence structure for a specific data set. So I selected 2 stock prices, transformed them into the returns and started to measure the dependency with R. So far it is no problem, I have a result for Bravais-Pearson, Kendall and Spearman. Additionally I plotted the regression model for this two values. I have read in many papers, according to Sklar's Theorem, that it is easy to get the copula function out of the distribution function, just by use the inverse. So my question is, if there is a possibility with R to plot the copula function (and density) just by having this data set (2 returns) or if I must first estimate the parameters to be able to plot this function. And how can I do this with R? I tried to search the answer in the handbook of the package "copula", but my search wasn't really helpful. Thanks in advance for your help!


title: How to find parameters in multivariate space efficiently?

text: I am trying to optimize sklearn.linear_model.SGDRegressor, and I was wondering if people could point me in which direction I should try to optimize? Personal experience and literature are both welcome. I think it would be much better if there is kind of a ranking in which way to search; should I start with finding a good alpha, or with a learning rate? Do you also have to do back checks at some point? What is the common way of searching this specific (or other sklearn/machine learning related) multivariate space? How dependent are these parameters on each other? I understand this matters per situation, but how to do this kind of optimizing in a smarter way than level wise search? I'm also never so sure when changing anything could have an effect on these paramaters. Basically if I would take a smaller sample, then I have no idea how solid the "optimal" values I found before, are. E.g. some possible parameters to tune:               loss : ‘squared_loss’, ‘huber’, ‘epsilon_insensitive’, or                    ‘squared_epsilon_insensitive’     alpha : float     l1_ratio: float between 0 (L2) and 1 (L1)     epsilon: float     learning_rate : optional     eta0 : double, optional     power_t : double, optional


title: How to check if there is any “dependence” between columns?

text: How to check if those values are "dependent" ? Take a look on the values in second column coz the difference between them is really huge so it's hard to choose the scale. Any propositions of checking "dependence" between them are welcome. I was thinking about correlation but I am new in statistic so please let me know step by step what should I do to analyse those data.               diff    abund_mean      0   3444804.79      1   847887.02      2   93654.19      0   721692.76      2   382711.04      1   428656.65      1   120933.91      0   157528.72      1   159650.70      0   124602.80      0   90844.33      2   501825.37      1   270592.56      I am learning R aswell so I can calculate evertyhing in R if you let me know how to do that.


title: Construct confidence interval of the mean for auto-correlated data

text: I feel like I'm missing something obvious, but here we go. I have auto- correlated data measured in triplicate for two (or more) treatments. Something like this:               t <- 3:20 #times in my real dataset are possibly not always equidistant     a <- structure(c(0.652492388457625, 0.905172522010166, 1.23437705454616,                       1.48003667490842, 1.77876898946135, 1.99175317367897, 2.31666502140984,                       2.43520651415548, 2.67903421794922, 2.84115747823017, 2.89693734873647,                       2.91199679761145, 2.85645436179354, 2.99371033437697, 2.99965220711105,                       2.84984814715963, 2.64275376547326, 2.64060469520379, 0.481029734912324,                       0.8466803252367, 1.31126162780809, 1.56745630574946, 1.74865844658142,                       1.80367117155375, 2.06688393210808, 2.24500095501872, 2.52978288460243,                       2.69073206006205, 2.89657418056785, 2.93759772556246, 2.99305951550274,                       2.89146932307489, 2.88890777189028, 2.7974672802907, 2.70933381639295,                       2.66799551352975, 0.624178180970784, 0.867127935268765, 1.09752295578438,                       1.35037796202753, 1.60094288950107, 1.97949255710341, 2.15496378191076,                       2.42556913246041, 2.54331160179646, 2.67440414122285, 2.84249532365163,                       2.95278639560433, 3.06192227561515, 3.03297885461444, 3.04101341059534,                       3.01736966686846, 2.80061410999215, 2.69852643323913),                     .Dim = c(18L, 3L), .Dimnames = list(NULL, c("a1", "a2", "a3")))     b <- structure(c(0.516527990622755, 0.84883434472028, 1.04202664437099,                       1.3100841689546, 1.48050413266838, 1.7824492800856, 1.96557179831706,                       2.17419105778186, 2.2453178060978, 2.35460428313729, 2.49308342865959,                       2.62343038370418, 2.70831189685371, 2.79459971623943, 2.94938536147398,                       3.04822554887815, 3.00287042052314, 2.91673487674283, 0.589490441973075,                       0.751768045201717, 0.917973959434798, 1.17617337222852, 1.39497560590896,                       1.65920945485901, 1.87749014780468, 2.11880355292648, 2.372755207219,                       2.46211141942227, 2.59688733749884, 2.72270421752644, 2.79848710425447,                       2.81134394947587, 2.75390203306788, 2.78499114431362, 2.86001341271914,                       2.95652300178809, 0.558662398944567, 0.834996005844121, 0.988238211915554,                       1.27569591423003, 1.38577342414377, 1.62664982549252, 1.83299700801392,                       2.04943560731628, 2.22950648854987, 2.38533269800646, 2.49845003387994,                       2.60036098089373, 2.61941602504858, 2.71298500309883, 2.78126388719353,                       3.04792375845498, 3.02691814463875, 3.06667590650438),                     .Dim = c(18L, 3L), .Dimnames = list(NULL, c("b1", "b2", "b3")))          matplot(t,a,pch=1,xlab="",ylab="",col="blue")     matlines(t,a,col="blue", lty=2)          matpoints(t,b,pch=16,col="red")     matlines(t,b,col="red", lty=2)      ![simulated auto-correlated data](http://i.stack.imgur.com/RbEPX.png) I would like to know in which time periods the treatments differ. I would like to avoid fitting any kind of model. (There are models for my kind of data from science, but they are known to be only an approximation for some ranges of my data and I'm afraid that model error might mask differences.) My idea is to calculate the mean and construct confidence intervals (using an assumption of normality) like this:               a_means <- apply(a,1,mean)     a_sds <- apply(a,1,sd)     a_lwr <- a_means-qt(0.975,3)*a_sds/sqrt(3)     a_upr <- a_means+qt(0.975,3)*a_sds/sqrt(3)          b_means <- apply(b,1,mean)     b_sds <- apply(b,1,sd)     b_lwr <- b_means-qt(0.975,3)*b_sds/sqrt(3)     b_upr <- b_means+qt(0.975,3)*b_sds/sqrt(3)          DF <- data.frame(treat=factor(rep(1:2, each=length(t))),                       time=rep(t, 2),                      mean=c(a_means,b_means),                      lwr=c(a_lwr,b_lwr),                      upr=c(a_upr,b_upr))          library(ggplot2)     p <- ggplot(DF, aes(x=time, y=mean, ymin=lwr, ymax=upr)) +       geom_ribbon(aes(fill=treat), alpha=0.3) +       geom_line(aes(color=treat))     print(p)      ![means and 95 % confidence intervals](http://i.stack.imgur.com/2SUyM.png) The way I'm constructing the confidence intervals obviously doesn't consider auto-correlation.   * Is there a way to construct some kind of "auto-correlated confidence interval"?   * Can I use the "un-correlated confidence interval"? Can I somehow estimate if it is too narrow or too wide in comparison to the auto-correlated confidence interval?   * Is there a better approach to my problem?


title: Simulating p-values as a function of sample size

text: We are trying to prove a very subtle effect occurring in cells after a certain treatment. Let's assume that the measurements are normally distributed. Let's also assume the untreated cells have $\mu = 1$ and $\sigma = 0.1$ and the treated cells have $\mu = 1.1$ and $\sigma = 0.22$. The question is: How large must the sample size be in order for the observed effect to be statistically significant ($\alpha = 0.05$)? I know that very subtile effects require a larger sample size than more apparent effects, but how many? I'm still learning statistics, so please be patient with me. I tried to perform a little simulation in R. Assuming that you randomly pick $n$ samples from a normal distribution, I tried to calculate the mean p-value as a function of $n$. ![enter image description here](http://i.stack.imgur.com/Fn5d3.png) Is this a correct way to do find the right sample size? Or am I completely off the track with this approach? Code:               library(ggplot2)          ctrl.mean <- 1     ctrl.sd <- 0.1     treated.mean <- 1.1     treated.sd <- 0.22          # Function that repeats t-test a number of times (rpt) with given sample size, means and sds.     # Returns a list of p-values from the test          tsim <- function(rpt, n, mean1, sd1, mean2, sd2) {       x <- 0       ppool <- NULL       while (x <= rpt) {         ppool <- c(ppool, t.test(rnorm(n,mean1,sd1), y = rnorm(n,mean2,sd2))$p.value)         x <- x + 1       }       return(ppool)     }          # Iterate through sample sizes and perform the function     # Returns data frame with list of mean p-values at a given sample size          i <- 2     num <- 50     res <- NULL          while (i <= num) {       sim <- tsim(1000, i, ctrl.mean, ctrl.sd, treated.mean, treated.sd)       res <- rbind(res, cbind(i, mean(sim), sd(sim)))       i <- i + 1     }          # Plot the result          res <- as.data.frame(res)          ggplot(res, aes(x=i, y=-log10(V2))) +       geom_line() +       geom_ribbon(aes(ymin=-log10(V2)-log10(V3), ymax=-log10(V2)+log10(V3)), alpha = 0.2) +       annotate("segment", x = 6, xend = num, y = -log10(0.05), yend = -log10(0.05), colour = "red", linetype = "dashed") +       annotate("text",  x = 0, y=-log10(0.05), label= "p = 0.05", hjust=0, size=3) +       annotate("segment", x = 6, xend = num, y = -log10(0.01), yend = -log10(0.01), colour = "red", linetype = "dashed") +       annotate("text",  x = 0, y=-log10(0.01), label= "p = 0.01", hjust=0, size=3) +       annotate("segment", x = 6, xend = num, y = -log10(0.001), yend = -log10(0.001), colour = "red", linetype = "dashed") +       annotate("text",  x = 0, y=-log10(0.001), label= "p = 0.001", hjust=0, size=3) +       xlab("Number of replicates") +       ylab("-log10(p-value)") +       theme_bw()


title: Not specifying the main effect of a term that is part of a tensor product interaction

text: Say you've got a model $$ y = f(X_1,X_2)+\epsilon $$ and you're OK with linear (or other parametric) functional forms. Say you think that the effect of $X_1$ on $y$ depends on $X_2$. The standard wisdom is that one should include the main effect of the interaction to capture the effect of $X_1$ on $y$ when $X_2$ = 0, and to render coefficients invariant to location changes. I am curious whether this applies to generalized additive models estimated by penalized splines, with interactions represented by tensor products. Intuitively, it seems like it would NOT apply. WRT the first issue, the tensor product smooth estimates the "height" of a surface defined by two variables, along a grid. Thus you've always got a function for one variable at a fixed value of the other. If you had a univariate smooth along with it, you'd have to simply add the two, at a cost of degrees of freedom (which would get penalized/reduced). WRT the second issue, this paper by Simon Wood lays out how the tensor products implemented in `gam` are scale invariant. Is this a reasonable interpretation? If I wanted to estimate the model above nonparametrically, is it unreasonable to leave out smooth functions of $X_1$ and $X_2$, and simply model $X_1 \otimes X_2$. If yes, how would you quickly and imply explain why to an audience that isn't used to GAMs?


title: Increasing power of test by transforming format of data - is this cheating?

text: I'm working on a health care dataset that contains one row per member, measuring total number of the member's hospital inpatient (IP) visits during 2012 and 2013. The initial research question of the study asks if the percentage of high IP utilizers (3+ visits per year) in 2013 is statistically significant different from 2012. However, our power is limited by small sample of members in either year who had 3+ IP visits. Rather than comparing the % of **members** who are high utilizers in 2012 and 2013, some colleagues have proposed we ask a slightly different research question which increases the power of the study by comparing the percentage of total hospital **visits** accounted for by the high utilizers. Now they're no longer interested in the % of members who are high utilizers, but rather the % of all IP visits which are associated with members who had 3+ visits per year. Therefore, Member A with 9 IP visits, Member B with 2 IP visits, and Member C with 1 IP visits in 2012 are now measured as **12 observations** rather than 3 observations. _We've quadrupled our sample size._ The % IP visits belonging to a high utilizer is 100% x (9/12) = 75%. My question: Is this "cheating" - increasing the power of the study by inflating the number of records by the number of visits per member?


title: How to deal with omitted dummy variables in a fixed effect model?

text: I am using a fixed effect model for my panel data (9 years, 1000+ obs), since my Hausman test indicates a value $(Pr>\chi^2)<0.05$. When I add dummy variables for industries that my firms included, they always get omitted. I know there is a big difference when it comes to the DV (disclosure index) among the different industry groups. But I am not able to get them in my model when using Stata. Any suggestions how to solve this? And why are they omitted?


title: How do you calculate the expected value of $e^{-X}$?

text: If $X$ is a random variable, I would like to be able to calculate something like $$E(e^{-X})$$ How can I do this? Thank you so much.


title: Cluster clickstream data

text: I've recently entered the realm of machine learning and a project I am working on requires me to cluster users based on the order they visited webpages on a website. I have data in the form of:               ['user_id', 1, 2, 4, 6, 3, 7, 3, 2, 4...]      Where each number is a category/page that the user visited. In addition the length of data for each user is not the same i.e. some users visit more pages than others. I realize this is really vague and defining similarity it hard. I tried following the example in this research paper and to be honest a lot of it went over my head. I need help in how to approach this problem and am open to new ideas and suggestions.


title: How to compute the number of games that must be played for the better team to win 70% of the time?

text: I've been stuck on this sample statistics problem for a couple days now: Assume that two baseball teams have win records of $58\%$ and $54\%$. Demonstrate how to compute or estimate the number of games that must be played in order for the better team to win $70 \%$ of the time, $80\%$ of the time, and $90\%$ of the time. I _think_ the question is asking how many games must each team play against each other in order for the better team (team with the $0.58$ winning percentage) to win $X\%$ of the time. If I use a binomial distribution, the Central Limit Theorem and the $Z$ score formula I can work backwards and get a formula for the number of games. However, this approach only uses the winning percentage of the better team (i.e. $0.58$) which makes me think that I'm missing something important. Perhaps I'm misinterpreting the question? Any thoughts or suggestions for how to approach this problem would be greatly appreciated!!!


title: Evidence on red-purple-blue graphs

text: For some reason I have seen a preponderance of graphs whose color schemes (to represent a continuous variable) transition from red through purple to blue. Here is but one example (from a blog post, but I've seen plenty of published examples): ![http://www.r-bloggers.com/bubble-plots- ggplot2/](http://i.stack.imgur.com/nwTfN.png) Personally, I find it impossible to clearly associate the color change with an increase in magnitude of the plotted variable, but I realize color is by its nature subjective and the very fact that I've seen it so much makes me suspect that there's a reason for that. Is there a clear consensus on red-purple-blue transitions in the literature?


title: How to quantify the plastic deformation point?

text: In the figure below, the stress-strain relationship is linear until they reach their corresponding plastic deformation points (A,B,C,D), beyond these points the relationship deviates from the linearity. The corresponding plastic deformation points shown in solid circles were actually determined by visual inspection (ball-parked). I am wondering if there is any better way to determine or quantify these points with any confidence level, say 90%. ![enter image description here](http://i.stack.imgur.com/lG8iB.png) (The above question is the second part of a previous question posted in this forum related to nonlinear fitting Here) UPDATE: In the following first Figure, for each data set, I took first two points and fitted with a straight line. What was interesting that for each data set, the measured strain values were increasing on the left side of the linear fit and then gradually decreasing to the other side. ![enter image description here](http://i.stack.imgur.com/fbwLd.png) Then I took the difference between the measured and the predicted value, and plotted against the stress values as shown in the following figure. The points where the slopes changing are marked as in the figure in the original post. ![enter image description here](http://i.stack.imgur.com/1ryj6.png) But the same problem prevails,   1. I had to chose the point based on my visual inspection   2. Could not figure out how to choose confidence interval   3. Not sure whether the procedure followed has any statistical significance Any further comments would be highly appreciated. Thanks!


title: Difference-in-difference more than two periods

text: I try to use a difference-in-difference design in my study. But I have three periods. So I would like to know how I could run my difference-in-difference analysis on Stata with more than two periods.


title: Variance and covariance notation: $\sigma^2 V_1$, $\sigma^2 V_2$

text: I am reading an internal paper that says: > Let $\sigma^2 V_1$ equal the variance of $\sum_{m\in M}Z_m - Z_0$ and > $\sigma^2V_2$ equal the covariance of $||M||^{-1}\sum_{m\in M}Z_m - Z_0$ and > $Z_m - Z_0$. I am unfamiliar with the notation of the $V$s. Can anyone explain what this is referring to or point me in the right direction?


title: Why do both the VIF and tolerance statistics exist, when the latter is just the reciprocal of the former?

text: Is taking the reciprocal helpful in some way, or is it just a matter of historical accident that there came to be two terms to describe the same thing? ![enter image description here](http://i.stack.imgur.com/oaPJo.png) (From the Wikipedia page.)


title: When is covariance a reasonable measure of correlation?

text: I am brushing up on statistics for a data science course. I saw an example of variables that are highly dependent, but still have 0 covariance (http://en.wikipedia.org/wiki/Covariance#Uncorrelatedness_and_independence) Now ... what I want to know is: when would covariance be a good measure of correlation, and when not. What kind of data is suitable for analysis using covariance, and what isn't?


title: Calculating probability from missing/sparse sample

text: Given               score(Q, R) = score obtained by person R for topic Q, where smaller value is better      In an ideal scenario, I am given a set of data (events are mutually exclusive, $r_n$ scoring well in $q_1$ does not affect scores in $q_2$)               score(1, A) = 40     score(1, B) = 50     score(1, C) = 60     score(1, D) = 20     score(1, E) = 10          score(2, A) = 20     score(2, B) = 40     score(2, C) = 60     score(2, D) = 90     score(2, E) = 90          score(3, A) = 30     score(3, B) = 20     score(3, C) = 70     score(3, D) = 60     score(3, E) = 20      Then I can derive a probability to see the chances of another similar person obtaining a good score in Q as follows               P(good score for Q) = $good_score_count / $people_count      for example (assuming good score standard is 50 and below)               P(good score for Q=1) = 4 / 5     P(good score for Q=2) = 2 / 5     P(good score for Q=3) = 3 / 5      so I can use for further analysis. However, in real life, we are only able to collect incomplete data for some practical reason, for example               score(1, A) = 40     score(1, C) = 60     score(1, D) = 20          score(2, A) = 20     score(2, B) = 40     score(2, C) = 60          score(3, A) = 30     score(3, B) = 20     score(3, C) = 70     score(3, E) = 20      In this case it doesn't seem to make sense to calculate P(good score for Q) with all the missing data. Is there anything I can do to somehow guess how well a new person ($r_n$) does in each of the topics based on performance of other people?


title: Smoothing algorithm for irregular time interval

text: I have various sets of irregular interval time series data to which I want to apply some sort of smoothing algorithm to produce a good fit. I have attempted various methods which all were unsatisfactory.   1. Loess - Too much of a tendency to overshoot/overreact to outliers   2. Moving Average - The lag is unacceptable Example Dataset: ![http://i.stack.imgur.com/sFkTe.png](http://i.stack.imgur.com/7Q6Sm.png) I have read about the "Improved Holt Method for Irregular Time Series", but the paper was too difficult for me to understand and implement in C#. Can someone point me to a good method / algorithm which produces good smoothing? The method must be able to calculate the smoothed point at time $t$, without requiring $t+1$, etc., data. It also must be capable of dealing with multiple $y$ values for a given $x$ time.


title: How to use Almon model (polynomial distributed lag) in R program

text: I am testing the influence of intramural R&D expenditure, life expectancy, HDP and tax on the GINI index, but I have a problem because I am trying to make my model better by using the Almon model of distributed lag on R&D since the effect of it comes in later years. I am using R to do my analyses. These are my data:               Gini          <- c(33, 33, 33, 34, 32, 30, 29, 30, 30.6, 31.5, 31.9, 31.9, 31.3,                         29.9, 28.8, 33.2, 33.5)     VaR           <- c(1.25, 1.3, 1.27, 1.23, 1.18, 1.11, 1.09, 1.1, 1.16, 1.23, 1.25,                         1.25, 1.29, 1.46, 1.76, 1.71, 1.72)     Vek           <- c(75.46, 75.83, 76, 76.2, 76.1, 76.6, 77.2, 77.9, 78.3, 78.9,                         79.4, 79.7, 79.7, 80.2, 80.2, 81, 80.6)     Danova.sadzba <- c(28.13, 28.48, 25.99, 24.94, 24.34, 27.44, 23.04, 22.18, 22.24,                         23.27, 23.56, 21.54, 19.33, 18.95, 21.34, 21.75, 18.9)     HDP           <- c(11.4, 9.3, 11.5, 8.8, 11.1, 10.7, 5.3, 5.6, 3.9, 4.4, 5.9, 5.4,                         5.4, -2.1, -5.5, -0.8, 1.4)      Can anybody help me with what to do in R to make this work?


title: Pareto Optimal - improving for one player, static for another

text: This is very basic - Something is pareto optimal if you can't improve the outcome for one player without worsening the outcome for another. But what about where one player gets the same but the other player gets more? Is one pareto optimal? E.g.                          p1, p2     outcome 1: 2   2     outcome 2: 2   4      The move from outcome 1 to outcome 2 improves for p2 but it doesn't make the outcome worse for the other player. I would say that outcome 2 is pareto optimal and that outcome 1-> outcome 2 is a pareto improvement.


title: Analyzing reflected and transformed variables

text: I want to perform multiple linear regression on my variables; since my dependent variable was negatively skewed, I reflected it and then I log transformed it. My independent variables was positively skewed, so I only log transformed them. My first question is: **after performing multiple linear regression on my variables, should I interpret results adversely?** I mean: if an independent variable has a positive beta, I say it has a negative effect on dependent variable; and if an independent variable has a negative beta, I say it has a positive effect on dependent variable? My second question is: **my adjusted $R^2$ is $0.35$, is that good enough to say that my model is acceptable?**


title: XOR gate with Restricted Boltzmann Machine

text: I'm new to this field and trying to perform some sort of supervised learning with Restricted Boltzmann Machines. I am curious to see whether an RBM can be trained to do the job of an XOR gate? This is a very simple example that can be used to test a typical neural network.


title: Emotion recognition lip zone processing

text: I want to detect emotions by comparing the cropped lip zone with different templates and output an emotion based on the score. So far I've applied:   1. Normalization   2. Discrete and fast Fourier transform   3. Gaussian blur   4. Lip-zone cropping > ![enter image description here](http://i.stack.imgur.com/gpwrp.png) ![enter > image description here](http://i.stack.imgur.com/aXqLT.png) ![enter image > description here](http://i.stack.imgur.com/fh9wE.png) All of the above done in Python with OpenCV, NumPy, matplotlib... The next step would be applying some other operation in order to get the lip zone more similar to these: **DESIRED OUTPUT** > ![enter image description here](http://i.stack.imgur.com/y9HJ9.png) ![enter > image description here](http://i.stack.imgur.com/39168.png) ![enter image > description here](http://i.stack.imgur.com/FyCtZ.png) ![enter image > description here](http://i.stack.imgur.com/ti3BP.png) **ACTUAL OUTPUT** > ![enter image description here](http://i.stack.imgur.com/WEsq5.jpg) ![enter > image description here](http://i.stack.imgur.com/e9mGa.jpg) What other operations should I apply on the images above (ACTUAL OUTPUT) in order to:   * reduce noise   * obtain a contour of the lips   * and after that, apply binarization in order to prepare the picture for comparison.


title: Find infinite moving average representation of seasonal model

text: ![enter image description here](http://i.stack.imgur.com/8DhTt.png) I am attempting previous exams while studying for my own final exam, and can't do this question. I did a-c fine, but am having a lot of trouble doing parts d,e and f. I found a coefficient of 0.5 while equating like terms of the infinite moving average, but don't really know where the summations came from. Any advice would be greatly appreciated.


title: Interpreting the drop1 output in R

text: In R, the `drop1`command outputs something neat.   These two commands should get you some output:   `example(step)#-> swiss`   `drop1(lm1, test="F")` Mine looks like this:               > drop1(lm1, test="F")     Single term deletions          Model:     Fertility ~ Agriculture + Examination + Education + Catholic +          Infant.Mortality                      Df Sum of Sq    RSS    AIC F value     Pr(F)         <none>                        2105.0 190.69                           Agriculture       1    307.72 2412.8 195.10  5.9934  0.018727 *       Examination       1     53.03 2158.1 189.86  1.0328  0.315462         Education         1   1162.56 3267.6 209.36 22.6432 2.431e-05 ***     Catholic          1    447.71 2552.8 197.75  8.7200  0.005190 **      Infant.Mortality  1    408.75 2513.8 197.03  7.9612  0.007336 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       What does all of this mean? I'm assuming that the "stars" help in deciding which input variables are to be kept. Looking at the output above, I want to throw away the "Examination" variable and focus on the "Education" variable, is interpretation this correct? Also, the AIC value, lower is better, yes? Ed. Please note the Community Wiki answer below and add to it if you see fit, to clarify this output.


title: How model time-dependent covariates?

text: I have already understood that modeling time dependent covariates should be analog to his:               cfit <- coxph(Surv(tstart, tstop, status) ~ treat + sex + age +     inherit + cluster(id), data=cgd)      In my dataset, however, the dependent variable has no apart category for 'start' en 'stop'. Just a list of random numbers. I'll copy the upper part:               12     9     27     15     14     2     8     ...      The dependent variable is time until infection. This should be written as a function of age, gender, time until prophalctic treatment an time until excision (the last two are time dependent) Does anyone know how I should model this?


title: The impact of rescaling a predictor on the standard error of the corresponding coefficent

text: I am trying to form a polynomial regression model using SVD linear model. As the predictor at a large degree goes too large, say _x_ 6, I first scale it down if the mean of _x_ 6 is over a threshold and then put it in the design matrix _X_. I then can easily get the actual coefficient of _x_ 6 by multiplying the same factor. Now the problem is how to get the actual standard error of the coefficient of _x_ 6. From "Standard errors for multiple regression coefficients?" for simple regression, _y = a + b'x'_ if _x'_ is scaled down by _m_. From ${\rm var}(b')=s^2/∑(X_i'−\bar X')^2$ we can easily see the actual   var( _b_ ) = var( _b'_ )/ _m_ ² and hence _SE(b) = SE(b')/m_ , and also the standard error of intercept is not impacted. However, for multivariate regression with the general form of covariance matrix   cov( _β_ ) = _s_ ²( _X'X_ )-1. I can't see how to get the actual standard error back. Another way to look at the problem is from SVD calculation. as _X_ = _UWV'_ from SVD, and   ( _X'X_ )-1 = _VW_ -2 _V'_ , as long as we know how the scaling impacts the V and W, we also can know the impact to the standard error. However, I also have difficulties knowing the impact on _V_ and _W_. Any suggestion is appreciated. I have been struggling with the issue for quite a while.


title: I need some examples of data where there is an interaction and with/without main effects

text: I had thought that if there is an interaction between 2 independent variables, then it is always the case that there will be main effect for both the independent variables. However, I found out that this is not the case after reading this site. Can I have some practical examples where there is interaction and:   1. no main effects for both independent variables.   2. no main effect for one of the independent variables only.


title: Test for Complete Spatial Randomness taking into account background distribution

text: I have a set of points on a grid which is a subset of a larger background set of points on a grid. I would like to test if the subset points are randomly distributed (or not). If it helps, here is the real-world example I am working on: http://funcgen.vectorbase.org/expression- map/Anopheles_gambiae/VB-2013-12/?search1=GO%3A0005549%0A%0A The length of the yellow sausages is proportional to the number of genes in the subset, while the area of the grey circles is proportional to the number of genes at each grid point belonging to the total/background set of genes. In this case the yellow genes are clearly non-randomly distributed. The spatial distribution is important, I believe... imagine just three yellow genes "landing" in three neighbouring grid squares - I'm pretty sure that considering each grid point in isolation there is nothing special about the sampling - however the fact that the three genes' grid nodes are so close on the grid of 500 nodes **is** special/interesting, right? I think the R spatstat package can help me here but I am not sure how to take into account the "background distribution". Here's what I have managed so far - with dummy random data:               library(spatstat)          # here is the "background" distribution of 100 points on a discrete 10x5 grid     # (yes there are duplicate points at the same coordinates)     bg = ppp(sample(0:9, 100, replace=T), sample(0:4, 100, replace=T), window=owin(c(0,9),c(0,4)))          # if you want to visualise the density of points:     # make a bitmap image version     bgim = as.im(bg, dimyx=c(5,10))     # plot(bgim)          # now make a random subset of the points     subset = bg[sample(1:100, 10, replace=F)]          # plot bg with plus signs     plot(bg, pch=3)     # add the fg points      plot(subset, add=T)      So I can do test for CSR as follows               mad.test(subset)      but I don't know how to include the background distribution. I've tried fitting models with ppm but am out of my depth. Can anyone help please? If possible I'd like the test to be fast (1-2s max on subsets up to a few thousand genes out of the total of 12,500 genes). Writing this maybe I have a rough idea how to tackle this... Should I load the entire set of genes into a ppp (point pattern) and a "mark" the subset with a label (e.g. selected/unselected)? Then maybe some of the simple tests will then work straight out of the box? Think of it like a distribution of trees (the forest kind), and some are marked as infected with a disease. You would want to know if the disease is localised with respect to the background distribution of trees. **Update: my solution thanks to Andy W.** Pass an expression which generates random subsets of the background point set (the same size as the subset of interest of course)               result = mad.test(subset,                simulate=expression(bg[sample(1:bg$n,subset$n,replace=F)]),                nsim=1000)


title: How does glmnet handle overdispersion?

text: I have a question about how to model text over count data, in particular how could I use the `lasso` technique to reduce features. Say I have N online articles and the count of pageviews for each article. I've extracted 1-grams and 2-grams for each article and I wanted to run a regression over the 1,2-grams. Since the features (1,2-grams) are way more than the number of observations, the lasso would be a nice method to reduce the number of features. Also, I've found `glmnet` is really handy to run lasso analysis. However, the count number of pageviews are overdispersed (variance > mean), but `glmnet` doesn't offer `quasipoisson` (explicitly) or `negative binomial` but `poisson` for count data. The solution I've thought of is to `log transform` the count data (a commonly used method among social scientists) and make the response variable roughly follow a normal distribution. As such, I could possibly model the data with the gaussian family using `glmnet`. So my question is: is it appropriate to do so? Or, shall I just use poisson for `glmnet` in case `glmnet` handles `quasipoisson`? Or, is there other R packages handle this situation? Thank you very much!


title: Discarding the violation of homogeneity of variance in ANOVA

text: If after a `bartlett.test(split(x,list(f1,f2)))` I get a very significant _p- value_ , but my sample size is huge -- millions of measurements -- is it justified to discard the homogeneity of variance in case of a significant ANOVA?


title: Pseudo R equal to one in logistic regression

text: In Stata, how can pseudo _R_ equal one? This is accompanied by blank spaces in standard error and other indicators. The analysis is binary logistic regression. What can cause this?               ----------------------------------------------------------------------------                F | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------            BESIM |   1.194625          .        .       .            .           .        KOMUNIKIM |   .2921695          .        .       .            .           .         KOOPERIM |   4.659727          .        .       .            .           .         ANGAZHIM |   1.470268          .        .       .            .           .         KONTROLL |   .5014514          .        .       .            .           .           BURIME |   23.79765          .        .       .            .           .       KOMPETENCA |   .0178707          .        .       .            .           .            KOSTO |   3.393804          .        .       .            .           .          AGJENCI |   .6703374          .        .       .            .           .             DIJE |   6.720107          .        .       .            .           .            _cons |   6.86e-69          .        .       .            .           .     ------------------------------------------------------------------------------     Note: 5 failures and 7 successes completely determined.


title: expected shortfall and value-at-risk

text: I once read a R example of computing `Value-at-Risk` and `expected shortfall` as follows               p = c(0.5,0.1,0.05,0.025,0.01,0.001)          VaR = qnorm(p)          ES = dnorm(qnorm(p))/p      I am not sure why the expected shortfall can be calculated this way?


title: Expectation maximization with variant length at observing data

text: Imagine one loaded dice. Based on EM algorithm, how could we compute how much it loaded if we introduced:   1. Variant length on each rolling attempt (look at first and second attempt below 1st one has 6 attempts and second one has 3);   2. Uncertainty on seeing each sides value (see-quality, The value in parentheses means we are 10% sure about seeing this side)? To be more clear here is the data for 5 series of rolling one of them :               1--> (10%)6(23%)4(15%)3(1%)2(100%)1(9%)1          2--> (33%)3(55%)2(11%)1          3--> (0%)4(3%)4(16%)3(11%)2(26%)2(76%)1(33%)4(65%)3(76%)6          4--> (91%)4(0%)2(11%)5(11%)1          5--> (18%)3(44%)1(96%)5(11%)1(96%)2(16%)2(83%)1      Any suggestions would be appropriated.


title: Type of time series to use

text: Consider the following time series plot: ![enter image description here](http://i.stack.imgur.com/hrUTb.png) I am trying to fit a time series to this data. I am using it for the purposes of prediction. What type of time series would work here? I would like to fit a time series to only the part of the data that doesn't spike immediately upward.


title: Multiple Comparisons help

text: I have already asked a similar question but did not explain my problem in entirety. Any help would be greatly appreciated! I am analyzing some data to which I think I need to make some kind of multiple comparisons correction. The data comes from an economic game (Dictator game) and is the form of donations made by subjects ranging from $ 0.00 to $ 0.50 (can be any $ 0.01 denomination). Subjects were primed with an experimental image or an control image I am primarily interested in the effect of which image they were shown on their donation. I did the following tests all on the same donation data; 1) I used a Mann- Whitney U test to test if there was a difference in the amount donated for each image shown, control or experimental prime (Main hypotheseis). 2)I tested If there was a difference in the amount donated by each gender (Mann-Whitney U test) 3)I tested the relationship between age of dictator and donation using Spearmans Rank 4) Income and Education were categorical, eg. one level of Education would be "Some High School" and one level of Income would be "$0-$15 000" for each there was 6 levels. I tested donation against the different levels of Income (Kruscall-Wallis) and did the same for Education (Kruscall- Wallis). 5)Although dictators could choose any value to donate they all choose one of 11 amounts. For each of these amounts I tested if they were over represented by dictators who saw one image compared to the other. For each amount I encoded the data as 1 if they donated this amount and 0 for any other amount and ran a fishers exact test for each amount (11 tests in total). The only one that was significant was $0.50, those that saw one image gave $0.50 significantly more often than those that saw the other. I realize I have gone overboard with tests. I am not too sure too which I should apply an adjustment. ie. If I Bonferroni adjust for the kruscall-wallis test do I just divide alpha by 6 as this is how many levels in income and Education or do I need to take into account the Mann-Whitney tests and Spearmans rank. Also, for the Fishers exact tests the data is encoded differently for each test. Therefore do I need to make a correction or not? What sort of correction should I make? Thanks


title: Bandits with mixed reward processes?

text: I am trying to model a sequential exploration-exploitation problem with learning as a multi-armed bandit, where the reward mixes a Markovian and a stochastic reward. I understand how to model a bandit problem with Markovian rewards. In a Markovian bandit problem, each time the Decision Maker (DM) pulls an arm, the state of the arm changes according to a Markov chain and it obtains the reward associated with the transition. One way to solve this kind of problem is the Gittins index. On the other hand, I also know different ways to solve a bandit problem with stochastic rewards. In a stochastic reward bandit problem, each time the DM pulls an arm, it obtains a reward drawn from a known distribution, whose parameters are unknown and the DM tries to learn. There are several index policies and inference models to apply (UCB, Knowledge Gradient, etc.). But **what if** the reward mixes both concepts? That is, each time the DM pulls an arm, it obtains a reward which has a factor that is dependent on a state transition, but at the same time, other factor of the reward is drawn from a known distribution (with unknown parameters), common to all states of the Markov chain? To clear things up: **the reward is something like A_st*X, where "A_st" depends on the state** (to make things easier, we could consider it is deterministic given a state), **and "X" ~ N(mu, sigma), common to all states?** The only model I can come up with would rely on integrating my belief about "X" into the state space and then trying to compute the Gittins index through some sort of ADP (which may not be possible). Bonus round: is it possible to go even further and formulate this as a superprocess? (instead of a Markov chain, having an MDP?). **EDIT tl;dr:** How do you deal with a multi-armed bandit problem in which the reward is NOT in the form of the typical conjugate distributions (like the exponential-gamma learning model)?


title: How to find the probability that a random sample of 10 points out of 52 points of data will have a mean close to the total mean?

text: I have a data set that consists of 52 weeks of data. I want to find the probability that any 10 weeks of data within the set of 52 weeks is close to the mean of the total 52 weeks. Close meaning whatever measure seems appropriate (whether it is within a standard deviation, or what have you). I assume this is a fairly obvious problem, but I'm having trouble trying to figure out how you would get this.


title: Choosing post-hoc test after Kruskal-Wallis

text: I have a time series (five time steps) of samples from a population of the same ~60 individuals, each sample being a haphazardly chosen (i.e. not completely random) subset of the 60. The sample sizes each time ranges from about 15 to up to 50, and the measurement I do for each individual in a sample is a continous variable that goes from 0 to 1. Since the sample sizes are different each time, are NOT paired (the same individual may have been measured more than once, but that's not tracked), have unequal variances, and does not conform to a normal distribution, I opted to use the Kruskal-Wallis test in JMP version 9 as a non-parametric test for the existence of differences between samples. And there are differences (p < 0.01). As far as I can tell, JMP provides me with two possible post-hoc tests to find _where_ the differences are: "Steel-Dwass All Pairs", and "Dunn All Pairs for Joint Ranks". According to the JMP documentation, Steel-Dwass All Pairs... > performs the Steel-Dwass test on each pair. This is the nonparametric > version of the All Pairs, Tukey HSD And Dunn All Pairs for Joint Ranks... > performs a comparison of each pair, similar to the Steel-Dwass All Pairs > option. The Dunn method is different in that it computes ranks on all the > data, not just the pair being compared My questions are: (1) was I right in running the Kruskal-Wallis test to begin with? (2) What are the differences between the two post-hoc tests I am presented with, and how should I choose between them? Thank you so much for your help! Please let me know if there is any other information I should provide. P.S. I have attached the results from running both tests, do I inteprete them correctly that the cases with p < 0.05 are the statistically significant differences? (my alpha = 0.05) ![enter image description here](http://i.stack.imgur.com/Ui6wv.png)


title: How do I partition variance among nested random effects for non-normal data and an unbalanced design?

text: I have a dataset of plant drought tolerance values (called TLP_DRY) that I would like to partition variance for among the nested levels Biome/Study site/Species to figure out whether most variation in TLP_DRY is explained by biome, study site, or species identity. These levels are random effects. The things that make this complicated are that: 1) TLP_DRY is non-normal, even when log or square-root transformed, and is better fit with a Gamma distribution than a normal distribution (according to the fitdistrplus R package), and 2) different Biome levels have different numbers of Study Sites, and different Study Sites have different numbers of Species, making this an unbalanced design. It looks like I can't use the lme4 package in R because its an unbalanced design and it doesn't have a fixed effect. What kind of model is appropriate to partition variance for this kind of design? What R package would I use to do this? Thanks very much for your help!


title: Time series analysis: Determine if trend is deterministic fluctuating/stable or stochastic

text: I am analysing sales data of certain products and need to determine if the demand trend is deterministic fluctuating or deterministic stable or stochastic. How do I do that in R / what approach is usually applied to do so? Quarterly values would be for example: 2424,3030,3207,2943,2631,2863,2711,2837,2839,2876,3075,2886,2520,3136,3235,3071,2923,3237, 3054,3321,2856,3389,3092,3325


title: How can I calculate a confidence interval for the percentage change between two means?

text: The percentage change between two measurements is usually calculated as $100 \cdot \frac{x_2 - x_1}{x_1}$. If I have confidence intervals $(a_1, b_1)$ for $x_1$ and $(a_2, b_2)$ for $x_2$, is the confidence interval for the percentage change given as $100 \cdot \left ( \frac{a_2 - b_1}{b_1}, \frac{b_2 - a_1}{a_1} \right )$? (The thought being that $a_2-b_1$ is the smallest difference, while $b_2-a_1$ is the greatest difference.) Additional questions:   * Is this valid for confidence intervals for variables other than the (sample) mean?   * Is this valid for confidence intervals for sample means that were obtained for different sample sizes?   * Is this valid regardless of the distribution of the (sample) mean?   * Is this valid regardless of the procedure used to obtain the confidence interval (assuming that all the assumptions of the procedure were met)? It would be nice to have references backing up every answer.


title: Regarding the R packages that share the similar functionalities of NLTK toolkit

text: Are there any R packages that share the similar functionalities of NLTK toolkit?


title: How to compare yearly categorical data?

text: I have location data from 5 years (>30,000 points). Each location is given a classification name (in my case vegetation classes). These vegetation classes were assigned by intersecting the locations with vegetation maps representing 5 different years. The vegetation categories changed due to clear cut practices. Basically, some old forest areas in 2007 became clear cut in 2008, and others in 2009 became clear cuts in 2010 etc. Also, some old clear cuts became young forests because they grew back. What I want to know is, if I used just 1 map (say 2007) - how many wrong classifications will I have per year after this? Because what was once forest in 2007 could be clear cut in 2009. Knowing that the location from year 2009 _is_ a clear cut, using the map from 2007 - was it misclassified? I will want to test this for each year. Basically, is there a map that minimizes miss-classification of points. In the data set for each XY location, I have Year location was collected, clearcut07, clearcut08, clearcut09..., Names07, Names08, Names09.... The clearcut07 is binary indicator of 1 or 0, if was really a clear cut that year. the Names07 are the categorical names for all locations (all 5 years of points - if they were to be classified from that year). So, all locations obtained are given the attributes from the 2007 map, then the 2008 map then 2009, etc. I need a summary of counts of locations within each category. A count of locations that were properly classified as clear cuts from each year of the data. A count of locations from data collected, for example, in 2008 that were classified as clear cuts from each map? And a count of locations that were incorrectly classified as clear cuts from each map when I know the correct classification.


title: lmPerm p-values and multiple testing

text: I've started using lmPerm in order to perform regressions in R. The equation I want to fit has the form:               out3 <- lmp(outcome ~  bin1 + bin2 + cont1 + cont2, perm="Exact")      Where "outcome" is a non-normally distributed continuous variable, and `bin*` and `cont` are binary and continuous regressors (similarly, they are non- normally distributed). Each variable has a length of approx. 110 cases. Here are my questions:   1. This code works fine, but every time it runs in R, different p-values appear for each regressor. Which p-value should be reported in my results? I've tried repeating the test several times (in a loop) and getting an estimate from that, but I'm not sure it works...    2. If some of the predictors are changed and (then) several models / hypothesis are tested, should Bonferroni corrections be used in the same way they are applied for ordinary regressions? Is `lmp` somehow robust to multiple testing procedures?


title: How would you explain Markov Chain Monte Carlo (MCMC) to a layperson?

text: Maybe the concept, why it's used, and an example.


title: How do you decide the sample size when polling a large population?

text: Australia is currently having an election and understandably the media reports new political poll results daily. In a country of 22 million what percentage of the population would need to be sampled to get a statistically valid result? Is it possible that using too large a sample could affect the results, or does statistical validity monotonically increase with sample size?


title: What is the equivalent of a standard deviation when considering a least squares fit line?

text: I am calculating a tolerance interval following http://www.itl.nist.gov/div898/handbook/prc/section2/prc253.htm but this says to multiply the k value by the standard deviation of the sample. I have a model with a fit line, so I would think I do not want to use the standard deviation, but rather some value that reflects the residuals, and instead of using the sample mean, I will use the predicted value from my linear regression model. Is that right? What value do I use instead of the standard deviation? I could (maybe should have) asked the question this way: Given a linear model, how do I compute a one-sided tolerance interval. I think a tolerance interval is the right thing for my problem based on this: http://www.kmjn.org/notes/tolerance_intervals.html Edit again: I found this formula for "Assuming linear function and no replicates, the standard deviation about the regression" (from here) ![standard deviation about the regression formula](http://i.stack.imgur.com/8yH6g.gif) Is this the right fomula to get a value to multiply by the k values?


title: Topic models for short documents

text: Inspired by this question, I'm wondering whether any work has been done on topic models for large collections of extremely short texts. My intuition is that Twitter should be a natural inspiration for such models. However, from some limited experimentation, it looks like standard topic models (LDA, etc) perform quite poorly on this kind of data. Does anyone out there know of any work which has been done in this area? This paper talks about applying LDA to Twitter, but I'm really interested in whether there are other algorithms which perform better in the short-document context.


title: Books on information theoretic approaches in R

text: For all the numerous statistics books out there, surprisingly few seem to cover information theoretic approaches to model selection. Furthermore I cannot find ANY books which cover how implement information theoretic approaches in R. Does any know of a book which provides a good summary of information theoretic approaches using R?


title: Why use factor graph for Bayesian inference?

text: I don't understand why converting a Bayesian network into a factor graph is good for Bayesian inference? My questions are:   1. What is the benefit of using factor graph in Bayesian reasoning?   2. What would happen if we don't use it?  Any concrete examples will be appreciated!


title: Probability of winning money with a pair of loaded dice

text: You go gambling with a pair of loaded dice. Because of this, your odds of winning are 53% on every throw. Assuming the game pays 2:1 and you keep betting the same amount, how many games do you need to play to ensure an 80% likelihood of winning money? I am lost with where to start on this problem. I would appreciate some help so I can figure it out. Thanks.


title: Transform multiclass classification to binary - benefits?

text: I have 400 instances which must be categorized into 4 classes. Using WEKA, I tried out a couple of multiclass classifiers like J48 and Random Forests, but never made it above Kappa 0.6 and ~65% correctly classified instances (10-fold X-V) Then I thought about transforming the problem into a 1-vs-all classification, which usually yields accuracies of ~90%. I would then remove the one "single" class and keep the merged ones. Then, again, having only instances with 3 classes, I would perform 1-vs-2 and remove the instances classified as belonging to the single class, ending up with a binary classification problem. As I said - I always have like 90% correctly classified instances, but I fear that the 10% incorrectly classified instances add up and propagate through the splitting and dataset reduction process --- so in the end I would maybe end up with the same garbage output I'd have when performing the original multiclass classification?! What's the stand on this approach? Does it have any benefits at all?


title: What statistics to use for testing hypothesis about regression coefficient when n is greater than 30

text: What statistics to use for testing hypothesis about regression coefficient when n is 120. I plan using a z score but am not sure.


title: Cosine Distance and Recommender Systems

text: Suppose that a computer is rated according to three numerical features: processor speed, disk size, main-memory size. Consider three computers $A,B$ and $C$ with values: $$A(3.06,500,6)$$ $$B(2.68, 320,4)$$ $$C(2.92,640,6)$$ The cosine distance emphasizes differences in disk size more than differences in other features. What is the interpretation of the angle between the normalized vectors? If we normalize the vectors, than the cosine distance can take values between $0^{\circ}$ and $180^{\circ}$. If we do not normalize the vectors, then the cosine distance can take values between $0^{\circ}$ and $90^{\circ}$.


title: optimal choice of smooth.spline parameter?

text: I'm analyzing a time series (terms of trades) on which I want to perform a trend estimation by nonparametric methods like the above mentioned. By the way, I'm a total beginner with R and using the help files is already cryptic enough. Using the default settings of R gives me a smoother which simply follows the curve of the original series, although I thought that by using the cross validation method (which is set to default in R) optimal parameters would be used, i.e. which are a perfect trade-off between smoothness and best- fit. By comparison, I set df=8 which was a much nicer curve, however randomly choosing values isn't really scientific. Can you tell me what's the correct procedure here? ![http://i.stack.imgur.com/1ANGe.jpg](http://i.stack.imgur.com/hMn3U.jpg)


title: What is the difference between ordinal regression and ranking?

text: In both ordinal regression and ranking you are learning from an ordered dependent variables, so my question is: What is the difference in formulation (if any) between the ordinal regression problem and a learning to rank problem?


title: Leave-one-subject-out cross validation in Caret

text: **Hi Dear Colleagues,** I wonder how to correctly setup a leave-one-subject-out cross validation (LOSO) for train() function in caret. Here is my example code:               dat <- as.data.frame(cbind(rnorm(1:500,1),rnorm(1:500,10),rnorm(1:500,5),               rnorm(1:500,100),c(rep('1',100),rep('2',100),rep('3',100),rep('4',100),                                     rep('5',100)),rep(c(rep('X0',50),rep('X1',50)),5)))          colnames(dat) <- c('var1', 'var2', 'var3', 'var4', 'subject','class')      What would you then specify in trainControl() instead of cv (LOGCV?)? How would you modify the following code?               svmFit <- train(dat,y,                     method = "svmRadial",                     preProc = c("center", "scale"),                     tuneGrid = MySVMTuneGrid,                     trControl = trainControl(method = "cv", number = 10, classProbs =  TRUE))      Thank you very much for your time!


title: Model to detect abnormalities in two-value data?

text: I have a question concerning data that oscillates between two primary value ranges, and how one might go about determining any performance metrics or abnormalities in that data. Example: I have a text file of values that cluster at approximately 6000, and approximately 12000; the low value, and the high value respectively. I wish to determine some useful metric to perform on this data over a time range, to determine times when the data was "abnormal." How would one go about determining variations in this data, plus or minus a certain percentage? My current tactic is to take an average just so I can display some kind of metric, but clearly this is not the way to go for this data set.


title: I want to know how independent two variables are

text: I have no formal training in statistics so please correct me if I use the wrong terms as I try to explain my problem. I have a set of data (less than 80 points) with essentially 1 single outcome (a float we will call `dcl`) that can potentially depends on 10 of other variables, most of them floats maybe one or two boolean. While I might ask some multi-variate regression question later, let's start with something simple. Historically, people in my field have focused on the strongest correlation between `dcl` and variable `J` and some of my data suggests some other dependence on a float `Re` which is I'm sure at least weakly correlated with 'J' as they share some variables in their respective expressions. So my first question is how do I test the correlation and/or the independence of 'Re' and 'J' on the outcome 'dlc'? Intuitively and physically, I expect 'dlc' to depend strongly on 'J' and weakly on 'Re'. How do I prove this with a statistical analysis? Here are a few graphs to illustrate the data: ![Re vs J](http://i.stack.imgur.com/dVRMa.png) ![dcl vs J](http://i.stack.imgur.com/hNPvb.png) ![dcl vs Re](http://i.stack.imgur.com/D3YJe.png) Final point, in terms of software, I have python and R installed, I'm fairly proficient in python but I just installed R and know pretty much nothing about it. EDIT 1: Following gung's suggestion, I ran my dataset through R:               Call:     lm(formula = dcl ~ J + I(J^2) + Re + I(Re^2), data = df)          Residuals:         Min      1Q  Median      3Q     Max      -9.0078 -3.7930 -0.4458  2.0869 21.2538           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.648e+01  1.232e+00  13.380  < 2e-16 ***     J           -2.662e+00  3.773e-01  -7.054 6.58e-10 ***     I(J^2)       1.096e-01  2.071e-02   5.293 1.10e-06 ***     Re           1.966e-06  1.621e-05   0.121    0.904         I(Re^2)      2.191e-11  3.441e-11   0.637    0.526         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 5.369 on 77 degrees of freedom     Multiple R-squared: 0.4831, Adjusted R-squared: 0.4562      F-statistic: 17.99 on 4 and 77 DF,  p-value: 1.818e-10      So now I need some help to decipher this (but I will look into R documentation too). I don't know if it's relevant but on physical grounds only, it's probable the dependency in J is $dcl \sim \frac{1}{\sqrt{J}}$. Can I put this directly into the model? Does that already tell me something about the dependency on $J$ vs $Re$?               Call:     lm(formula = dcl ~ J + I(J^(-0.5)) + Re + I(Re^(-0.5)), data = df)          Residuals:         Min      1Q  Median      3Q     Max      -7.8119 -3.0097 -0.8504  1.8506 12.1439           Coefficients:                    Estimate Std. Error t value Pr(>|t|)         (Intercept)   8.175e-01  1.634e+00   0.500   0.6184         J            -2.946e-01  1.258e-01  -2.343   0.0217 *       I(J^(-0.5))   4.516e+00  7.053e-01   6.403 1.09e-08 ***     Re            3.332e-05  6.684e-06   4.985 3.72e-06 ***     I(Re^(-0.5))  6.009e+02  1.354e+02   4.438 2.98e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 4.426 on 77 degrees of freedom     Multiple R-squared: 0.6487, Adjusted R-squared: 0.6305      F-statistic: 35.55 on 4 and 77 DF,  p-value: < 2.2e-16          > model = lm(dcl ~ J+I(J^(-0.5)) + Re+I(Re^(-0.5)), data=df)     > summary(model)      **EDIT 2** : OK I'm starting to understand things better. Also, again based on physical grounds, I would think that the dependency is more something like $dcl ~ \frac{1}{\sqrt{J}} Re^{n}$ with possibly other variables in that product that I ignore. So when I enter such model in R (can I still use `lm` for something non-linear?), here is what I get:               Call:     lm(formula = dcl ~ I(J^(-0.5)) * I(Re^(-0.1)), data = df)          Residuals:          Min       1Q   Median       3Q      Max      -10.5363  -3.0192  -0.2556   1.4373  17.1494           Coefficients:                              Estimate Std. Error t value Pr(>|t|)         (Intercept)               -43.220      9.164  -4.716 1.03e-05 ***     I(J^(-0.5))                63.813     11.088   5.755 1.62e-07 ***     I(Re^(-0.1))              124.245     24.038   5.169 1.77e-06 ***     I(J^(-0.5)):I(Re^(-0.1)) -142.744     27.269  -5.235 1.36e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 4.668 on 78 degrees of freedom     Multiple R-squared: 0.6042, Adjusted R-squared: 0.5889      F-statistic: 39.68 on 3 and 78 DF,  p-value: 1.122e-15      Does that 4th line tell me something about the dependence between $J$ and $Re$? What kind of tools could I use to get an estimate on the exponent on Re? Because right now I'm just trying a few different numbers empirically to see how the errors evolve. Next for me is to plot the dcl vs the new model and see how well the data collapses visually... EDIT 3: In the end, I used `nls` to explore the possible exponents of my fit. I also removed some outliers in my data that used a different experimental method. I settled on a fit that gave me decent Pr(>|t|) and residuals and which visually produce a decent fit. The last outlier is actually another experiment with a different setup, but one that I trust. So in a sense it's good that it shows up as an outlier as it hints at other parameters that need to be explored. Thank you gung, I accept your answer as I believe it guided me in the right direction.               > model = nls(L.D ~ C*I(J^(c1))*I(Re_s^(c2)), start=list(C=10,c1=-0.25,c2=-0.25),data=df)     > summary(model)     Formula: L.D ~ C * I(J^(c1)) * I(Re_s^(c2))     Parameters:        Estimate Std. Error t value Pr(>|t|)         C  57.20389   26.40011   2.167   0.0337 *       c1 -0.27721    0.05901  -4.698 1.27e-05 ***     c2 -0.16424    0.04936  -3.327   0.0014 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Residual standard error: 4 on 70 degrees of freedom     Number of iterations to convergence: 8      Achieved convergence tolerance: 2.91e-06       ![enter image description here](http://i.stack.imgur.com/dfgWX.png)


title: Choosing c such that log(x + c) would remove skew from the population

text: I have data for which I would like to take the log transformation before doing OLS. The data include zeros. Thus, I want to do a log(x + c). I know a traditional c to choose is 1. I am wondering whether there is a way to have the data choose c such that there would be no skew anymore using features like the sample mean or variance? Is there a formaula?


title: How to measure uncertainty of a parameter when false positives exist?

text: The main goal of my research is to measure the percentage of brown dwarf stars in the Pleiades star forming cluster that are actually double stars (i.e. the brown dwarf star has a companion brown dwarf orbiting it). To do this, I need to estimate what is the faintest brown dwarf star detectable in our observations using the Hubble Space Telescope Wide Field Camera 3. My question is this: Say you measure a parameter $X$ 100 times ($X_1, X_2,...,X_{100}$). There are only true positive detections and false positive detections.   1. 80 of those measurements have a 99% chance of being a true positive.    2. 20 measurements have a 50% chance of being a true positive (50% chance of a false positive).  If there was no false positives, the uncertainty $\sigma_X$ = standard deviation of all measurements. How do you measure $\sigma_X$ when false positives exist as shown above? To elaborate: each measurement of parameter $(X_1,X_2,...,X_{100})$ has an assigned chance of being a false positive $(P_{f1},P_{f2},...,P_{f100})$. So how to measure $\sigma_X$ is complicated by the existence of these false positives.