title: Time Series Data with Varying Start and End Dates

text: I've done some research and it doesn't seem this question has been directly answered. I have a situation where the outcome variable is measured monthly from 1975-current and then multiple predictors with varying start dates and end dates. My question is what is the best approach to take with this data? I'm working with another statistician who has more experience with time series data in these situations and he tells me the common practice is to delete the data up to the shortest predictor. For example, if predictor 1 is measured from 1975-current and predictor 2 is measured from 2001-current, he tells me the best approach and only option is to use data for both predictors from 2001-current. So, essentially deleting all the data for predictor 1 from 1975-2000. He also said this often results in better prediction. I'm a bit skeptical of all this and was wondering if an ML or Bayesian approach in R or other common methods in R that might be able to handle this situation and use all available data. One thought is to handle this as a "variable occasion design" (see Chapter on Longitudinal Data in Snijders & Bosker, "Multilevel Analysis").


title: Use Bayesian hierarchical model to predict new data points

text: I have a data set $(n_i,y_i),i=0,...,10$. I modeled it as a Bayesian hierarchical beta-binomial model. $y_i∼Binomial(n_i,p_i)$ and $p_i∼Beta(\alpha,\beta)$. I have used MCMC to estimate $\alpha$ and $\beta$ (use the median as the estimated value). Given a new $(n_j, y_j)$, I want to judge if it comes from the same model as the training set $(i=0,...,10)$, namely $p_j$ is significantly greater than $\theta$ (shrinkage $p_i$ from all group). How can I use the estimated $\alpha$ and $\beta$ to compute it?


title: Test if two binomial distributions are statistically different from each other

text: I have three groups of data, each with a binomial distribution (i.e. each group has elements that are either success or failure). I do not have a predicted probability of success, but instead can only rely on the success rate of each as an approximation for the true success rate. I have only found this question, which is close but does not seem to exactly deal with the this scenario. To simplify down the test, let's just say that I have 2 groups (3 can be extended from this base case).   * Group 1 trials: $n_1$ = 2455   * Group 2 trials: $n_2$ = 2730 * * *   * Group 1 success: $k_1$ = 1556   * Group 2 success: $k_2$ = 1671 I don't have an expected success probability, only what I know from the samples. So my implied success rate for the two groups is:   * Group 1 success rate: $p_1$ = 1556/2455 = 63.4%    * Group 2 success rate: $p_2$ = 1671/2730 = 61.2% The success rate of each of the sample is fairly close. However my sample sizes are also quite large. If I check the CDF of the binomial distribution to see how different it is from the first (where I'm assuming the first is the null test) I get a very small probability that the second could be achieved. In Excel: > 1-BINOM.DIST(1556,2455,61.2%,TRUE) = 0.012 However, this does not take into account any variance of the first result, it just assumes the first result is the test probability. Is there a better way to test if these two samples of data are actually statistically different from one another?


title: Probability of observing a false correlation and confidence limits

text: In oil and gas exploration/development it is common to use acustic impedance derived from reflection seismic surveys to predict the porosity measured in wells drilled in the reservoir. I often use tables such as the one below (from a paper) to test for spuriousness of correlation:               0.87    0.78    0.72    0.67    0.63    0.57    0.49    0.39    0.32     0.75    0.58    0.47    0.40    0.34    0.25    0.16    0.09    0.05     0.62    0.40    0.28    0.20    0.15    0.08    0.03    0.00    0.00     0.50    0.25    0.14    0.08    0.05    0.02    0.00    0.00    0.00     0.39    0.14    0.06    0.02    0.01    0.00    0.00    0.00    0.00     0.28    0.07    0.02    0.01    0.00    0.00    0.00    0.00    0.00     0.19    0.02    0.00    0.00    0.00    0.00    0.00    0.00    0.00     0.10    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00     0.04    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00      where values in the table give the probability of observing the absolute value of the sample correlation coefficient, r, being greater than some constant R, given the true correlation (ρ) is zero, in other words the probability of a spurious correlation. These values are calculated with the expression (in both papers):               p=Pr(|r|≥R)=|t|≥((R√(n-2))/√(1-R^2 ))      where n is the sample size, or the number of locations (wells) where both reservoir property (porosity) and seismic attribute (acoustic impedance) are available, and t is distributed as a Student's t- critical value, with n-2 degrees of freedom. For the columns in this table n is, respectively:               5   10  15  20  25  35  50  75  100      For the rows, R (the magnitude of the spurious sample correlation) is, respectively:               0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9      This table is used to assess the chance that the sample correlation, r, is actually false or uncorrelated with the reservoir property being predicted. Quoting again:               For example, given 5 wells, and an r = 0.7, there is,a 19% probability that     the correlation is false.      I've used this method for years and (I thought) I understood well the theory and application. However I recently read this statement on a paper that expands on the original paper where the table was published: > ...however there is another aspect of the correlation coefficient that > should be considered — the confidence limits of the true correlation > coefficient. For this example, the 95% confidence limits are from a minimum > r of -0.48 (P97.5) and a maximum r of 0.98 (P2.5). Because the minimum r is > negative, we cannot say with confidence that there is any correlation and we > should reject this attribute as a predictor. Considering the same sample > correlation of 0.7, we need 9 samples before the minimum r is positive, but > its value is only 0.07, with a 4% chance that the correlation is false. This is my question: where is this coming from. Neither the original paper nor the recent one one published the data at each well, just the tables, so how can the author of the latter estimate the 95% confidence limits? All I could think of is bootstrapping the 95% confidence interval around the mean r ... except that even for that they'd need at least the one sample (the 7 wells) to get the mean. Is there any other way to get at that just using the values in the table? NB I originally posted the question in here: http://math.stackexchange.com/questions/912230/probability-of-observing-a- false-correlation-and-confidence-limits


title: Estimating probabilities using Bayes rule?

text: I am working on a past exam paper. I am given a data set as follows: Hair {brown, red} = {B,R}, Height {tall, short} = {T,S} and Country {UK, Italy} = {U,I} (B,T,U) (B,T,U) (B,T,I) (R,T,U) (R,T,U) (B,T,I) (R,T,U) {R,T,U) (B,T,I) (R,S,U) (R,S,U) (R,S,I) Question: Estimate the probabilities P(B,T|U), P(B|U), P(T|U), P(U) and P(I) As the question states estimate, I am guessing that I don't need to calculate any values. Is it just a case of adding up how many times P(B,T|U) occurs over the whole data set e.g. (2/12) = 16%. Then would the probability of P(U) be 0?


title: How do I test the difference between two dependent correlation coefficients whilst controlling for other variables?

text: How do I test the difference between r(X,Y) and r(X,Z) [Williams Test?] whilst controlling for demographic variables that were measured in addition to X and Y and Z [Partial Correlation]? Thanks, Moimeme


title: Does randomness imply independence?

text: By looking at a lag plot, it seems that my dataset is random (no autocorrelation whatsoever). Should this be sufficient to infer independence of observations? As a matter of fact, I _know_ that the observations are dependent and hence randomness should not guarantee that. But intuitively it feels otherwise. Given that randomness does not determine independence. Is there some statistical measure to test independence of observations? edit: Please note that the data is **not** of timeseries nature, it is the _average_ bandwidth usage of a customer across the month. ![enter image description here](http://i.stack.imgur.com/fuWId.png)


title: Pearson r with percentages across time

text: Data sources: CDC National Survey of Family Growth, CDC vital stats.   Population: American women aged 15-44 in 1982, 95, 2002, 2006   Sample sizes: several thousand Can I use Pearson r to show linear correlation between two percentage variables? That is, how % of American women who do X correlates with a percentage of events that are categorized a certain way? An example would be percent of Americans who report having driven drunk in the past month in various years compared to percent of motor vehicle collisions that were alcohol-related in same years. The sample is not the same over time but the population almost is, except that it's within an age range so for example the women who were 44 in 1982 will not be included in 1995 and so on. The reason I want to use percent is because the numbers from the 2 arms of the CDC don't exactly correspond. The National Survey of Family Growth asks women in the survey about births that they have experienced within the past 5 years and applies that to the whole American population, and that number is similar, but not identical, to the number of actual recorded births in the appropriate 5 year interval in the CDC Vital Statistics. I want to combine data from both sources into one aggregate category. Is that possible?


title: Can independent variables with low correlation with dependent variable be significant predictors?

text: I have eight independent variables and one dependent. I have run a correlation matrix, and 5 of them have a low correlation with the DV. I have then run a stepwise multiple regression to see whether any/all of the IVs can predict the DV. The regression showed that only two IVs can predict the DV (can only account for about 20% of the variance though), and SPSS removed the rest from the model. My supervisor reckons that I have not run the regression correctly, as due to the strength of the correlations, I should have found more predictors in the regression model. But the correlations were tiny, so my question is: if IVs and the DV hardly correlate, can IVs still be good predictors of the DV?


title: Add species characteristics, alongside count data, when using ordination

text: I have count data by site for 15 species, however I also have average size by species by site. I'm looking to use ordination techniques (non-metric multidimensional scaling) to reduce dimensionality and visualize differences in community composition. How can I incorporate the average size by species by site into the NMDS? I'm using the `metaMDS` function in `{vegan}` with R. Here is a simple representation of a similar type of data that I have               data(dune)     sizes <- matrix(runif(20*30), dim(dune)[1], dim(dune)[2]) #sizes corresponding to each species in each site          metaMDS(dune)      Can you somehow incorporate size data (or other attribute data for that matter) into an ordination technique like NMDS?


title: Multiple Logistic regression with binary random variables

text: I am attempting to do logistic regression on a data set of patients. The dependent variable is disease status (binary, yes or no). There are multiple independent variables some binary such as susceptibility gene (y/n) and others continuous such as patient age. When I have looked up logistic regression, in all the books I have read the examples given always have all the independent variables being continuous. Therefore is it mathematically valid for me to do logistic regression on my data set with continuous and binary random variables? Or is there an alternative regression model I should use?


title: How to do pre-election polling with statistics? (discrete results)

text: In statistics, it is taught about samples that can have some range of numbers as their values. However, in election like presidential election, there are only several candidates to choose. In that case, how can we assign assumed statistical distribution (such as normal distribution) to the poll and assign confidence interval?


title: Probability density and moment generating

text: ![enter image description here](http://i.stack.imgur.com/xsnKk.jpg) Please show some working when substituting the equation. Thanks


title: Mahalanobis distance for a multivariate normal distribution before and after uncorrelation

text: I have two questions: Suppose we uncorrelate variables of a multivariate normal distribution using Cholesky transformation. Then:   1. What is the relation between Mahalanobis distances before and after this transformation?    2. Can I use this method to calculate the Mahalanobis distance between the means of two multivariate normal distributions with different means and different covariance matrices? In other word, can I calculate the Mahalanobis distance by uncorrelating each distribution and then using Euclidean distance?


title: Fitting models in R with time restriction on coefficients

text: How should I define a model formula in "R", when one (or more) exact linear restrictions binding the coefficients is available. Equation: y = b1*x1 + b2*x1 where y = b1*x1 for t < t1 and y = b2*x1 for t > t1 Thank You..!


title: Estimating Poisson process intensity using GLM

text: Suppose I want to build an explanatory model for events generated by an inhomogeneous Poisson process with unknown intensity $\lambda$. Each entry in my dataset represents the registration of an event: time $t_i$, location $(x_i,y_i)$ and some explanatory variables $z_{i1},...,z_{ip}, i=1,...,n$. The entries are in chronological order. I want to estimate the effect of the explanantory variables on the intensity of the Poisson process using GLM. The entries cannot be grouped, because this inevitably leads to loss of explanatory information. Hence a count model with poisson regression does not work in this case. However, if we assume that the sojourn times $s_i$ (which we can compute in this dataset) are exponentially i.i.d., we can apply Gamma regression. Please ignore for the moment that we lose one observation due to differencing. Using a log link this becomes $E[s_i|\textbf{z}_i] = exp(\beta_0 + \beta_1 z_{i1} + ... + \beta_p z_{ip}) = \lambda(z_{i1},...,z_{ip})$ (Edit) Since the explanatory data analysis reveals that the process is inhomogeneous both in time and two-dimensional space, I also include the location coordinates in the linear component: $E[s_i|\textbf{z}_i,x_i,y_i] = exp(\beta_0 + \beta_1 z_{i1} + ... + \beta_p z_{ip} + \gamma_1 x_{i} + \gamma_2 y_{i}) = \lambda(z_{i1},...,z_{ip},x_i,y_i)$ Questions:   1. Is this approach valid?   2. Can I assume the sojourn times to be exponentially i.i.d. given that the Poisson process is inhomogeneous?   3. Does $E[s_i|\textbf{z}_i,x_i,y_i]$ indeed model the intensity of the process?   4. Am I modeling a one-dimensional (time), three-dimensional (time and location) or a $p+2$ dimensional Poisson process?


title: Plot directed acyclic graph with scaled edge length

text: I am trying to design a network (more precisely a directed acyclic graph) with specific edge lengths. The data is on the form of an edge list, and for each edge, there is an associated length. It might take the form of an edgelist matrix and a vector containing edge length, such as :               > edgelist = matrix( c("A","B","B","A","A","B","C","D","F","G"), ncol = 2 )     > length = c(1, 8, 7, 6, 5)     > edgelist          [,1] [,2]     [1,] "A"  "B"      [2,] "B"  "C"      [3,] "B"  "D"      [4,] "A"  "F"      [5,] "A"  "G"      > length     [1] 1 8 7 6 5      I can plot the network using               > network = as.network(edgelist, matrix.type = "edgelist", directed = TRUE)     > plot(network)      The network is directed and, it won't have directed cycles (it can have, if we ignore direction, some sort of simple cycles - they will be very rare in practice), and it shouldn't be impossible from a geometric point of view to scale the edge length to the value of the third column in 2D. It is to be noted that many nodes will have a degree > 3 . My question is : is there an algorithm, implemented in R (for the moment, I use the `network` package, but I am not limited to it and can use for instance `igrap` if necessary) that will allow me to draw such a graph ? I have looked into various functions (from network and other packages), but haven't found any. PS: This is a partly network analysis partly R-coding question, I hope it is on the right site here (rather than Stackoverflow).


title: R neural network model with target vector as output containing survival predictions

text: ## Overview I want to simulate the survival prediction using neural networks described in this paper entitled " _Application of Artificial Neural Network-Based Survival Analysis on Two Breast Cancer Datasets_ " by Chi, Street and Wolberg where a target vector of 1's representing time points are evaluated and returns a vector of the probabilities. Unfortunately they did not provide the software used to create the method. For example, a dataset with `time`, `status` (indicating alive or dead) and other predictive variables are needed to be trained in a neural network. There will be a data structure to hold the time points the research needs (say a `{1,1,1,1,1,1}`, six-month intervals in three years). The output of this model, after taking into account the aforementioned predictive variables will look like this: `{0.98761,0.91111,0.82710,0.70003,0.64253,0.47181}` that corresponds to:                         6      12      18      24      30      36     1   0.98761 0.91111 0.82710 0.70003 0.64253 0.47181     2   ...      These data can be used for future records that will be fed into the network. * * * ## R neural network packages I have been searching for ways to implement this in the `nnet`, `neuralnet` and `rminer` packages and unfortunately my limited knowledge can't modify them to suit my needs. I only know that these predict nominal and numeric values but now how to do vectors. So far the functions relating to predicting in those packages do not give a hint on the usage of a vector as output.   1. `nnet`'s `predict()`'s description clearly states below. > _Predict new examples by a trained neural net._              # X1, X2, and the rest of predictive variables     model.nnet <- nnet(Surv(time,status)~X1+X2, data=data.train, size=1, maxit=500)        2. `neuralnet`'s `compute()` does not support a similar target vector for support. > _Computation of a given neural network for given covariate vectors._   3. `rminer`'s `predict()` and `lforecast()` prove to have potential (based from the examples) but I have no idea how to transform them into what I want to do.   4. the `survnnet` package is said to support the common `predictSurvProb` function of the `pec` package but I stayed away from it because of the poor documentation and support here in the Internet. > _Predict new examples by a trained survival neural net._              model.survnet <- survnnet(Surv(time,status)~X1+X2, data=dat,         model='llog', decay=0.1, bias.decay=25, size=1,         skip=T, alpha=0.1)       predictions   <- predict(model.survnet, data.train, type="raw")      * * * ## Cox PH (current known method) The closest method I've got so far was `coxph` and `cph` applied with the `predictSurvProb` function where a `times` variable is declared with the numeric points of interest. Usually I do this:               data.train  <- SimSurv(300)     model.coxph <- cph(Surv(time,status)~X2,data=dat,surv=TRUE,x=TRUE,y=TRUE)          # declare target times as 25,50,75,100,150 for probabilities     predictions <- predictSurvProb(coxph12, newdata=data.train, times=c(25,50,75,100,150))     round(predictions, digits=6)       and I would get an output of probabilities per time period I specified:                         25       50       75      100      150     1   0.648268 0.509353 0.460196 0.425917 0.324364     2   0.648268 0.509353 0.460196 0.425917 0.324364     3   0.756732 0.648020 0.607077 0.577596 0.484789     4   0.648268 0.509353 0.460196 0.425917 0.324364     5   0.648268 0.509353 0.460196 0.425917 0.324364      These are now ready to be integrated as new variables into the dataset for other purposes. * * * I now want to implement this in a neural network with the target vector or a similar implementation to the Cox PH process like the example above. Unfortunately I can not find a straightforward package or tutorial online that states if this is possible to do, as I said before. I checked into the data mining software WEKA's `MultilayerPerceptron` implementation but it requires that the outcome variable (or in my case a vector of probabilities) to be existing first, which takes me back here into R. **Is this method possible in these R packages or the only way is to create my own? All help and suggestions would be greatly appreciated.**


title: F-statistic value is too high for my model

text: > Number of obs = 501 , Method: Fixed-effects regression   >  Number of groups = 101 > >> F( 10, 100) = 3422.31 >> >>> , Prob > F = 0.0000 , within R-squared = 0.6325 My model (panel data with fixed effects) has six variables and four dummies for five years and as you see above, the F-value is too high with a value of 3422. Is this normal? Or what could be the problem? I should also mention that when I add one variable to the six, the f stat turn from 3422 to 86. Thank you in advance


title: How to report results from Friedman's test including post-hoc results

text: I have performed comparison of repeated measures for a group using Friedman’s test (K-related-samples test) and additional post-hoc analysis with Wilcoxon signed-rank test and Bonferroni correction. I have been told that this would be correct to do when looking for a change between the measurements and your measures are not normally distributed. **Now I'm not sure how to report my results correctly. Is it correct to a report mean with standard deviation or median with 25th to 75th percentile?** Example with mean: Comparison of the repeated measures was performed using Friedman’s test showing a statistically significant decrease in concentration, χ2(3) = 24.877, p < 0.001. Post-hoc analysis with Wilcoxon signed-rank test was conducted with a Bonferroni correction applied, resulting in a significance level set at p < 0.0083. The mean concentration (±SD) was 0.045 (0.012) at baseline, 0.032 (0.012) at time2, 0.027 (0.007) at time 3 and 0.024 (0.006) at time4. A significant decrease was seen between time2 and baseline (Z = 2.900, p = 0.004), between time3 and baseline (Z = 3.180, p = 0.001) and between time4 and baseline (Z = 3.180, p = 0.001). Could that work? Thanks for helping! Best, Rebecka.


title: The probability of getting one variation of consecutive outcomes in Bernoulli trials before another

text: Consider tossing a fair coin and writing down the outcomes in a sequence of symbols. "H" stands for head and "T" for tail. Let A , B and C be the words HTHH , HHTH and THHH respectively. What is the probability of encountering each A, B and C as a subword in your sequence of outcomes before both others? **Edit:** My question is relateted to this question. While its answers intuitively explain why some patterns hit sooner than others, none of them include a way of calculating those probabilities. Which is what I'm asking. **Edit 2:** I have already recieved an answer. Having read the wikipedia articles, this question still stumps me. I'm starting a bounty on it.


title: Contradictory xtunitroot results

text: I am using unit root tests in Stata's xtunitroot command. The data are annual wage replacement rates for unemployment benefits in each of the 50 states for 1977 to 2006. (They vary noticably from state to state and year to year) I am getting contradictory results for tests of a common unit root. In all cases the tests specify both demean and detrend. the Leven Lui Chin result indicates that I can reject the null of "panels contain unit roots. p<.001 The HT version indicates I cannot reject the null "panels contain unit roots" (p<1.00) The Breitung version tells me that I cannot (quite) reject the null "panel contains unit roots (p<.08) Any advice?


title: Is there a formula to determine the cumulative distribution associated with multiple different kinds of events?

text: Similar to this question, but with more than one different sided die. Is there a way to calculate probability distributions associated with dice, where the probability distributions for each die are not all the same? **Example:** Using dice, find the probability of rolling 1 once in a set when rolling 2d4 and 1d6. As I understand it, the long way is to calculate each event separately and then add up the events where a 1 is rolled once: ` # `` result ``#of 1s` `1: ``|{1,1,1}``   3  `   `2: ``|{1,1,2}``   2  `   `3: ``|{1,1,3}``   2  `   `4: ``|{1,1,4}``   2  `   `5: ``|{1,1,5}``   2  `   `6: ``|{1,1,6}``   2  `   `7: ``|{1,2,1}``   2  `   `8: ``|{1,2,2}``   1  `   `9: ``|{1,2,3}``   1  `   `10:``|{1,2,4}``   1  `   `11:``|{1,2,5}``   1  `   `12:``|{1,2,6}``   1  `   `13:``|{1,3,1}``   2  `   `14:``|{1,3,2}``   1  `   `15:``|{1,3,3}``   1  `   `16:``|{1,3,4}``   1  `   `17:``|{1,3,5}``   1  `   `18:``|{1,3,6}``   1  `   `19:``|{1,4,1}``   2  `   `20:``|{1,4,2}``   1  `   `21:``|{1,4,3}``   1  `   `22:``|{1,4,4}``   1  `   `23:``|{1,4,5}``   1  `   `24:``|{1,4,6}``   1  `   `25:``|{2,1,1}``   2  `   `26:``|{2,1,2}``   1  `   `27:``|{2,1,3}``   1  `   `28:``|{2,1,4}``   1  `   `29:``|{2,1,5}``   1  `   `30:``|{2,1,6}``   1  `   `31:``|{2,2,1}``   1  `   `32:``|{2,2,2}``   0  `   `33:``|{2,2,3}``   0  `   `34:``|{2,2,4}``   0  `   `35:``|{2,2,5}``   0  `   `36:``|{2,2,6}``   0  `   `37:``|{2,3,1}``   1  `   `38:``|{2,3,2}``   0  `   `39:``|{2,3,3}``   0  `   `40:``|{2,3,4}``   0  `   `41:``|{2,3,5}``   0  `   `42:``|{2,3,6}``   0  `   `43:``|{2,4,1}``   1  `   `44:``|{2,4,2}``   0  `   `45:``|{2,4,3}``   0  `   `46:``|{2,4,4}``   0  `   `47:``|{2,4,5}``   0  `   `48:``|{2,4,6}``   0  `   `49:``|{3,1,1}``   2  `   `50:``|{3,1,2}``   1  `   `51:``|{3,1,3}``   1  `   `52:``|{3,1,4}``   1  `   `53:``|{3,1,5}``   1  `   `54:``|{3,1,6}``   1  `   `55:``|{3,2,1}``   1  `   `56:``|{3,2,2}``   0  `   `57:``|{3,2,3}``   0  `   `58:``|{3,2,4}``   0  `   `59:``|{3,2,5}``   0  `   `60:``|{3,2,6}``   0  `   `61:``|{3,3,1}``   1  `   `62:``|{3,3,2}``   0  `   `63:``|{3,3,3}``   0  `   `64:``|{3,3,4}``   0  `   `65:``|{3,3,5}``   0  `   `66:``|{3,3,6}``   0  `   `67:``|{3,4,1}``   1  `   `68:``|{3,4,2}``   0  `   `69:``|{3,4,3}``   0  `   `70:``|{3,4,4}``   0  `   `71:``|{3,4,5}``   0  `   `72:``|{3,4,6}``   0  `   `73:``|{4,1,1}``   2  `   `74:``|{4,1,2}``   1  `   `75:``|{4,1,3}``   1  `   `76:``|{4,1,4}``   1  `   `77:``|{4,1,5}``   1  `   `78:``|{4,1,6}``   1  `   `79:``|{4,2,1}``   1  `   `80:``|{4,2,2}``   0  `   `81:``|{4,2,3}``   0  `   `82:``|{4,2,4}``   0  `   `83:``|{4,2,5}``   0  `   `84:``|{4,2,6}``   0  `   `85:``|{4,3,1}``   1  `   `86:``|{4,3,2}``   0  `   `87:``|{4,3,3}``   0  `   `88:``|{4,3,4}``   0  `   `89:``|{4,3,5}``   0  `   `90:``|{4,3,6}``   0  `   `91:``|{4,4,1}``   1  `   `92:``|{4,4,2}``   0  `   `93:``|{4,4,3}``   0  `   `94:``|{4,4,4}``   0  `   `95:``|{4,4,5}``   0  `   `96:``|{4,4,6}``   0  `   Which totals 39 out of 96 possibilities that a 1 will be rolled only once, or 40.625% But I know there has to be formula that does this because if I can tediously count it out manually (as above), then there must be some solid math function behind it. Through the past few days of research, I've learned that this can also be done using the cumulative distribution function. However, so far as I understand it, the cdf requires you to have the same probability for all three rolls - so it would work if I did one formula for 2d4 and then one for 1d6 but not cumulatively. Delving deeper, it seems that the Poisson distribution may provide the probability distributions that I want, but honestly, my capacity for learning this from a formula on wikipedia is rapidly approaching zero. If this is the correct formula (or if the answer lies with a different one), could you provide a proof using my example from above? Similarly, could you provide an example where there are several different probabilities (ex: 1d2, 2d4, 1d6, 3d10)?


title: Calculating ICC for random-effects logistic regression

text: I'm running a logistic regression model in the form:               lmer(response~1+(1|site), family=binomial, REML = FALSE)      Normally I would calculate the ICC from the intercept and residual variances, but the summary of the model doesn't include residual variance. How do I calculate this?


title: Does unequal sample size matter for single response

text: I am doing a corpus analysis on the relationship between the direction of signs and their valence ratings (from another corpus). I coded each sign as either upward, downward or non-vertical, and has one valence rating. I performed a one-way ANOVA with to determine if direction predicts the valence ratings.   1. Should I be worried that there are unequal sample sizes (e.g., 4 times more downward than upward signs)? I was told that unequal sample sizes is not a problem because I am not aggregating the valence ratings in each direction (e.g., not comparing sample means). Is this right and why? Thank you.


title: fit and cross-validate categorical sample data formed from observations

text: I am working with the following sample dataset:               t1  t2  ntrial  nsuccess      1       4    1000       4      1       8    1000       8      2       4    1000       4      2       8    1000       8      Predictors (features) t1 and t2 are categorical: t1 has categories 1 and 2, and t2 has categories 4 and 8. There is no ordering to these categories. In reality, for each combination of t1 and t2, I observe either a success, or a failure. For example, in first row (t1, t2) = (1, 4) and I recorded 4 successes from 1000 observations. So really the dataset has 4000 rows (call this unrolled binary data), which can be compressed for the purposes of, say, Logistic Regression. Obviously unrolling this data is memory-inefficient, because unrolled data will look something like this:               t1  t2  success      1  4    1      1  4    1      1  4    1      1  4    1      1  4    0       ... 996 zeros      1  4    0      1  8    1      1  8    1      ... and so on      In this toy example, I only have 1000 observations per row, so the data can be unrolled, but in reality my dataset has 10^9 observations, and many features/predictors, so unrolling it is not really feasible. I would like to know:   1. Is there an algorithm (in say, R, or Java, etc) which can operate on this dataset directly, without me having to unroll the rows? The algorithm has to be relatively fast and simple to train and predict.   2. Is there a cross-validation routine which can also work on this dataset without unrolling the rows? Here's some starter R code which I tried so far. I fit weighted Logistic Regression (the more observations we have, the lower the variance of predicted quantity, so we have weighted regression on our hands) and Regression Tree.               library(boot)     set.seed(1)     input_file = 'data\\test\\test.txt'     # number of cross-validation folds     K = 10     input <-        read.csv(input_file, header=TRUE, sep='\t', quote="",                 colClasses=c(rep('factor', 2), 'numeric', 'numeric'))     # change the contents of the frame     input$nfail = pmax(0, input$ntrial - input$nsuccess)         # compute the probability of success         input$prob = input$nsuccess/input$ntrial          # fit the main model     glm.model =         glm(cbind(input$nsuccess, input$nfail) ~              input$t1 + input$t2, family = binomial, weights = input$ntrial)          cost <-        function(y,yhat) sum(input$ntrial*((y-yhat)^2))/sum(input$ntrial)     # this produces an error, cv.glm doesn't realize rows have to be unrolled     cv.glm(input, glm.model, cost = cost, K=K)$delta[1]          library(tree)     control = tree.control(sum(input$ntrial))         tree.model <-            tree(input$prob ~ input$t1 + input$t2,             weights = input$ntrial, control = control)     plot(tree.model)     text(tree.model, pretty=0)     # this is too big to fit in memory, R fails     cv.tree(tree.model, K = K)      For example, if I want to perform 10-fold cross-validation, I have to sample each fold randomly from each of 4000 rows, either sampling a 1 or a 0 (success or failure), and then for each fold I add up the numbers of successes and failures for each particular combination of t1 and t2.


title: How to specify a Bayesian binomial model with shrinkage to the population?

text: # Problem I’m currently working on a problem where I have count data for $n$ items in the following form:                   Item   Count   Total         --------------------         1          1      30         2         10     340         3          0       0         4          1      10                  ...      I want to calculate the posterior of the success probability $p_i$ for each item. Additionally, I suspect that all the items have a similar or even the same success probability but I want the model to determine that. My expectations for the resulting model are:   * Estimate the underlying global parameter of all the items _and_   * Shrink the success probabilities $p_i$ towards this global parameter in order to deal with items which have few observations.   * Items with observations indicating a different success probability should be able to “override” the global success probability. # Stan code This is the Stan code I implemented so far:                   data {           int<lower=0> J; // amount of items           int count[J];           int total[J];         }         parameters {           real<lower=0,upper=1> mu;           real<lower=0.001> ss;           real<lower=0,upper=1> prob[J];         }         model {           mu ~ uniform(0, 1);           ss ~ gamma(0.01, 0.01);           prob ~ beta(mu*ss, (1-mu)*ss);           count ~ binomial(total, prob);         }      Here I achieve the shrinkage by assuming that the success probabilities of the items are sampled from a beta distribution (parametrized using mean and sample size) where the mean is the global success probability. The problem occurs when specifying the prior of the sample size. If I choose a very vague prior, Stan requires a lot of iterations until the effective sample size is large enough. Also I’m not sure if this is the right approach since it is not clear if the amount of shrinkage this model provides is always appropriate. I’m interested if there is an elegant way to implement my expected behavior or if there is some flaw in my overall model.


title: Time interval as the independent variable

text: ![time interval as the independent variable for linear regression](http://i.stack.imgur.com/8jFuh.png) How do I find the $a$ and $b$ in $Y=a+bX$? The way I have been doing it so far are by the formulas with tons of $\Sigma(X)$ and $\Sigma(Y)$


title: Missing values on weekends in transportation forecasting

text: I am currently analyzing transport volume time series (private consumption goods) on a daily basis. I am trying to fit an auto regressive model an than also include exogenous factors (as weather, web search activity etc.). There are two obvious seasonal patterns. During the week and during the year. However, all values at Saturday/Sunday are zero (trucks are not driving). Just ignoring weekends is not feasible, because the exogenous variables in the model occurring at the weekend influence transportation volumes during the week, e.g. web search activity on saturday influences transport volume on monday. How can this missing values be handled appropriately?


title: Can we determine the dimensions left after implementing PCA?

text: I have a data set in which I need to reduce the dimensions. Each dimension is labelled say $x_1, x_2, x_3, \cdots, x_n$. Now if I implement Principal Components Analysis (PCA) on it, I get some $p$ dimensions where $p$ is less than $n$. Now how do I know which $p$ dimensions are these among the initial $n$ dimensions. Can we even determine that?


title: How to interpret variation explained by principal coordinates?

text: I have recently seen a couple of Principal Coordinates Analysis (PCoA) projection plots which show "percentage variation explained" by the respective principal coordinates. Given that the analysis is not done on a co-variance matrix (it is usually done on some ecological metric) it seems wrong to imply that the magnitude of the eigenvalues is in any way related to the variance explained. My question is, what would the best interpretation of the eigenvalues be?


title: What is the 'same distribution' mean?

text: Say if I have two random variable X and Y and they have the same distribution, what is that suppose to mean? Is that mean they have same mean and variance?


title: multi-armed bandit with seasonality

text: I'm working on implementing a multi-armed-bandit-like approach for determining the best price to offer for a product. Our goal is to optimize profit, meaning, we want to find the price where (price-cost) x number of purchases is greatest. The problem is that there is a known pattern of seasonality within the week. So, if I set the price to the same thing on a Monday vs. a Wednesday, they might have dramatically different volume. That's my specific situation, but it started me wondering more broadly about using a multi-armed bandit when there are external factors (aside from which "arm" of the bandit you pull) that affects the results you get. Is there a way of taking these factors into account? Or-- perhaps there are results that suggest that if you run the algorithm long enough, it converges in spite of the other factors?


title: Nested mixed model

text: I have a data set where we have collected fish for behavioural experiments from different areas and within each area we have selected a clear water site and a eutrophied site. We then brought the fish to lab - all fish were used only once, they built a nest and that was it. What we wanted to know was if habitat type (clear or eutrophied) affected the nest building. This is the original model (for example for nest weight, male weight as a covariate):               m1 <- lme(NestW ~ Habitat + MaleW, random = (~1|Area/Site))      The trouble is `dfs` \- they are way too high and reflecting the amount of fish used, not the areas, which we would want to be the unit of replication here. So, how to fix this? (adding Habitat to the random effect does not help)


title: What is the adjusted R-squared formula in lm in R and how should it be interpreted?

text: **What is the exact formula used in R lm() for the Adjusted R-squared? How can I interpret it?** ### Adjusted r-squared formulas There seem to exist several formula's to calculate Adjusted R-squared.   * Wherry’s formula `[1-(1-R2)·(n-1)/(n-v)]`   * McNemar’s formula `[1-(1-R2)·(n-1)/(n-v-1)]`   * Lord’s formula `[1-(1-R2)(n+v-1)/(n-v-1)]`   * Stein `[1-(n-1/n-k-1)(n-2)/n-k-2) (n+1/n)](1-R-squared)` ### Textbook descriptions   * According to the textbook Field, Discovering Statistics Using R (2012, p. 273) R uses Wherry's equation which _"tells us how much variance in Y would be accounted for if the model had been derived from th. population from which the sample was taken"_. He does not give the formula for Wherry. He recommends using Stein's formula (by hand) to check how well the model cross-validates.   * Kleiber/Zeileis, Applied Econometrics with R (2008,p. 59) claim it's _"Theil's adjusted R-squared"_ and don't say exactly how its interpretation varies from the multiple R-squared.   * Dalgaard, Introductory Statistics with R (2008, p.113) writes that _"if you multiply [adjusted R-squared] by 100%, it can be interpreted as '% variance reduction'"_. He does not say to which formula this corresponds. I had previously thought, and read widely, that R-squared penalizes for adding additional variables to the model. Now the use of these different formulas seems to call for different interpretations. I also looked at this question on differences between r-square and adjusted r-square, and this site ### Questions   * **Which formula is used for adjusted r-square by R lm()?**   * **How can I interpret it?**


title: How are the classes determined in a growth mixture model?

text: I am trying to fit a growth mixture model to a set of trajectory data. I am using a model that on a generic level is the same as this one. I do understand how the latent intercept, linear, and quadratic variables are determined. However, I do not understand how the latent class variable is determined. Do I have to assign cases manually to different classes (as if I were doing a multi-group analysis when doing Structural Equation Modeling) or is there some automatic way of assigning cases to the classes? If so, how does this work?


title: Choosing which variable to sample to get a better model

text: Apologies if this question is long-winded and vague, but I'm really not familiar with the field and I'm having a hard time finding references. We have a random variable $X$ which we assume is governed by a distribution with unknown parameters $\theta^*$ in some parameter space $\Theta$. We will write $p_\theta(x)$ the probability of outcome $x$ in the distribution with parameters $\theta$. We have a _loss function_ $L : \Theta \times \Theta \rightarrow \mathbb{R}$ such that $L(\theta^*, \tilde\theta)$ is the cost of having predicted $\tilde\theta$ as the underlying model of $X$ whereas the true model was $\theta^*$. We assume that we have a collections of independent samples $S$ drawn according to $\theta^*$ by sampling $X$. Those samples define a _likelihood function_ $\mathcal{L}_{S} : \Theta \rightarrow [0, 1]$: $\mathcal{L}_{S}(\theta)$ is the probability of the observed $S$ for parameter $\theta$. From there, I would like to define the _risk (or expected loss) of a prediction $\tilde\theta$_ (relative to the probability distribution inferred from $S$) as $e_{S}(\tilde\theta) = \int_{\theta \in \Theta} \mathcal{L}_{S}(\theta) \cdot L(\theta, \tilde\theta)$. Since the likelihood function is not a pdf over $\Theta$, this integral may not even be well- defined, so by this point I'm already confused. The _best prediction_ is $\widehat{\theta_{S}} = \textrm{argmin}_\theta e_{S}(\theta)$ and the _best risk_ is the risk of the best prediction $\widehat{e_{S}} = e_{S}(\widehat{\theta_{S}})$. I now look at what would happen if I am allowed to draw a new sample for $X$. The _next sample distribution_ is $P_{S}(x) = \int_{\theta \in \Theta} \mathcal{L}_{S}(\theta) p_\theta(x)$, our current estimate of the probability distribution on $X$ according to $S$. The _next best risk_ is $\int_{\theta \in \Theta} \mathcal{L}_{S}(\theta) \int_x p_\theta(x) \widehat{e_{S \cup \\{x\\}}}$: it is the best risk when drawing one more sample from $X$ and adding it to $S$, where we estimate the distribution of samples from $X$ using the next sample distribution (inferred from the currently observed samples). Now, my situation is that I have different variables, each of them having an unknown parameter, I have independent samples for them, and I have to choose for which one of them I want to request a new sample. The goal is to choose the variable for which a new sample is expected to achieve the best risk decrease. My questions are:   1. Is there a way to ensure that my definition of risk (and the subsequent definition) make sense?   2. Is this approach known? What are the standard names of the concepts that I define? Are people using methods of this kind in practice?   3. Do some of the above expressions simplify for some reasonable choices of the probability distributions (e.g. normal distributions) and the loss function (e.g. the KL-divergence)? For example, can it be the case that the next sample distribution (as I defined it above) is actially the distribution given by the maximum likelihood estimation $\theta = \mathrm{argmax}_{\theta \in \Theta} \mathcal{L}_S(\theta)$? Can it be the case that the best prediction is actually the maximum likelihood estimation?   4. Can we prove that the next best risk is always less than the current best risk under reasonable hypotheses on the probability distribution and the loss function?


title: Approximation of Pr(X > a), with X a multivariate normal rv

text: Let $X = (X_1, ..., X_p)$ a random variable with a $N(\mu, \Sigma)$ distribution. $$ $$ $$ \Pr(X_1 > a_1, ..., X_p > a_p) \\\ =\int_{a_1}^\infty ... \int_{a_p}^\infty (2\pi)^{-p/2} (\det(\Sigma))^{-1/2} \exp\left\\{-\frac{1}{2} (x-\mu)' \Sigma^{-1} (x-\mu)\right\\} dx_1 ... dx_p $$ $$ $$ Does it exist an (analytical?) approximation of this multivariate integral?


title: Characterizing relationship between two datasets

text: I have two data sets of electricity prices in a given region, one with data at 5 minute intervals and one with just the hourly values (which are almost, but not quite, averages of the 5 minute values within the hour). The hourly average data set goes back much further in time. I would like to study the variances of the 5-minute data vs. the hourly so I can make estimations or predictions as to the chance of the 5-minute price exceeding a given value if I only have the hourly average. How would you go about this? My Stats 101 approach has been to look at the standard deviation of either the ratio or the difference between the 5 minute and hourly data. But, now I'm realizing that that was assuming a normal distribution for the ratio/differences, which turns out not to be true. Given a large data set of mappings of 1 hourly value to 12 5-minute values, how can I estimate the chance that any hourly value maps to one or more 5-minute values >= some threshold?


title: Multiplicative error in survey data

text: I'm working on a panel survey data where each individual's income was multiplied by a individual-specific random number (each random number is evenly distributed from 0.5 to 1.5) to avoid any participant being recognized from the data. Basically, I cannot directly compare two individual's income using the data. For this type of multiplicative error, I was wondering do I need to use any particular method to correct it? Any advice will be appreciated! Problem solved. Please see Nestor's answer.


title: Statistics and R: Determine if N-dimensional vectors are all similar

text: I need to see if all 35 vectors of 120 dimensions are similar, to prove the hypothesis that they are not. I used euclidean distance and cosine similarity and by learning through this forum I got a similarity matrix. (its a lot of data so I will summarize this point in R code)               matrix<-cbind( 35 vectors) #matrix is then 120x35     sim1<-cosine(matrix)     dist1<-as.matrix(dist(t(matrix), diag= TRUE, method="euclidean"))      But, the thing is that cosine compares vector_i with vector_j, and I would prefer to see vector_i with the rest of the vectors and evaluate a p-value of some sort. For me, it would be like:               similarity(vector_i, other vectors)>0.90 --> vector_i is similar to the other vectors      With similarity hopefully being like cosine similarity, as it is more about direction than magnitude. Does anyone have a recommendation or know a similarity criteria for these problems? I am starting to use R, so if you know a test that is implemented there, i would be very grateful. PD1>Statistics was a class I took a loong time ago, sorry if there are misunderstandings) PD2>the value of each dimension can go between 0 and 1, and many times it is 0, so I think that chi squared might not be ok for this...


title: Should the Shapiro-Wilk test and QQ-Plot always be combined?

text: I have read multiple places that Shapiro-Wilk test should always be added with a QQ-plot, but no one has given a reason, and I do not see the intuition behind this. Can anyone explain why one need to confirm a Shapiro-Wilk test with QQ-plot?


title: Is It Correct to Talk About the Probability of an Event that Has Already Occurred and Is Not Repeatable?

text: I have a limited understanding of probability, PDFs, etc. When the Viking landers went to Mars, there was talk of the probability of life on Mars. Now we read of the probability of life beyond the solar system. It seems to me though, that there either is life beyond the solar system or there isn't. It's not like putting one's hand into a jar filled with black and white marbles in some proportion and taking one out. One might repeat that experiment many times and try to estimate the true proportion of black vs. white, or one might know the proportion beforehand and calculate the probability of pulling out 10 black and 5 white, or 15 all white, etc. I just have a hard time understanding what the "probability of life elsewhere in the universe" actually means.


title: Correlation between two variables measured in separate experiments in R

text: I'm trying to show a correlation between growth in a petri dish of some fungi and its effect on a plant. I have ten strains of fungi which I tested in the plant and in petri dishes. I can put data from both experiments into linear models (lm) to get estimated means and variances. If I run a regression on the estimated means for each strain I find a significant correlation, but this doesn't take uncertainty in the strain means into account. So far I've tried a parametric bootstrap, but I'm having a hard time figuring out the ultimate test statistic. I've included the parameter estimates and my R code. At the moment it generates simulated strain means from the estimated parameters. I need a p value for the covariance between growth and virulence. Any help would be greatly appreciated! My specific questions are: 1) Is a parametric boostrap the right way to go about analyzing this? 2) How would one go about doing this in R?               Strain  MeanVirulence MeanGrowth    VirulenceVarGrowthVar       1  -5.26064 0.066716    0.67834 0.053247       2   -4.05482    -0.055524   0.68385 0.047111       3   -5.47282    0.029047    0.68385 0.046739       4   -3.50632    -0.161811   0.68385 0.047083       5   -4.94051    -0.224949   0.68385 0.04727       6   -4.04982    -0.062938   0.68385 0.047647       7   -4.53178    -0.142985   0.68385 0.04788       8   -3.01697    -0.199349   0.68385 0.047255       9   -3.81093    -0.254793   0.68385 0.047255       10  -1.61882    -0.325289   0.68385 0.0469      And here is my R code:               gendata<-function(par,npar=TRUE,print=TRUE){         n     = 10         k     = 2         x=matrix(data=NA, nrow=n, ncol=k)         for(i in 1:n){           x[i,1] = rnorm(1,mean=par[i,2],sd=par[i,4])           x[i,2] = rnorm(1,mean=par[i,3],sd=par[i,5])         }         return(x)     }          lmp <- function (modelobject) {         f <- summary(modelobject)$fstatistic         p <- pf(f[1],f[2],f[3],lower.tail=F)         attributes(p) <- NULL         return(p)     }               samp=20000     rescor=matrix(data=NA, samp)     resvar=matrix(data=NA, samp)     pvals=matrix(data=NA, samp)     numsig = 0     numnotsig = 0     for (i in 1:samp){         x<-gendata(parameters)         rescor[i]<-cor(x[,1],x[,2], method = "pearson")         resvar[i]<-var(x[,1])         a<-lm(x[,1]~x[,2])         pvals[i]<-lmp(a)         if (pvals[i] < 0.05){             numsig = numsig + 1         }         if (pvals[i] > 0.05){             numnotsig = numnotsig + 1         }         }          Strain Plate Growth       1      1     200        1      2     210       1      3     190       2      1     150        2      2     130       2      3     140       ...            Strain Plant Growth       1      1     70       1      2     40       1      3     50       1      4     45       2      1     80       2      2     90       2      3     85       2      4     75       ...


title: Is there a better term for frequency distribution of a term frequency?

text: I am wondering if there is a better way of formulating this and also if this is a correct way of formulating it. What I'm trying to represent in the frequency distribution is the count of documents that contain a certain count of the term I'm interested in. For example, the document sample contains 50 documents with 1 occurrence of the term, 30 documents with 2 occurrences, 10 with 3, and so on. I'm looking to answer three questions: 1. is _frequency distribution of a term frequency_ a correct way to formulate what I described, 2. if yes, is it clear what I mean, and 3. is there a better way to formulate it?


title: Why prediction of a predicted variable from a discriminant analysis is imperfect

text: I am puzzled by something I found using Linear Discriminant Analysis. Here is the problem - I first ran the Discriminant analysis using 20 or so independent variables to predict 5 segments. Among the outputs, I asked for the Predicted Segments, which are the same as the original segments for around 80% of the cases. Then I ran again the Discriminant Analysis with the same independent variables, but now trying to predict the Predicted Segments. I was expecting I would get 100% of correct classification rate, but that did not happen and I am not sure why. It seems to me that if the Discriminant Analysis cannot predict with 100% accuracy it own predicted segments then somehow it is not a optimum procedure since a rule exist that will get 100% accuracy. I am missing something? Note - This situation seems to be similar to that in Linear Regression Analysis. If you fit the model $y = a + bX + \text{error}$ and use the estimated equation with the same data you will get $\hat{y}$ [$= \hat{a} + \hat{b}X$]. Now if you estimate the model $\hat{y} = \hat{a} + \hat{b}X + \text{error}$, you will find the same $\hat{a}$ and $\hat{b}$ as before, no error, and R2 = 100% (perfect fit). I though this would also happen with Linear Discriminant Analysis, but it does not. Note 2 - I run this test with Discriminant Analysis in SPSS.


title: Unbalanced design effect

text: In my project, I am looking at the relationship between X and Y by the grouping factor "gender". This is at two stages. The first is to see how XY correlate in each subgroup (Pearsons r) and the second is to see if there is a gender difference in the XY correlation (general linear model). I am using SPSS and the scatter plot shows a linear relationship. My sample comprises 60 males and 95 females. Should I be correcting the effect of unbalanced design (i.e. different numbers in each subgroup) for my analyses? I know Sum of Squares Type III corrects it for GLM but I am unsure about the bivariate correlation. (Does it matter for a correlation analysis?) **EDIT** In response to the comments below:   1. I have used Probability Proportional to Size (PPS) cluster sampling and have taken into consideration the population distribution. For example, if the chosen entity has 700 members and half of them are females, I have strived to have this reflected in my sample.   2. I have inadvertently simplified the question above. I am testing two hypotheses. The first is to ascertain the relationship between X and Y. It is a simple bivariate correlation analysis. (Please ignore the first stage mentioned in the question above!). The second is to ascertain the effect of gender on the X-Y relationship.   3. I have created the interaction term and I have included the main effects in my model. My question again:   1. Should I 'correct for' (unsure if this is the right terminology) for the different subgroup sizes in the X-Y correlation analysis? This is a simple bivariate correlation analysis. (i.e it has nothing to do with subgroup analysis by gender)   2. Should I correct for the different subgroup sample sizes when doing the subgroup analysis by gender?


title: Finding the standard error for an estimate of the change in mean z-score

text: I have 2 normally distributed, standardized (mean = 0; SD = 1) measures of an outcome taken at 2 time points. To estimate the standard error for the mean change between the 2 timepoints, a colleague has suggested the following: If the data are collected in a different sample at each timepoint, the variance of the mean change in standard deviation score (SDS) is 4/n, where n = sample size, and n/2 measurements are taken at each time point. If the data are collected in the same sample at each time point, the variance of the mean change in SDS is ((4 * (1 - r)) / n), where r is the correlation between the 2 measures, and where n = sample size, and n/2 measurements are taken at each time point. I can't get the calculation for the longitindal case to match their results when I try to replicate them; I can't get in touch with the person who has suggested these; and I'm not clever enought to figure out of they are simply incorrect, or if I am making a mistake somewhere. Edit: I can replicate the standard errors from my colleague with the following, but am not convinced this is correct, particuarly since the variance of the average of two standardized normal variables with correlation r should be (1/2) + (1/2 * r). Longitudinal data: sqrt(4 * (1 - r)) / 2 * sqrt(n/2) Cross-sectional data:   sqrt(4) / sqrt(n)


title: Heuristics for optimizing ν-SVM?

text: Do you know any good heuristics for finding optimal value of ν in case of ν-SVM classification? In this particular problem I have a radial basis kernel, if it helps.


title: Smallest Spectral Norm / Deviation Inequality

text: Consider $A_{m \times n}$ be an i.i.d. random matrix with finite first to fourth moments. There is a good number of asymptotic and non-asymptotic results regarding the spectral norm of $A$, $\|A\|_2$, e.g., [1, Proposition 2.4] yields a deviation result on the **largest** spectral norm, and also reported in [1, Theorem 2.5] there is a very nice result by R. Latala on the expectation $\mathbf{E}[\|A\|_2]$ under suitable assumptions. I would like to ask whether anyone knows of a result on the **smallest** spectral norm, so that I could assume with very high probability that $\|A\|_2 \geq \tau$ where e.g. $\tau = \mathbf{E}[\|A\|_2] - t$ for some $t > 0$. Looking at the proof in [1, Proposition 2.4] it seems to be non-trivial to find a value for which the previous holds with very high probability, at least for someone who is not familiar with this type of argument. Is there any way to prove an inequality of the type: $$P[\|A\|_2 < \mathbf{E}[\|A\|_2] - t] > 1-e^{-c t^2}$$ for some random matrix so that $\min_A \|A\|_2 \approx \mathbf{E}[\|A\|_2] - t$ for some $t > 0$, either for finite or infinite dimensions provided that $m/n \rightarrow \gamma$? [1] Rudelson, Mark, and Roman Vershynin. "Non-asymptotic theory of random matrices: extreme singular values." arXiv preprint arXiv:1003.2990 (2010).


title: When should one avoid using Durbin-Watson test for autocorrelation?

text: Over 60 years ago Durbin and Watson suggested a testing procedure for assessing autocorrelation in regression relationships.The test is known to not work in the presence of lagged dependent variables, in which case Breusch- Godfrey test from the late 1970s applies. Any other known situations where the Durbin-Watson test should not be used?


title: Benefits of using genetic algorithm

text: Can anyone explain to me the benefits of the genetic algorithm compared to other traditional search and optimization methods?


title: Ordinary Least Squares method: why are my regression results insignificant?

text: I have a problem in my thesis results of OLS regression being insignificant. I have 3 sectors and each sector has 130 observations. Is this sample size is sufficient or not ? Can anyone suggest some reason why my results are insignificant?


title: How to add prediction intervals to describe the variability in monthly samples of data?

text: The data analyzed here is a sample of individuals collected on a monthly basis. What would be the best way to compute "prediction intervals" for the monthly sample means, in order to indicate the natural variability that can be expected from this value month after month? Here is an example:   * All individuals are sampled with replacement from the same population whose mean is expected be the same across all months;   * $y_i$ for individual $i$ is the variable of interest, with $ y_i \sim N(y,\sigma)$, both $y$ and $\sigma$ are unknown;   * $N_m$ is the sample size for month $m$;   * $\bar y_m$ is the sample mean for month $m$. **What would be the best way to compute a 95% prediction interval for $\bar y_m$?** I thought of doing the following:   * $\bar y $ is the sample mean across the whole sample;   * $\bar \sigma $ is the standard deviation across the whole sample;   * For each month $j$, I compute the interval as $[ \bar y - 1.96 \frac{\bar \sigma}{ \sqrt{N_m}} ; \bar y + 1.96 \frac{\bar \sigma}{ \sqrt{N_m}} ]$. I'm afraid this is not statistically correct. Simulations I've been doing show that this interval contains the monthly value more than 95 times out of 100. I assume this is caused by the fact that I do not consider the uncertainty in the estimation of $\bar y$. **What would be the correct way to proceed?** Thanks.


title: What test can I use to compare slopes from two or more regression models?

text: I would like to test the difference in response of two variables to one predictor. Here is a minimal reproducible example.               library(nlme)      ## gls is used in the application; lm would suffice for this example     m.set <- gls(Sepal.Length ~ Petal.Width, data = iris,                     subset = Species == "setosa")     m.vir <- gls(Sepal.Length ~ Petal.Width, data = iris,                     subset = Species == "virginica")     m.ver <- gls(Sepal.Length ~ Petal.Width, data = iris,                     subset = Species == "versicolor")      I can see that the slope coefficients are different:               m.set$coefficients     (Intercept) Petal.Width        4.7771775   0.9301727     m.vir$coefficients     (Intercept) Petal.Width        5.2694172   0.6508306      m.ver$coefficients     (Intercept) Petal.Width         4.044640    1.426365       I have three questions:   1. How can I test the difference between slopes?   2. How can I test the difference between residual variances?   3. What is a simple, effective way to present these comparisons?  A related question, Method to compare variable coefficient in two regression models, suggests re-running the model with a dummy variable to differentiate the slopes, are there options that would allow the use of independent data sets?


title: How to generating MAR data with fixed proportion of missing?

text: To generate data with missing at random (MAR) mechanism, usually we can first generate a complete data set then model the missing probability for variable Y, i.e. $Pr(Y=missing|{\bf X})$, using logistical model, i.e. $Pr(Y_i=1|{\bf X_i})=\frac{exp(\beta'{\bf X_i})}{1+exp(\beta'{\bf X_i})}$. Based on this probability we can generate the binary variable indicating if $Y$ is missing: $I_i$ = $rbinom(1,p_i)$. If $I_i=1$ then we delete corresponding value of $Y$. My questions, if we would like to control the proportion of missing of $Y$, how can we do that? Since by using above method the proportion of missing depends the value of $X's$ then the generation of indicator variable is random so it's hard to control the overall proportion of missing.


title: dummy variables with overlapping categories?

text: I'm quite familiar with traditional dummy variable coding - code 1 for presence of the attribute and 0 for absence. A multi-category variable is then represented by a series of dummy variables while omitting 1 category as the reference, so for a variable with n categories I would include n-1 dummy variables. Simple. But what happens if I have overlap in my categories? Here's a simple (slightly contrived) example to illustrate. Let's say I'm looking at the effect of different sports on injury (a dichotomous outcome). There are 6 sports - football, baseball, basketball, soccer, lacrosse, and hockey. Now, I know what you're thinking, these ARE mutually exclusive, there is no overlap. True, I could represent these sports with 5 dummy variables and use one, say football, as the reference. But instead I want to look at some facet of the sport that is related to injury. Put differently, it's not the 'sport' per se, but the actions involved in playing each sport. Some sports involve the same actions, so there is overlap. I would like to have dummies such as the following:   1. 'ball' is hard (baseball, hockey, lacrosse)   2. all players wear helmets (lacrosse, hockey, football)   3. floor/ground is hard (hockey, basketball) Now, I think I can do this so long as the dummy variables are not highly collinear. That would be tantamount to the so called 'dummy variable trap'. Right? How would I check this? VIFs for the dummies? Condition number? Is there anything else I need to look out for? Anything I'm missing? In my actual application I'm thinking of around 5 'facets' and there are well over 50 different categories. I can collapse these categories down into 5 or so catchall categories but I'd rather not do that for theoretical reasons that we don't need to get into at this point. I could let the machine chose the 'dimensions' or 'facets' via exploratory factor analysis, but I have a very specific set of theoretical 'facets' that I wish to test, hence the preference for dummy variables of my choosing.


title: Confidence intervals for two-part regression model

text: I am working with a two part regression model for semi-continuous data slightly modified from Duan et al. The Duan et al. model is used to predict medical expenses over the course of a year. There is a fraction of people who have zero expenses, and the rest have positive expenses that are log-normal distributed. The Equations From Duan _et al._ $I_i = x_i\delta_1+\eta_{1i}$,       $\eta_{1i} \sim N(0,1)$ Where $\mbox{MED} > 0$ if $I \ge 0$, and $\mbox{MED} = 0$ otherwise. The second equation is a linear model on the log scale for positive expenses: $\mbox{log}(\mbox{MED}_i|I_i > 0) = x_i\delta_2 +\eta_{2i}$, $\eta_{2i} \sim N(0,\sigma^2)$ I'm using a similar model but allowing for negative values, so use a normal distribution for the values rather than a log-normal. I'm fitting the model using the BayesGLM package in R with normal priors on the regression coefficients, which is equivalent to L2/Ridge Regression. The model has two equations the first is probit for whether or not the event is zero, the other is a linear model. I would like to be able to find the mean effect, and confidence intervals around the mean. I know that for each part of the model, I can calculate this but I'm not sure how to calculate it for the combined model. For the linear model I can look at the estimate of the effect and the standard error to get the confidence intervals, but how do I layer in the probit model? For example if a variable has a positive effect in the linear componenent and a negative effect in the probit component, how can I figure out if the effect is overall positive or not? > Duan, N., Jr, W. M., & Morris, C. (1983). A comparison of alternative models > for the demand for medical care. _Journal of Business & Economic Statistics_


title: Poisson Rate Regression: Offset?

text: I'm comparing two groups on the # of hospitalizations they have in a 5 year period. Length of hospitalization stay is highly variable. For example, a patient might have 3 hospitalizations that in total lasted 100 days, where as another patient might have had 50 visits, each lasting a single day. Therefore, It would seem natural to factor in length of stay (in days) into the picture. Right now I'm thinking of a Poisson Rate Regression, with a single categorical predictor (Group) and # of hospitalizations as the DV. I'd like to include Total Length of hospitalization stays as an offset, which is in days. My question is, I have some patients who did not have any hospitalizations, therefore their length of stay is effectively 0. In order to hit the Poisson Rate Regression model (in SAS anyway), I must log the offset term. I clearly can't take the log of 0, so I'm wondering what can I do in this situation? Can we add a constant of 1 to both Hospitalization and Total Length of Hospitalization to remedy this? Any input is appreciated, Best (cross posted @ Talkstats.com)


title: Mean square error or mean squared error

text: As a non-native english speaker I was wondering which of the **square** or **squared** expression I should use. For instance in mean **square** error or mean **squared** error. According to the Internet, it seems both forms are used indistinctly. Is one expression more square than the other ?


title: Can I combine Standard errors of coefficients with an unbalanced data set?

text: I have just read the answers to the question ...... Adding coefficients to obtain interaction effects - what to do with SEs? .......and found this really helpful! I am looking to do something similar in R, but I have slightly different sample sizes between groups. Can I still use the same equation to combine the two errors to give me the new Std. error? Thank you in advance for any help Crystal


title: Mixture model for dependent observations with additive group-level confounders

text: I'm looking for a special type of mixture model (described below) and I'm hoping to get some hints with regards to relevant literature to look at or names to be searching for. On the general level, I'd like to do clustering using a mixture of Gaussians to infer the class label for each observation (cell types in my case). The complication I'm facing is that some of the individual observations are dependent because they are taken from the same subject. I'd like to treat the subjects as a random effect that adds an additive offset to each observation. Thus, the generative model looks something like $x_{ijk} = \mu_i + \nu_j + \epsilon_{ijk}$ where $\mu_i$ is the mean for each cell type, $\nu_j$ is the random offset for each subject, and $\epsilon_{ijk} \sim N(0, C_i)$ is multivariate normal noise. The class means $\mu_i$ and noise covariance matrices $C_i$ are model parameters, while the additive offset for each subject $\nu_j$ and the class memberships are latent variables. We can assume $\nu_j$ to be Gaussian with a low-rank covariance. Any hints?


title: How to implement this classification/labelling problem?

text: I have a dataset that has 12 explanatory variables for every 1 observation. I hypothesize that the data is generated by underlying process which undergoes several different phases, and that the 12 explanatory variables would somehow help in identifying clusters. I want to write a machine learning algorithm that can:   1. At first pass, identify the clusters within the data   2. The data is dynamic in that new (additional) data is generated periodically. After classification (step 1 above), I want to be able to correctly label (i.e. classify) any new data not previously seen into one of the previously identified classes/clusters (or fail gracefully). I assume that Y(c) ~ X(c) + error where:   * Y(c) is a nx1 vector of observations belonging to cluster C   * X(c) is a nx12 vector of explanatory factors that 'belong' to cluster C Observations in different clusters will differ from each other by having different distribution shapes. That is to say observations WITHIN a cluster will have a different distribution shape compared to observations from another cluster. I am relatively new to machine learning, and would like some guidance on how to implement such an algorithm (or perhaps one already exists?) I would be particularly interested in finding out how to 'classify' observations based on determining the _shape_ of the empirical distribution of the observed data. Last but not the least, I would also appreciate some advice on whether to implement this in R, or Octave (and why).


title: What is the meaning of this formula in R?

text: I have rows of data with columns `age`, `sex`, `education` and `income`. I am doing homework that asks me to predict income with naive Bayes in R with the formula   `income ~ age + sex + educ` I know formulas tell R the shape of the model to apply in the method. I also know that formulas don't refer to any specific variables in a dataset -- they only give the shape of the data. So with all that said, what does the formula `income ~ age + sex + educ` mean? I am guessing that it means               p(income|age, sex, income) = P(y|x1,x2,x3) = P(x1|y) * P(x2|y) * P(x3|y) * P(y)


title: Compare group difference with replicates?

text: I have a measurement for two groups of subjects where each subject contains 3 replicates. For example, the data is like below: > Group Subject Replicate Value   >  1 A 1 0.2   >  1 A 2 0.3   >  1 A 3 0.25   >  1 B 1 0.4   >  1 B 2 0.3   >  1 B 3 0.7   >  2 C 1 0.2   >  2 C 2 0.1   >  2 C 3 0.3   >  2 D 1 0.4   >  2 D 2 0.2   >  2 D 3 0.25   > I am thinking of using ANOVA where the subject is the error term,e.g. aov(Value~Group+Error(Subject)) in R. Is it correct? Thanks! Gim


title: Asymptotic normal distribution via the central limit theorem

text: I have a sample $n = 100$ with two "successes" (Two kids having a disease among 100). So we obviously have a binomial distribution. First I had to compute the maximum likelihood (ML) estimator $\hat{p}$. I got $\hat{p} = \frac{k}{n}$. Now, I have to derive asymptotic normal distribution for $\hat{p}$ via the central limit theorem (CLT). I know that the expected value of $\hat{p}$ is not infinite and also variance is not infinite, so I know it will be normally distributed. I have to know expected value and variance of $\hat{p}$ to get the asymptotic normal distribution, right? I know that expected value is $\frac{k}{n}$. But what is variance?


title: Why is the variance of $X-Y$ equal to the sum of the variances when $X,Y$ are independent?

text: ![enter image description here](http://i.stack.imgur.com/nXT1s.png) I have one question about this. I know that if we have $\mathrm{X}_1,\mathrm{X}_2,\ldots,\mathrm{X}_n$ independent and normally distributed random variables, then the sum $\mathrm{X}_1+\mathrm{X}_2+\ldots+\mathrm{X}_n$ has the normal distribution with mean $M_1+M_2+..+M_n$ and variance $\sigma^2_1 + \ldots + \sigma^2_n$. Why is in this problem the difference $W-M$ the mean obtained by subtraction and variance obtained by addition? Thank you.


title: Methods to find the relationships between independent and dependent variables instead of regression

text: Which statistical method can be used to find the relationships between independent and dependent variables instead of regression? And their advantages over regression. It should be mentioned I am not able to do regression due to less number of observations that I have 5 independent variables 3 categorical and two continuous.


title: How to compare AUCs using the boot package in R?

text: I would like to compare AUCs for two treatments using the boot package in R. However, something wrong with my R code and I don't know how to fix it. Can someone please give me a hint? many thanks in advance               require(boot)     #Create Data     set.seed(123)     N = 150     y1 = rnorm(150)     y2 = rnorm(150,1,0.2)     y3 = rnorm(150,2,0.5)     time = rep(1:3, each = 50)     trt = rep(rep(0:1, each = 25), times=3)     Data = data.frame(trt,time,y1,y2,y3)          #Trapezoid Rule     trap.rule = function(x,y) sum(diff(x)*(y[-1]+y[-length(y)]))/2          # Sum y1-y3 for each trt at each time point     SUM = function(data){          set1 = subset(x = data, data$time==1)          set2 = subset(x = data, data$time==2)          set3 = subset(x = data, data$time==3)          s1 = by(set1[,c("y1","y2","y3")], set1$trt, sum)          s2 = by(set1[,c("y1","y2","y3")], set2$trt, sum)          s3 = by(set3[,c("y1","y2","y3")], set3$trt, sum)          result = matrix(data = c(s1[1],s2[1],s3[1],                                   s1[2],s2[2],s3[2]), ncol=2)          colnames(result) = c("trt = 0", "trt = 1")          output = cbind(time = 1:3,result)          return(output)     }          # Statistic for boot     AREA = function(data,i){             d = data[i,]             MQ = SUM(data = d)             A.0 = trap.rule(x = 1:3, y = MQ[,2])             A.1 = trap.rule(x = 1:3, y = MQ[,3])             area = A.1 - A.0             result = c(A.0, A.1, area)             return(result)     }          area.boot = boot(data = Data, statistic = AREA, R = 1000)     area.boot      However, there is an error message               Error in tapply(1L:49L, list(`set2$trt` = c(0L, 1L, 1L, 1L, 0L, 1L, 0L,  :        arguments must have same length


title: How to sample from this dirichlet distribution with an L1 prior?

text: I'd like to draw a sample from a distribution with p.d.f $$f(p,q,r,s) \sim \mathrm{e}^{-w(|p+q-r-s|+|p-q-r+s|+|p-q+r-s|)}p^aq^br^cs^d \mathbb{1}_{p+q+r+s=1}$$ $w > 0$ is a free parameter (which could be very large or very small) and $a,b,c,d$ are all non negative integers. **What I'm currently thinking**   1. rejection sampling from the Dirichlet distribution      * too slow if w is large   2. drawing a sample from a Dirichlet distribution use the likelihood weighted sample to infer better parameters from the D.d. Use as a proposal for MCMC.     * a bit Voodo, hard to check for convergence.   3. if I separate the 8 quadrants, the absolute value disappear and the integral becomes somewhat tractable.     * I suppose I can then partition the simplex where my distribution lives and if I can calculate integrals I can randomly zoom in on a region. It doesn't look trivial.   4. somehow figure out the marginal distribution of pairs and use Gibbs sampling.     * Ah but the L1 constraint means that everything can be pretty correlated and Gibbs sampling may utterly fail.   5. reparametrize with $x=p-q-r-s,~ y=p-q-r+s,~ z=p-q+r-s$     * Ok great, but now the correlation comes from the Dirichlet part of my distribution. Gibbs sampling might still be a problem. **What I'm looking for would be**   * fast   * unbiased (or at least converging to something unbiased fast enough with means to check convergence)


title: ARIMA Representation using polynomial long division

text: I have 2 questions from Tsay's book, dealing with ARIMA representation. I. How does the division of the polynomial work - what are the steps that lead to the solution they show? They show a "simple" example but I am unsure how this math works?: ![enter image description here](http://i.stack.imgur.com/iQ8UA.jpg) II. They also show the following, can anyone explain what theta(1) in the denominator of the first equation represents? ![enter image description here](http://i.stack.imgur.com/RgSeW.jpg)


title: How to combine two fitted variables

text: I have fitted two (and probably more) sets of data to a model using least- squares regression. So I've got the parameter values together with their errors and a reduced chi-squared for both. The data is 2D surface brightness versus radius. How do I go about combining these parameters and chi-squareds to give me an average answer with its own error?


title: In linear regression, how to test if doubling of x leads to a doubling of log-transformed y?

text: In a simple linear regression with a log-transformed y variable, how do I test if a doubling of the x variable leads to a doubling of the y variable?


title: Compute a probability using stochastic simulation

text: My problem is to solve the following optimisation problem using GA (Genetic Algorithm)and stochastic simulation. The goal is to solve the maximisation problem : \begin{equation*} \begin{aligned} \text{Max}_{x_{1},\cdots,x_{n}} \qquad &\mathbb{E}[\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}]\\\ \text{s.t.} \qquad &\left\\{ \begin{array}{rcr} & \mathbb{P}r\\{(0.4-(\xi_{1}x_{1}+\xi_{2}x_{2}+\cdots+\xi_{n}x_{n}))\geqslant r\\}\leqslant \alpha(r)\\\ & \forall r\geqslant 0\\\ &x_{1}+x_{2}+\cdots+x_{n}=1\\\ &x_{1} \qquad \forall i \in \\{1,2,\cdots,n\\}.\\\ \end{array} \right. \end{aligned} \end{equation*} where $ n=5 $, $ \ \xi_{i}\sim N(\mu_{i},\sigma_{i}) $ with $\mu_{i}\in\\{0.1, 0.2, 0.3,0.4,0.5\\}$, $\sigma_{i}=1$ and $ \ \alpha(r)=\dfrac{1}{(r+1.1)^{4}} $, $\ r\geqslant 0 $.


title: Mixture of beta distributions: full example

text: I have a vector of measured proportions and want to   1. test if these proportions follow a bimodal distribution    2. characterize the two underlying distributions   3. determine for each data point which distribution it belongs to   4. overaly the density distributions on the histogram I've been using "betamix" from the betareg package but haven't figured out steps 2 & 3\. I've searched stackoverflow for solutions but haven't found a clear answer. Here's my current code:               # parameters of distribution #1      alpha1 <- 10     beta1 <- 30          # parameters of distribution #2      alpha2 <- 30     beta2 <- 10          # Generate bimodal data     set.seed(0)     d <- data.frame(y = c(rbeta(100, alpha1, beta1), rbeta(50, alpha2, beta2)))     # correction recommended in cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf, cf first paragraph in section 2, page 3):     n <- length(d$y)         d$yc <- (d$y* (n-1)+0.5)/n          # histogram      par(mfrow=c(1,2))     hist(d$yc, 50)          # fitting mixtures of beta distributions     uni.modal <- betamix(yc ~ 1 | 1, data = d, k = 1)     bi.modal <- betamix(yc ~ 1 | 1, data = d, k = 2)          # 1) test for bimodality:     lrtest(uni.modal, bi.modal)          # 3) determine for each data point which distribution it belongs     d$group <- (posterior(bi.modal)[,1] <= posterior(bi.modal)[,2])+1         plot(d$y, col=d$group)               # 2) characterize the two underlying distributions     # betamix uses a different parametrization than dbeta (cf cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf).      # instead of alpha and beta, betareg parametrizes the beta distribution using mu and phi, where     # mu=alpha/(alpha+beta)     # phi= alpha + beta          # ERROR: converting mu and phi to alpha and beta     mu1 <- coef(bi.modal)[1,1]     phi1 <- coef(bi.modal)[1,2]          (a1 <- mu1*phi1) # output: 4.61028     (b1 <- (1-mu1)*phi1) # ouput:-0.7240308               # 4) overaly the density distributions on the histogram      ![enter image description here](http://i.stack.imgur.com/KbAQM.png)


title: What should I use to run a pre and post after intervention with a control group 15 participants?

text: I want to know the effect of reflective video analysis on teacher intentions to stay or leave the teaching position, teacher efficacy, and teacher job satisfaction in rural middle schools. I will have a control group (without random assignment) and a non control group. Both groups will complete a pre and post likert scale survey measuring teachers intention to stay or leave the position, teacher efficacy,and teacher job satisfaction.


title: Asynchronous (irregular) Time Series Analysis

text: I am trying to analyze the lead-lag between time series of two stock prices. In regular time series analysis, we can do Cross Correlaton, VECM (Granger Causality). However how does one handle the same in irregularly spaced time series. The hypothesis is that one of the instruments leads the other. I have data for both symbols to the microseconds. I have looked at RTAQ package and also tried applying VECM. RTAQ is more on a univariate time series while VECM is not significant on these timescales.               > dput(STOCKS[,]))     structure(c(29979, 29980, 29980, 29980, 29981, 29981, 29991,      29992, 29993, 29991, 29990, 29992), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("Pair_Bid", "Calc_Bid" )), index = structure(c(1340686178.55163, 1340686181.40801, 1340686187.2642,      1340686187.52668, 1340686187.78777, 1340686189.36693), class = c("POSIXct", "POSIXt"), tzone = ""), class = "zoo")


title: Estimation of parameters as a mode of posterior distribution

text: I was referring to this http://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation. Here it is mentioned that estimating the parameters $\theta$ is actually finding the mode of the posterior distribution. I didn't get how this is true.


title: plot data on a map

text: I’ve a data set contains three variables: longitude, latitude, population (in USA). I need to plot this data on a map. I use the following R code (it works fine): map("state", xlim = range(lon), ylim = range(lat)) text(lon, lat, population) box() But the thing I’m facing is that the map looks really ugly, and It’s hard to read since data displayed on the map are mixed with each other. I was wondering if there is a better way to plot these data on a map. Thanks.


title: Normalize sequence lengths

text: I'm implementing a Markov Chain Model for web site access sequences. Some sequences are very lengthy and some are very short containing a single transition. Lengthy sequences generate a very low probability while shorter ones are having a high probability. Is it possible to normalize sequence lengths? Please let me know any idea, tutorial, conference article etc. Thanks.


title: Fitting (multlple) linear models by group in R

text: I'm trying to (somewhat) elegantly fit 3 models (linear, exponential and quadratic) to a dataset with classes/factors and save p-values and R2 for each model and class/factor. Simple dataset with 3 variables: x,y, and class. What I can't figure out is how to force each of the 3 models to fit to each of the 3 classes. What I have now fits each model to the complete dataset. The next question is how I then output p-values & R2 to a table, for each model+class My code looks like:               set.seed(100)     library(plyr)       #create datast     nit <- within(data.frame(x = 3:32),                      {                        class <- rep(1:3, each = 10)                        y <- 0.5 * x* (1:10) + rnorm(30)                        class <- factor(class)                 # convert to a factor                      }                      )     x2<-nit$x*nit$x #for quadratic model     forms<- paste(c("y ~ x", "y ~ x+x2", "log(y) ~ x"), sep = "") # create 3 models     names(forms) <- paste("Model", LETTERS[1:length(forms)])     models <- llply(forms, lm, data = nit)     models                # shows coefficients for each of the 3 models


title: histograms/ kernel density estimation for proportion data

text: I am looking at digital advertising data with the following columns: site_id | impressions_count| clicks_count. I am interested in a histogram of the click through rate (clicks/impressions) across my population of sites : to answer the question on 50% of the sites the CTR is greater than x% etc. Now the sampling error of CTR on each site will depend on the number of samples ( =impressions) and population site CTR. Is there any way of taking this into account? I could see three ways a) try to use the site binomial standard error in gaussian kernel (or use beta distribution?) \- this seems to confuse uncertainty with spread of values b) calculate upper and lower confidence intervals of the CTR for each site, and separately make histograms of these. c) use glen_bs answer in Histogram of uncertain data to add error bars to the histogram of site CTRs. Any pointers in the literature would be welcome- it would seem a standard problem.


title: How to find which time series is trending more?

text: Let us say I have two sets of time varying series as shown below:               a <- c(895,0,0,0,0,832,0,3084,0,434,0,0,0,0,853,0,3727,0,513,0,0,0,0,1248,0,3704,0,646,0,0,0,0,2066,0,5500,0,424,0,0,0,0,1069,0)     b <- c(234,0,0,0,0,521,0,1683,0,152,0,0,0,0,740,0,3242,0,355,0,0,0,0,1117,0,3047,0,443,0,0,0,0,832,0,3213,0,326,0,0,0,0,1235,0)      Here, each value represents the number of times a person has listened to a song. Also, the data is daily which means if the first value (a[1] = 895) is for 10th Feb then the last value will be for 25th Mar. I want to find out whether a or b is trending more and by what margin, so if I have let say another series c then I can arrange them in the trending order for example c > a > b . Is there a statistical approach to do this.


title: Error bars using median absolute deviation

text: I have tried to find a solution to my question on Google, but I can’t seem to find much information about error bars and median absolute deviation and I do not know much about statistical error analysis so any help would be greatly appreciated. I am creating a semi-log plot for my astronomy research that splits the data into 5 equally spaced bins in log base 10 (x-axis) and plots the average value for each bin (y-axis). Thus each bin has a larger sample size as the x-axis increases. Since it’s a statistical study there is a lot of potential for uncertainty and there is a high presence of outliers (non-Gaussian error distributions). As a result, I would like to use median absolute deviation (MAD) error bars since MAD is less susceptible to outliers. So my question is: Much like standard deviation and standard error, where standard error is just stdev/sqrt(N) and N is the size of the sample – is there an analog to this for median absolute deviation? The error bars I get when I use MAD/sqrt(N) look correct, but I am not confident in explaining my reasoning for using this approach. Is dividing MAD by the square root of the sample size an acceptable means to produce error bars? And if not, do you have any other suggestions? Also, from research I’ve done it seems that median absolute deviation is a better estimator than mean absolute deviation…would you agree with this? UPDATE: Here is my original plot with error bars from median absolute deviation/sqrt(N) (NOTE: title should read median, not mean!) The first bin contains 39 samples, second contains 146, third 454, fourth 1287, and fifth 2371 samples. It looks nice, but the method for producing error bars does not seem very accurate. ![enter image description here](http://i.stack.imgur.com/mg7fx.png) I've created error bars using the bootstrap method as suggested by @Glen_b (attached below). It looks good to me, but my advisor feels that the error bars have been overestimated in this method. ![enter image description here](http://i.stack.imgur.com/bmW6F.png)


title: Finding statistical significance between two binary tests

text: I have two sets of data. One with a sample size of 82 with 53 "hits" and 29 "Misses." And a second sample of 105 with 67 "hits" and 38 "misses" Given that the second set of data is a control, is there a way to show that the results of the first data set are or are not statistically significant? Thanks in advance, It's been a while since I took a stats course and googling for the answer left me confused.


title: Is the Pearson goodness-of-fit test parametric or non-parametric?

text: In this question there was an interesting discussion concerning the Pearson goodness-of-fit test going on which was far from conclusive: Is there any statistical test that is parametric and non-parametric? Also, when you surf around you find arguments for and against its being of either sort - obviously because the line between parametric and non-parametric statistics seems to be rather blurred?!? But I ask this question because many people still differentiate in this way: Is the Pearson goodness-of-fit test parametric or non-parametric (or even some kind of combination)? (NB: There is a similar but more general question - which was _not_ answered here: Are Goodness-of-fit tests parametric or nonparametric?)


title: What algorithms are available for logistic regression?

text: I am trying to implement a logistic regression function in c++, and not sure what algorithm to use. So far I have heard of these:   1. Newton-Raphson   2. IRLS   3. Gradient descent Are there other algorithms available? What are their pros and cons? Is there an algorithm generally recognized as superior to the others?


title: The proper scaling for generating uniformly distributed points in the n-d ball of *non-unit* radius?

text: The question of generating uniformly distributed points on the surface of n-dimensional unit ball has been already posted here a dozen of times. What I'm interested in is the proper scaling when we want to generalize the method to uniform distribution on the surface of _non-unit_ radius n-ball. Specifically, I'm looking for the scaling for the "projecting Gaussian on the ball" method. I know this is a fairly easy question, but I cannot figure out the math and I appreciate any comments or references.


title: Weight Size in Neural Networks

text: I am using a Neuronal Network in combination with Reinforcement Learning. The network should learn the values of three actions in given states. The reward from the environment is scaled to [-0.9,0.9]. The network consists of one input layer with 30 nodes of which not more than three are activated at the same time. The hidden layer consists of 30 nodes as well and tanh is used as activation function. The output layer consists of said three nodes, which also are activated using tanh. When training my network, it often finds solutions for the given task, but often enough the values returned by the network are getting closer and closer to -1 and 1, although the reward is normalized to [-0.9,0.9]. Through debugging I found out that the weight from the hidden to the output layer are getting bigger and bigger. Since the weight size is used for the backpropagation, a huge error is passed through the network, causing even bigger weights and so on. Is there a way to prevent the network from getting this vicious circle?


title: Non-parametric paired analysis of variance with unequal data

text: I have tried to search, but couldn't find an answer for my problem. I have data from different subjects that are exposed to different conditions. But they are exposed to the different conditions unequal times (sometimes even zero times). Like this                          Condition 1               Condition 2                      Condition 3             -------------------------------------------------------------------------------------------          Subj 1 |   Meas1 Meas2 Meas3         Meas1 Meas2 Meas3 Meas4 Meas5    Meas1 Meas2 Meas3          Subj 2 |   Meas1                     Meas1 Meas2 Meas3                Meas1           Subj 3 |   Meas1 Meas2 Meas3 Meas4   Meas1                            Meas1 Meas2 Meas3          Subj 4 |    ---                      Meas1 Meas2 Meas3 Meas4          Meas1           Subj 5 |   Meas1 Meas2 Meas3         Meas1 Meas2 Meas3                Meas1 Meas2           Etc...      I want to test (non-parametrically), if the measurements under the conditions are different, like in an Kruskal-Wallis or Friedman test, but I want to keep the information from the pairing by subjects. And I don't want to have to discard measurements just because there were no measurements for one of the conditions. What kind of test am I looking for? Also, if it can be done in Matlab, I would be really happy!


title: p-value as a distance?

text: Can p-values between multiple pair-wise tests be considered as a similarity/distance measure and multidimensional scaling be applied over a pair-wise matrix of p-values to reduce the dimensionality? This is a soft- question, but what would be the biggest issue here, and how could that be best overcome? (ex: triangular inequality?)


title: How can I update a Bayesian network model given new data on only a subset of the variables in the original model?

text: There are several methods for inferring network structure in Bayesian networks, given data. In my case I have a Bayesian network model built from old data, and I have a new source of data that I want to use to update the model, both in terms of structure and parameters. This new data source has a number of observations, much higher than the past data, making it ideal for structure inference. The problem is that the new data covers only a subset of the variables in the original model. Has anyone heard of a way to update BN model network structure (and parameters) given only a subset of the variables in the model?


title: Rugarch and rolling estimation

text: I use Rugarch for a long time in order to calibrate GARCH models on FX rates time series and perform simulations. I am trying to understand the ugarchroll method. However even if I can find plenty of examples and indications on how this method has to be used in the documentation as well as in the vignette, I still do not get what rolling estimation exactly is and what is its purpose. Could you give me some explanations ? Thanks.